text:
# Theory of the Mercury’s spin-orbit motion and analysis of its main librations

## 1 Introduction

Before 1965, the rotational motion of Mercury was assumed to be synchronous with its orbital motion. In 1965, Pettengill & Dyce discovered a 3:2 spin-orbit resonance state by using Earth-based radar observations (the Mercury’s rotation period is $`P_\varphi `$ = $`58.646`$ days while the orbital one is $`P_\lambda `$ = $`87.969`$ days). This surprising resonance results from a non-zero eccentricity and a permanent asymmetry in the equatorial plane of the planet. In addition, the 3:2 resonance strongly depends on the functional dependance of the tidal torque on the rate of the libration in longitude. Moreover the 3:2 resonance state is preserved by the tidal torque (Colombo & Shapiro 1966). The main dynamical features of Mercury have been established during the 1960’s in some pioneer works such as Colombo (1965), Colombo & Shapiro (1966), Goldreich & Peale (1966) and Peale (1969). Goldreich & Peale (1966) notably studied the probability of resonance capture and showed that the 3:2 ratio is the only possible one for a significant probability of capture. In addition, in a tidally evolved system, the spin pole is expected to be trapped in a Cassini state (Colombo 1966; Peale 1969, 1973). The orbital and rotational parameters are indeed matched in such a way that the spin pole, the orbit pole, and the solar system invariable pole remain coplanar while the spin and orbital poles precess. The reader may find in Balogh & Campieri (2002) a review report on the present knowledge of Mercury whose the interest is nowadays renewed by two upcoming missions: MESSENGER (NASA, Solomon et al. 2001) and BepiColombo (ESA, ISAS, Anselin & Scoon 2001).

Our work deals with the physical and dynamical causes that contribute to induce librations around an equilibrium state defined by a Cassini state. In order to wholly analyze the spin-orbit motion of Mercury and its rotational librations, we used a gravitational model of the solar system including the Moon’s spin-orbit motion. The framework of the model has been previously constructed by Bois, Journet & Vokrouhlický (BJV model) in accordance with the requirements of Lunar Laser Ranging (LLR thereafter) observational accuracy (see for instance a review report by Bois 2000). The approach of the model consists in integrating the N-body problem on the basis of the gravitation description given by the Einstein’s general relativity theory according to a formalism derived from the first post-Newtonian approximation level. The model is solved by modular numerical integration and controlled in function of the different physical contributions and parameters taken into account. We have extended this model to the integration of the rotational motions of the terrestrial planets (Mercury, Venus, Earth, and Mars) including their spin-orbit couplings. The updated model is then called SONYR (acronym of Spin-Orbit N-bodY Relativistic model). As a consequence, using SONYR, the N-body problem for the solar system and the spin motion of Mercury are simultaneously integrated. Consequently we may analyze and identify the different families of Hermean rotational librations with the choice of the contributions at our disposal.

Starting with the basic spin-orbit problem according to Goldreich & Peale (1966), we have computed a surface of section for the Mercury’s rotation showing its very regular behavior. We have calculated again the proper frequency for the spin-orbit resonance state of Mercury. Using our model, an important part of the present study deals with the main perturbations acting on the spin-orbit motion of Mercury such as the gravitational figure of the planet as well as the planetary effects and their hierarchy. A detailed analysis of the resulting rotational librations due to these effects is presented and described in the paper. A new determination of the Hermean mean obliquity is also proposed. Moreover, we identify in the Hermean librations the impact of the variation of the greatest principal moment of inertia on the instantaneous obliquity and on the libration in longitude. Such a signature gives noticeable constraints on the internal structure of Mercury.

## 2 Geometry of the spin-orbit coupling problem

According to Goldreich & Peale (1966), we consider the spin-orbit motion of Mercury with its spin axis normal to the orbital plane. The orbit is assumed to be fixed and unvariable (semi-major axis $`a`$ and its eccentricity $`e`$). The position of Mercury is determined by its instantaneous radius $`r`$ while its rotational orientation is specified by the angle $`\theta `$. The orbital longitude is specified by the true anomaly $`f`$ while the angle $`\theta f`$ measures the angle between the axis of least moment of inertia of Mercury and the Sun-to-Mercury line (see Fig. 1). According to these assumptions, the dynamical problem of the spin-orbit motion of Mercury is reduced to a one-dimensional pendulum-like equation as follows :

$$C\ddot{\theta }+\frac{3}{2}(BA)\frac{GM_{}}{r^3}\mathrm{sin}2(\theta f)=0$$
(1)
where $`G`$ is the gravitational constant, $`M_{}`$ the solar mass, and $`A,B`$, and $`C`$ the principal moments of inertia of Mercury. This equation of motion has only a single degree of freedom for the spin-orbit coupling, the characteristic angle of rotation $`\theta `$, but depends explicitly on time through the distance $`r`$ to the planet and the non-uniform Keplerian motion of the true anomaly $`f`$. As a consequence, it is a problem not reducible by quadrature and the equation (1) is non-integrable (Wisdom $`1987`$).

In order to know the structure of the phase space of the Mercury’s rotation, a surface of section for its spin-orbit coupling is very useful. Let be $`\alpha =\sqrt{\frac{3(BA)}{C}}`$ the asphericity of the Mercury’s dynamical figure combining the principal moments of inertia $`A,B`$, and $`C`$. The equation (1) becomes :

$$\ddot{\theta }+\frac{1}{2}n^2\left(\frac{a}{r}\right)^3\alpha ^2\mathrm{sin}2(\theta f)=0$$
(2)
where $`n`$ is the orbital mean motion, $`f`$ and $`\theta `$ as defined above. The numerical integration of the equation (2) yields to the Figure 2 showing stroboscopically, i.e. once point per orbit at the periapse, various trajectories of the phase portrait of Mercury ($`\theta `$ and $`1/nd\theta /dt`$ on the axes $`x`$ and $`y`$ respectively). Figure (2.a) illustrates the different quasiperiodic librations around the 1:2, 1:1 (synchronous), $`3`$:$`2`$, and 2:1 rotation states while Figure (2.b) shows a zoom of the present $`3`$:$`2`$ resonance state for Mercury. Let us notice that in the 1:2 rotation state, the two libration areas are shifted by $`\pi /2`$ with respect to the other rotation states. The former rotation rate of Mercury was higher and the planet could have experienced large and chaotic variations in obliquity at some time of the past (Laskar & Robutel 1993). Due to its slowing down by tidal despin, the Mercury’s rotation has then encountered different resonance states, from higher orders to the one trapped in the present state. Goldreich & Peale $`(1966)`$ have shown that this present state (the $`3`$:$`2`$ resonance) was the first one with a substantial probability of capture. This one indeed depends crucially on the functional dependence of the tidal torque acting on the spin-orbit motion. This dependence is expressed by the variation rate of the longitude libration angle $`\gamma `$ as defined below. In the end, the stabilization of Mercury at the $`3`$:$`2`$ spin-orbit resonance is due to permanent asymmetry on the equatorial plane, as well as its non-zero eccentricity equal to $`0.206`$ (Colombo 1966, Colombo & Shapiro 1966).

According to the Chirikov resonance overlap criterion (1979), the chaotic behavior appears when the asphericity of the body is larger than the following critical value :

$$\alpha ^{cr}=\frac{1}{2+\sqrt{14e}}$$
(3)
where $`e`$ is the orbital eccentricity. In the case of Mercury, $`\alpha =0.0187`$ is lower than the critical value $`\alpha ^{cr}=0.2701`$. As a consequence, the spin-orbit behavior is regular. The zoom in Figure (2. b) shows indeed that the separatrix surrounding the $`3`$:$`2`$ resonance state is very small. The width of the associated chaotic zone is then estimated to the order of $`10^{43}`$ (Wisdom et al. 1984). Consequently, tidal friction pulls Mercury accross the chaotic separatrix in a single libration period.

From the equation (1), it is possible to obtain an integrable approximated equation using the spin-orbit resonance, the spin rate $`\dot{\theta }`$ being commensurable with the mean orbital motion $`n`$. Following Murray & Dermott ($`2000`$), by introducing a new variable $`\gamma =\theta pn`$ where $`p`$ parametrizes the resonance ratio ($`p=\frac{3}{2}`$ in the case of Mercury), one may expand the equation in form-like Poisson series. Taking into account that $`\dot{\gamma }n`$, one averages all the terms over one orbital period, and finally obtain the following equation :

$$\ddot{\gamma }+\frac{3}{2}n^2\frac{BA}{C}H(p,e)\mathrm{sin}2\gamma =0$$
(4)
where $`H(p,e)`$ is a power series in eccentricity. In the case of Mercury, this expression is written as follows :

$$H(\frac{3}{2},e)=\frac{7}{2}e\frac{123}{16}e^3$$
(5)
Finally, the proper frequency of the Mercury’s spin-orbit motion is :

$$\omega _0=n\left[3\frac{BA}{C}|H(p,e)|\right]^{\frac{1}{2}}$$
(6)
which by using the values listed in Table 2 gives the proper period of $`15.830`$ years.

Balogh & Giamperi (2002) developed the equation (2) and obtained the following expression :

$$\ddot{\gamma }+\alpha _0\underset{q}{}G_{20q}(e)\mathrm{sin}[2\gamma +(1q)M]=0$$
(7)
where the $`G_{20q}`$ coefficients are eccentricity functions defined by Kaula (1966). The $`G_{201}`$ coefficient is equal to $`H(\frac{3}{2},e)`$. $`\alpha _0=\frac{3}{2}\frac{BA}{C}`$ = 1.76 $`10^4`$ for the Mercury case. Figure 3 presents then the numerical solution of the formula (7) (thanks to G. Giampieri, private communication). The angle $`\gamma `$ describes a periodic behavior with a period equal to the revolution of Mercury. The amplitude of 42 arcseconds (as) depends on the value of $`\alpha _0`$. The behavior of the $`\gamma `$ angle is nicely matched by the approximate formula of Jhen & Corral (2003) :

$$\gamma =\varphi _0\mathrm{sin}M+\varphi _1\mathrm{sin}2M$$
(8)
where $`\varphi _0=\alpha _0(G_{200}G_{202})`$ and $`\varphi _1=\frac{\alpha _0}{4}(G_{201}G_{203})`$. The authors noticed that the ratio of the two amplitudes ($`\varphi _0`$ and $`\varphi _1`$) does not depend on the $`\alpha _0`$ parameter, and $`K=\frac{\varphi _0}{\varphi _1}=9.483`$.

However, the above equations describe a spin-orbit motion of Mercury where the spin axis is normal to the orbital plane while the orbital motion is Keplerian.

## 3 The extended BJV model : SONYR

### 3.1 Theory

In order to wholly analyze the spin-orbit motion of Mercury and its rotational librations, we have enlarged a gravitational model (called BJV) of the solar system including the Moon’s spin-orbit motion. The accurate theory of the Moon’s spin-orbit motion, related to this BJV model, was constructed by Bois, Journet & Vokrouhlický in accordance with the high accuracy of the LLR observations (see previous papers : Bois et al. 1992; Bois & Journet 1993; Bois & Vokrouhlický 1995; Bois et al. 1996; Bois & Girard 1999). The approach of the BJV model consists in integrating the N-body problem (including translational and rotational motions) on the basis of the gravitation description given by the Einstein’s general relativity theory. The equations have been developped in the DSX formalism presented in a series of papers by Damour, Soffel & Xu (Damour et al. 1991, 1992, 1993, 1994). For purposes of celestial mechanics, to our knowledge, it is the most suitable formulation of the post-Newtonian (PN) theory of motion for a system of $`N`$ arbitrary extended, weakly self-graviting, rotating and deformable bodies in mutual interactions. The DSX formalism, derived from the first post-Newtonian approximation level, gives the post-Newtonian representation of the translational motions of the bodies as well as their rotational ones with respect to the locally transported frames with the bodies.

Gravitational fields of the extended bodies are parameterized in multipole moment expansions : $`(M_L^\mathrm{A},S_L^\mathrm{A})`$ define the mass and spin Blanchet-Damour multipoles characterizing the PN gravitational field of the extended bodies while $`(G_L^\mathrm{A},H_L^\mathrm{A})`$ are tidal gravitoelectric and gravitomagnetic PN fields. Because we do not dispose of dynamical equations for the quadrupole moments $`M_{ab}^\mathrm{A}`$, and although the notion of rigidity faces conceptual problems in the theory of relativity, we have adopted the ‘rigid-multipole’ model of the extended bodies as known from the Newtonian approach. Practically this is acceptable since the relativistic quadrupole contributions are very small. Consequently and because it is conventional in geodynamical research to use spherical harmonics analysis of the gravitational fields with the corresponding notion of harmonic coefficients $`(C_{lm}^\mathrm{A},S_{lm}^\mathrm{A})`$, the quadrupole moments $`M_{ab}^\mathrm{A}`$ have been expressed in those terms, according to reasons and assumptions given in Bois & Vokrouhlický (1995). Gravitational figures as well as the figure-figure interactions of the bodies are then represented by expansions in spherical harmonics (Borderies 1978, Shutz 1981). Moreover, internal structures of solid deformable bodies, homogeneous or with core-mantle interfaces, are represented by several terms and parameters arising from tidal deformations of the bodies (both elastic and anelastic). More details and references on these topics are given in the above quoted papers related to our works concerning the theory of the Moon’s spin-orbit motion.

The BJV model, as described above, has been extended to the spin-orbit integration of the terrestrial planets (Mercury, Venus, Earth, and Mars). This new model is henceforth called SONYR (for Spin-Orbit N-bodY Relativistic model). In the present paper framework, the SONYR model is devoted to the detailed analysis of the complete spin-orbit motion of Mercury.

The simultaneous integration of the solar system, including the Mercury’s spin-orbit motion, uses a global reference system given by the solar system barycenter. Nevertheless, let us recall that local dynamically non-rotating frames show a slow (de Sitter) rotation with respect to the kinematically non-rotating frames. As a consequence, the reference frame for the Mercury’s rotation is affected by a slow precession of its axes transported with the translational motion of Mercury. In the Earth’s case, the de Sitter secular precession of the Earth reference frame is close to 1.92 as/cy (see Fukushima 1991, Bizouard et al. 1992, Bois & Vokrouhlický 1995). Consequently, the real rotation of Mercury has not to be expressed in an inertial system fixed in space, but in a local dynamically non-rotating frame fallen down in the gravitational field of the Sun. Because of the proximity of Mercury to the Sun, its de Sitter precession may be expected quite significant.

In the end, the SONYR model and its analysis method take into account (i) the experience in post-Newtonian gravitation in the definition of reference frames required to deal with rotational motions combined with translational ones, and (ii) the modern knowledge of dynamical systems for studiing librations as quasi-periodic solutions according to the axiomatic presented in Bois (1995). We can state that the model is not Newtonian but rather ‘Newtonian-like’, resulting from truncation of the fully post-Newtonian (DSX) framework. In the present paper, we deal with the Newtonian-like librations (classical physical librations), while the formally relativistic contributions (relativistic librations and de Sitter precession of the Mercury’s reference frame) will be analyzed in a forthcoming paper.

### 3.2 Method

The model is solved by modular numerical integration and controlled in function of the different physical contributions and parameters taken into account. The N-body problem (for the translational motions), the rotational motions, the figure-figure and tidal interactions between the required bodies are simultaneously integrated with the choice of the contributions and truncations at our disposal. For instance, the upper limits of the extended figure expansions and mutual interactions may be chosen as follows : up to $`l=5`$ in the Moon case, $`4`$ for the Earth, $`2`$ for the Sun while only the Earth-Moon quadrupole-octupole interaction is taken into account (see previous papers). The model has been especially built to favor a systematic analysis of all the effects and contributions. In particular, it permits the separation of various families of librations in the rotational motions of the bodies.

The non-linearity features of the differential equations, the degree of correlation of the studied effect with respect to its neighbors (in the Fourier space) and the spin-orbit resonances (in the Moon and Mercury’s cases), make it hardly possible to speak about ‘pure’ effects with their proper behavior (even after fitting of the initial conditions). The effects are not absolutely de-correlated but relatively isolated. However, the used technique (modular and controlled numerical integration, differentiation method, mean least-squares and frequency analysis) gives the right qualitative behavior of an effect and a good quantification for this effect relative to its neighbors. In the case of the particular status of the purely relativistic effetcs, their quantitative behaviors are beyond the scope of the present paper and will be discussed in a forthcoming work. When a rotational effect is simply periodic, a fit of the initial conditions for a set of given parameters only refines without really changing the effect’s behavior. The amplitudes of librations plotted in Figures 11 and 12 are then slightly upper bounds.

The precision of the model is related to the one required by the theory of the Moon. One of the aims in building the BJV model (at present included inside SONYR) was to take into account all phenomena up to the precision level resulting from the LLR data (i.e. at least $`1`$ cm for the Earth-Moon distance, $`1`$ milliarcsecond (mas) for the librations). For reasons of consistency, several phenomena capable of producing effects of at least $`0.1`$ mas had been also modeled (the resulting libration may be at the observational accuracy level). Moreover, in order to justify consistence of the Moon’s theory, this one had been adjusted to the JPL ephemeris on the first $`1.5`$ years up to a level of a few centimeter residuals. In the other hand, the internal precision of the model is only limited by the numerical accuracy of the integration. Thus, in order to avoid numerical divergence at the level of our tests for Mercury, computations have been performed in quadrupole precision ($`32`$ significant figures, integration at a $`10^{14}`$ internal tolerance).

### 3.3 Terminology

In order to de-correlate the different librations of Mercury, we use the terminology proposed in Bois (1995), which is suitable for a general and comparative classification of rotational motions of the celestial solid bodies. This terminology derives from a necessary re-arrangement of the lunar libration families due to both progress in the Moon’s motion observations (LLR) and modern knowledge of dynamical systems.

Traditionally, the libration mode called physical libration is split up according to the conventional dualism ”forced-free”. The forced physical librations are generally related to gravitational causes while the free librations would be departures of the angular position from an equilibrium state. These cuttings out contain ambiguities and redundancies discussed in previous papers (Bois 1995, Bois 2000). Formally, the free librations are periodic solutions of a dynamical system artificially integrable (by a convention of writing related to specific rates of the spin-orbit resonance, for instance 1:1), whereas the forced librations express, in space phase, quasi-periodic solutions around a fixed point (the system is no longer integrable). Moreover, any stable perturbed rotation of celestial solid body contains imbricate librations of different nature, and those are too strongly overlapped to keep the traditional classification.

In the present terminology, the libration nature, its cause and its designation are linked up. Two great libration families serve to define the physical librations, namely the potential librations and the kinetical librations. They simply correspond to a variation energy, potential or kinetical respectively. For libration sub-classes, the designation method is extensive to any identified mechanism (see more details in Bois 1995). The terminology permits easily the separation of various families (see the Moon’s case described in a set of previous papers). These librations are called direct when they are produced by torques acting on the body’s rotation. They are called indirect when they are produced by perturbations acting on the orbital motion of the body. Indirect librations derive from spin-orbit couplings.

A specificity of the SONYR model with its method of analysis is to isolate the signature of a given perturbation. The SONYR model allows indeed the identification of relationships between causes and effects including interactions between physics and dynamics, such as the dynamical signature of a core-mantle interaction (called centrifugal librations).

### 3.4 Parameters and initial conditions

In the computations presented in the paper, the required dynamical parameters and general initial conditions come from the JPL DE405 ephemeris (Standish 1998). However, concerning the parameters related to the Mercury’s rotation (second-degree spherical harmonics $`C_{20}`$ and $`C_{22}`$), which are not included in the JPL ephemeris, our model uses those given by Anderson et al. (1987) (see Table 2). Besides, up to now it does not exist any ephemeris of the Mercury’s rotation. As a consequence, to build initial conditions for the Hermean rotation (described by an Eulerian sequence of angles $`\psi `$, $`\theta `$ and $`\phi `$ defined below in Section 4.1), we use the following principle : Assuming the polar axis of Mercury normal to its orbital plane, we obtain $`\psi =\mathrm{\Omega }`$ and $`\theta =i`$ where $`\mathrm{\Omega }`$ and $`i`$ are respectively the ascending node and the inclination of the orbit of Mercury on the Earth equatorial plane (which is the reference frame used in the DE405 ephemeris). The long axis of Mercury being pointed towards the Sun at its periapse allows to fix the $`\phi `$ angle of polar rotation. The value of $`\dot{\phi }`$ is found in Seidelmann et al. (2002). We use at last $`\dot{\psi }=0`$ and $`\dot{\theta }=0`$; these two variables reach to mean values generated by the complete spin-orbit problem of Mercury : $`0.616\times 10^7`$ deg/day and $`0.267\times 10^8`$ deg/day respectively. The numerical integrations presented in the paper start from these initial conditions related to the planar problem for Mercury; they are listed in Table 1. Departure from the planar case is understood as the integration of physics included in SONYR.

In the other hand, for the computations carried out in this paper, the global reference frame $`O^{}X^{}Y^{}Z^{}`$ is given by a reference system centered on the solar system barycenter, fixed on the ecliptic plane, and oriented at the equinox $`J2000`$. The rotational motion of Mercury is evaluated from a coordinate axis system centered on the Mercury’s center of mass $`Oxyz`$ relative to a local dynamically non-rotating reference frame, $`OXYZ`$, whose axes are initially co-linear to those of $`O^{}X^{}Y^{}Z^{}`$. In the framework of the present paper without purely relativistic contributions, let us note that axes of $`OXYZ`$ remain parallel to those of $`O^{}X^{}Y^{}Z^{}`$.

The N-body problem for the planets of the solar system and the Mercury’s spin-orbit motion are simultaneously integrated in the SONYR model. Concerning the rotational equations written in a relativistic framework, the reader may refer to Bois & Vokrouhlický (1995). In a Newtonian approach, these equations amount to the classical Euler-Liouville equations of the solid rotation (see e.g. Goldstein 1981). We follow the formalism and the axiomatic expanded in Bois & Journet (1993) and Bois (1995) for the definition of the different rotational elements as well as the used terminology. Let us simply precise that l is the angular momentum expressed in $`Oxyz`$ and is related to the instantaneous rotation vector $`\omega `$ as follows :

$$\text{l}=\left(I\right)\omega $$
(9)
where $`\left(I\right)`$ is the tensor of inertia for the body. According to classical assumptions, $`\left(I\right)`$ is reduced to three principal moment of inertia $`A,B`$, and $`C`$. The gravity field of Mercury is essentially unknown. The tracking data from the three fly-by of Mariner 10 in 1974-75 have been re-analyzed by Anderson et al. (1987) to give a low accuracy estimate of the normalized coefficients $`C_{20}`$ and $`C_{22}`$ (see Table 2, values are expressed in the body-fixed frame of the principal axes of inertia). The principal moments of inertia $`A`$ and $`B`$ are then infered from $`C_{20}`$ and $`C_{22}`$ by the following formulae (Ferrari et al. 1980) :

$`{\displaystyle \frac{A}{MR^2}}`$ $`=`$ $`C_{20}2C_{22}+{\displaystyle \frac{C}{MR^2}}`$ (10)
$`{\displaystyle \frac{B}{MR^2}}`$ $`=`$ $`C_{20}+2C_{22}+{\displaystyle \frac{C}{MR^2}}`$
Let us note that the parameter $`C/MR^2`$ is not constrained within these relations. We present in Section 4.2 the variations of $`C/MR^2`$ by analyzing its signature into the librations of Mercury.

## 4 The librations of Mercury

### 4.1 Planetary perturbations

For Mercury assumed to be a rigid body reduced to three oblateness coefficients, the general expression for a torque coming from a disturbing point mass $`m`$ at the vectorial position $`r𝐮`$ ($`r`$ is the instantaneous distance) is written as follow :

$$𝐍_2=\left(\frac{3Gm}{a^3}\right)\left(\frac{a}{r}\right)^3𝐮\times (I)𝐮$$
(11)
where $`a`$ is the mean distance between the two bodies. Such a torque due to the Sun is the dominant one acting on the rotational motion of Mercury. The solar torque is indeed responsible of the global dynamical behavior of the rotation of Mercury. Due to Venus, the resulting torque is about of the order of $`10^5`$ times smaller, that is quite negligible.

Figure 4 presents the rotational motion of Mercury including only the solar torque in the rotational equations and taking into account simultaneously the whole N-body problem for the Sun and the planets (the planetary interactions inducing indirect effects on the rotation of Mercury). In this Figure, the Euler angles $`\psi ,\theta ,\phi `$ related to the 3-1-3 angular sequence describe the evolution of the body-fixed axes $`Oxyz`$ with respect to the axes of the local reference frame $`OXYZ`$. Let us recall the definition used for these angles : $`\psi `$ is the precession angle of the polar axis $`Oz`$ around the reference axis $`OZ`$, $`\theta `$ is the nutation angle representing the inclination of $`Oz`$ with respect to $`OZ`$, and $`\phi `$ is the rotation around $`Oz`$ and conventionally understood as the rotation of the greatest energy (it is generally called the proper rotation). The axis of inertia around which is applied the proper rotation is called the axis of figure and defines the North pole of the rotation (Bois 1992). Let us remark that in the Figure 4 (in other Figures involving $`\phi `$ as well) plotted over 3000 years, we have removed the mean rotation of $`58.646`$ days in the $`\phi `$ angle in order to better distinguish the librations. We may then clearly identify the first proper frequency of $`15.847`$ years (to be compared to the analytical determination, namely 15.830 years given in Section 2).

The $`\psi `$ angle expresses the nodal precession of the equatorial plane of Mercury with respect to the ecliptic plane. It splits up in a periodic term with a period $`\mathrm{\Psi }=1066`$ years and a secular one $`\mathrm{\Pi }`$ = 278 898 years. $`\mathrm{\Psi }`$ is the second proper frequency of the Mercury’s spin-orbit coupling. It can be analytically approximated by the following formula (used in the Earth’s case by Goldstein $`1981`$) :

$$2\pi /\left(\frac{3}{2}\frac{n^2}{\omega _z}\frac{C\overline{A}}{C}\mathrm{cos}\theta \right)$$
(12)
This analytical expression is suitable for an axis symmetric body. It is not the case of Mercury. However, assuming Mercury as a symmetric top rotating about its smallest axis of inertia, with an average equatorial moment of inertia $`\overline{A}=\frac{A+B}{2}`$, one finds for this period 1300 years (by using the values given in Table 2). The difference between the two values permits to appreciate the departure of Mercury from a symmetric body. In the other hand, the dynamical behavior of $`\psi `$ coming from SONYR is due to the direct effect of the solar torque by the way of the true dynamical figure of Mercury. In the 2-body problem, Sun-Mercury, $`\mathrm{\Psi }=1066`$ years and $`\mathrm{\Pi }=0`$ years. Related to the planetary interactions, $`\mathrm{\Pi }`$ = 278 898 years expresses the departure from the 2-body problem.

The dynamical behaviors of $`\psi `$ and $`\mathrm{\Omega }`$ (the ascending node of the orbit) are quite superimposed as shown in Figure 5 (top panel). As it is mentioned by Béletski (1986), a second synchronism is generally involved in a Cassini state. However, in the Mercury’s case, it is not a periodic synchronism as it is the case for the Moon. Using SONYR, we make easily in evidence the periodic 18.6 year synchronism in the lunar spin-orbit resonance while it is about a secular behavior in the Hermean spin-orbit resonance. The 3:2 resonance of Mercury introduces a mechanism of spin-orbit secular resonance qualitatively analogous to the orbital secular resonances combined with mean motion resonances (as for instance in the 2:1 case, see Bois et al. 2003). As it is shown in Figure 5 (bottom panel), the spin-orbit secular resonance variable $`\xi =\psi \mathrm{\Omega }`$ does not present any secular term. As a consequence, $`\psi `$ and $`\mathrm{\Omega }`$ on average precess at the same rate equal to $`\mathrm{\Pi }`$, confirming then the mechanism of spin-orbit secular resonance. We find that $`\xi `$ librates with the particular frequency of 1066 years.

In addition, in order to give a detailed inspection of the short periods involved in the rotational motion of Mercury, Figure 6 presents the solution plotted over 500 days. The $`P_\phi `$ rotation period of 58.646 days appears in the $`\psi `$ and $`\theta `$ angles. Whereas the mean rotation of 58.646 days is removed in the $`\phi `$ angle (as in Figure 4), the signature of the $`P_\lambda `$ orbital period of 87.969 days is clearly visible (this angle is called libration in longitude of 88 days in literature). A third period appears in the $`\psi `$ and $`\theta `$ angles, namely 176.1 days. This one results from the $`3`$:$`2`$ spin-orbit resonance ($`\stackrel{~}{P}=2P_\lambda =3P_\phi `$).

Figure 7 presents the planetary interactions acting on the rotational motion of Mercury by the way of its spin-orbit couplings (i.e. indirect effects of the planets on the Mercury’s rotation). In the black line case the problem is reduced to the Sun and Mercury. In this 2-body problem, the orbital plane does not precess as it is clear in the $`\psi `$ and $`\theta `$ angles without secular terms. The secular variations rise up from the departure of the 2-body problem (as it is visible with the broken, dots, and cross line cases in Figure 7). In the broken line case, the interactions between the Sun, Mercury and Venus are taken into account. With the dots line, the later case includes Jupiter in addition. The whole planetary interactions are integrated in the cross line case (except for Pluto). We show that Venus is the planet which induces the greatest secular term. After Venus, the role of Jupiter is dominant, and this 4-body problem (Sun, Mercury, Venus, and Jupiter) defines our ”standard” case used in our Section 4.3 for the analysis of the Hermean librations. The rate of secular variations in the Mercury’s rotation between all planetary interactions (cross lines) and our standard case (dots lines) is $`11.8`$ as/cy (as : arcseconds) in the $`\theta `$ nutation angle and $`1.9`$ amin/cy (amin : arcminutes) in the $`\psi `$ precession angle. These values should be used as corrective terms in analytical theories of the rotational motion of Mercury. Let us emphasize that the spin-orbit motion of Mercury coming from our standard case is sufficient for preserving the $`3`$:$`2`$ resonance ratio between the two modes of motion.

Starting with the initial conditions defined in Section 3.4 (where in particular the initial obliquity of Mercury is equal to zero), the SONYR model permits obtaining the dynamical behavior of the Hermean obliquity by its simultaneous spin-orbit integration. The variables $`i,\theta ,\mathrm{\Omega },\psi `$ from SONYR substituted in the following relation :

$$\mathrm{cos}\eta =\mathrm{cos}i\mathrm{cos}\theta +\mathrm{sin}i\mathrm{sin}\theta \mathrm{cos}(\mathrm{\Omega }\psi )$$
(13)
produce the instantaneous obliquity $`\eta `$ plotted over 3000 years as presented in Figure 8. Such a behavior for $`\eta `$ gives easily a mean obliquity $`\eta _0=1.6`$ amin. We show at present that this mean value is quite consistent with the Cassini state of Mercury. Let be the following equation established by Peale (1988) and reformulated by Wu et al. (1997) and coming from the Cassini laws :

$$\frac{MR^2}{C}=\frac{\mu }{n}\frac{\mathrm{sin}(i+\eta _0)}{\mathrm{sin}\eta _0[(1+\mathrm{cos}\eta _0)G_{201}C_{22}\mathrm{cos}\eta _0G_{210}C_{20}]}$$
(14)
where $`\mu `$ is the precessional angular velocity of the Hermean orbit, namely $`2\pi /\mathrm{\Pi }`$ $`years^1`$, while

$$G_{201}=\frac{7}{2}e\frac{123}{16}e^3\text{ and }G_{210}=(1e^2)^{\frac{3}{2}}$$

are eccentricity functions defined by Kaula (1966). $`i`$ is the inclination of the orbital plane of Mercury relative to a reference system precessing with the orbit. This inclination varies between from $`5^{}`$ to $`10^{}`$ while the eccentricity varies from 0.11 to 0.24, over $`10^6`$ years (Peale 1988). As a consequence, $`\eta _0`$ obtained with (14) belongs to \[1.33, 2.65\] amin. Let us note that a conventional value of $`7`$ amin is often given in literature. Such a value, outside the interval of possible values, is very probably incorrect as already claimed by Wu et al. (1997).

### 4.2 Principal figure librations

#### 4.2.1 Signature of the $`C/MR^2`$ coefficient on the rotational motion of Mercury

Let us consider at present the disturbing torques acting on the rotational motion of Mercury and as a consequence inducing direct librations. This section focuses on the librations related to the dynamical figure of the planet. Such librations are called principal figure librations (Bois 1995). We assume the Sun reduced to a point mass while the gravity field of Mercury is expanded in spherical harmonics up to the degree 2. We express the solar torque acting on the figure of Mercury according to the equation (11).

The first coefficients of the Hermean gravity field have been determined with the Mariner $`10`$ probe (Anderson et al. 1987). We use these values for $`C_{20}=J_2`$ and $`C_{22}`$ given in Table 2. In order to complete the Hermean tensor of inertia (coefficients $`A,B,C`$), the $`C/MR^2`$ principal moment of inertia is required (see equation (3.4)). Its value is related to the internal density distribution of the planet according to the polar axis of Mercury (rotation of greatest energy about the smallest principal axis of inertia). For an homogeneous planet, such a normalized value is equal to 0.4. We use a nominal value of 0.34 (Table 2) used by Milani et al. (2001) and coming from an internal structure model of Mercury including three layers (crust, mantle and core).

Figure 9 presents the rotational behavior of Mercury computed over $`\mathrm{10\hspace{0.17em}000}`$ years in our standard case with two different values of its greatest principal moment of inertia : (i) the computation with $`C/MR^2`$=0.4 is plotted with the dashed lines (homogeneous planet); (ii) the broad lines are obtained with $`C/MR^2`$=0.34 (three layer model of Mercury). On the bottom panel ($`\phi `$ angle), the dashed lines are shifted away 1 degree in order to distinguish them from the broad lines. This Figure 9 (especially the $`\theta `$ angle) shows how the $`C/MR^2`$ coefficient value and the constant of precession $`\mathrm{\Psi }`$ (2<sup>nd</sup> proper frequency) are related. With $`C/MR^2`$=0.4, $`\mathrm{\Psi }=1254.01`$ years while with $`C/MR^2`$=0.34, $`\mathrm{\Psi }=1066`$ years. Besides, $`C/MR^2`$ and $`\mathrm{\Phi }`$ (1<sup>st</sup> proper frequency) are also linked and for evaluating this relation, Figure 10 shows the variations $`\mathrm{\Delta }\phi `$ obtained by differentiation : (i) on the top panel $`\mathrm{\Delta }C=1\%`$, (ii) on the middle panel $`\mathrm{\Delta }C=2\%`$, (iii) on the bottom panel $`\mathrm{\Delta }C=15\%`$ ($`0.34+15\%(0.34)=0.4`$). The beats are signatures related to the variations in $`\mathrm{\Phi }`$.

#### 4.2.2 Signature of the $`C/MR^2`$ coefficient on the obliquity and on the libration angle in longitude

One of the main objectives of the BepiColombo and MESSENGER missions is to measure the rotation state of Mercury, up to an accuracy allowing to constrain the size and physical state of the planet’s core (Milani et al. 2001, Solomon et al. 2001). Consequently, the two missions have to determine the four following parameters : $`C_{20}`$, $`C_{22}`$, $`\eta `$, and $`\phi `$ that are sufficient to determine the size and state of the Mercury’s core (see Peale 1988, 1997). Combining $`C_{20}`$, $`C_{22}`$, and $`\eta `$, one obtains the $`C/MR^2`$ coefficient while with $`C_{22}`$ and $`\phi `$, one obtains $`C_m/MR^2`$ (i.e. the $`C/MR^2`$ coefficient for the mantle). The validity condition of the first combination is that the dynamical behavior of the core has to follow the one of the mantle over a period of time at least the one of $`\mathrm{\Pi }`$ (assertion 1). The validity condition of the second combination is that the dynamical behavior of the core has not to be coupled to the one of the mantle over a period of 88 days ($`P_\lambda `$) (assertion 2). These two conditions linked together imply some constraints on the nature of the core-mantle interface (Peale 1997). In order to reach such an objective, the BepiColombo mission has to obtain a value on the $`C/MR^2`$ coefficient with an accuracy of 0.003, i.e. $`1\%`$ and therefore foresees measuring the libration angle and the obliquity with an accuracy of 3.2 and 3.7 as respectively (Milani et al. 2001).

Our SONYR model gives (i) the true relation between the three parameters ($`C/MR^2`$, $`\eta `$, $`\phi `$), and (ii) the upper bounds of the impact of $`C_m/MR^2`$ on the $`\phi `$ angle. Figure 11 presents the impact of the $`C/MR^2`$ coefficient on the instantaneous obliquity $`\eta `$. In these plots (Figures 11 and 12), the spin-orbit motion of Mercury is integrated within the whole solar system with an initial obliquity of 1.6 amin, which is the mean obliquity of Mercury evaluated in Section 4.1. The top panel of Figure 11 expresses the dynamical evolution of $`\eta `$ computed over 500 days with $`C/MR^2`$=0.34 (black lines) and $`C/MR^2`$=0.3434 (dashed lines). Dashed lines are shifted from 0.001 amin in order to distinguish the two different kinds of lines. Figure 11 shows also how the instantaneous obliquity of Mercury differs from its $`1.6`$ amin nominal value. The bottom panel shows by differentiation the signature of the $`1\%`$ variation of $`C/MR^2`$ on $`\eta `$. The maximal amplitude crest to crest is of the order of 2.3 mas within the characteristic period of 176.1 days related to the $`3`$:$`2`$ ratio of the Mercury’s spin-orbit resonance.

Figure 12 presents the signature of the $`C/MR^2`$ coefficient on the $`\phi `$ libration angle in longitude. The top panel expresses the behavior of $`\phi `$ computed over 500 days with $`C/MR^2`$=0.34 (black lines) and $`C/MR^2`$=0.3434 (dashed lines). One may compare this Figure 12 to the Figure 3 resulting from the usual analytical resolution of the Eulerian equation 4 (thanks to G. Giampieri, private communication). Let us note that the angle $`\gamma `$ defined in Figure 3 is equal to the angle $`\phi `$ plotted in Figure 12. The later only gives a simple double sine curve with an amplitude of 42 as while the solution of the SONYR model includes the couplings between the three rotational variables as well as the indirect couplings due to planetary interactions (we notice that in the two Figures 12 and 3 the amplitude of libration is of the order of 40 as). Let us note that the Figure 12 corresponds to the libration related to the $`C/MR^2`$ coefficient of the planet without core-mantle couplings. Let us add that in Peale (1972), the amplitude of $`\phi `$ is related to the $`C_m`$ coefficient by assuming that the assertion 2 quoted upper is true. On the contrary, our first results on this topic make in evidence the existence of a faint coupling. This core-mantle coupling will be presented in a forthcoming paper. The bottom panel of Figure 12 shows by differentiation the signature of the $`1\%`$ variation of $`C/MR^2`$ on $`\phi `$. The maximal amplitude within the period of about 88 days (i.e. the signature of $`P_\lambda `$) is of the order of $`0.45`$ as.

In conclusion, signatures of the indeterminacy of $`1\%`$ in $`C/MR^2`$ on the obliquity and on the libration in longitude are 2.4 mas and 0.45 as respectively. What is very faint (may be too much) with respect to the expected accuracy forecasted in the BepiColombo mission.

### 4.3 variing the obliquity

Because the initial obliquity value is unknown, we test in this last section the impact of the indeterminacy of this value on the spin-orbit motion of Mercury. The results are presented in Figure 13 plotted over 3000 years; top panel: the effect on the nutation angle $`\theta `$, middle panel: the effect on the orbital inclination $`i`$, and bottom panel: the effect on the instantaneous obliquity $`\eta `$. On each panel, three curves are related to three different initial values of $`\eta `$, namely 0 amin (black lines), 1 amin (dashed lines), and 2 amin (dot lines). In the bottom panel, the amplitudes of these librations are of the order of 1.4 amin with a period of 1066 years. For any initial value of $`\eta [0,3.2]`$ amin, the mean value of $`\eta `$, let be $`\eta _0`$, is equal to 1.6 amin, which is in good agreement with the determination of $`\eta _0`$ in a previous section. We may claim that $`\eta _0`$=1.6 amin.

For obtaining such a mean obliquity by measurements, let us underline that the theoretical behavior of $`\eta `$ points out to fit the observations by a sine function taking into account the long period $`\mathrm{\Psi }=1066`$ years with an amplitude of 1.6 amin.

## 5 conclusion

The 3:2 spin-orbit resonance between the rotational and orbital motions of Mercury results from a functional dependance of the tidal friction adding to a non-zero eccentricity and a permanent asymmetry in the equatorial plane of the planet. The upcoming space missions, MESSENGER and BepiColombo with onboard instrumentation capable of measuring the rotational parameters stimulate the objective to reach an accurate theory of the rotational motion of Mercury.

Starting from our BJV relativistic model of solar system integration including the coupled spin-orbit motion of the Moon, we have obtained a model generalizing the spin-orbit couplings to the terrestrial planets (Mercury, Venus, Earth, and Mars). The updated model is called SONYR (acronym of Spin-Orbit N-BodY Relativistic model). It permits to analyze and identify the different families of rotational librations. This work has been carried out for Mercury in the present paper.

The spin-orbit motion of Mercury is characterized by two proper frequencies (namely $`\mathrm{\Phi }`$ = 15.847 and $`\mathrm{\Psi }`$ = 1066 years) and its 3:2 resonance presents a second synchronism which can be understood as a spin-orbit secular resonance, ($`\mathrm{\Pi }`$ = 278 898 years). A new determination of the mean obliquity has been proposed in the paper. By using the SONYR model, we have found a mean obliquity of 1.6 amin. This value is consistent with the Cassini state of Mercury. Besides, we have identified in the Hermean librations the impact of the uncertainty of the greatest principal moment of inertia ($`C/MR^2`$) on the obliquity and on the libration in longitude (2.3 mas and 0.45 as respectively for an increase of 1$`\%`$ on the $`C/MR^2`$ value). These determinations prove to be suitable for providing constraints on the internal structure of Mercury. The direct core-mantle interactions will be presented in a forthcoming paper.

###### Acknowledgements.

The authors thank Aleksei Pavlov for his help in the Poincaré cross-section computations and Jean Brillet for providing his efficient method of mean last squares useful for accurate determinations of periods in our data files.

## Appendix A free rotation

Let us assumed Mercury isolated in space; in this sense, its rotation is free and the Euler-Liouville equations for its rotation are written without right hand side, i.e. without any external disturbing torques. If we add the assumption of a rigid Mercury, we are in the Euler-Poinsot motion case (whose solutions are the well-known Eulerian oscillations). Without explicitly integrating such equations, the assurance of integrability in the Poincaré sense can be obtained by some theoretical simple considerations. Indeed, whatever being the triplet of generalized coordinates used to describe the spatial attitude of a solid body in a fixed frame, one knows that there exists four independant integrals of motion : the hamiltonian $`H`$, and the three components $`L_X,L_Y,L_Z`$ of the angular momentum (in $`OXYZ`$). Four integrals of motion for three degrees of freedom, the problem is then integrable and even over-integrable. One does not lose the generality of the problem choosing for instance $`ABC`$. The choice $`AB<C`$ makes possible to write the general solution of the system under a form involving the elliptical functions of Jacobi (Landau & Lifchitz 1969). By convention, let us adopt that the resulting oscillations in space be called the Eulerian oscillations, expressing exclusively the oscillations of the non-perturbed rotation of the rigid body. From this resolution, we obtain the Eulerian frequencies :

$$\mathrm{\Omega }_f=\sqrt{\frac{L^22AH}{C(CA)}}$$
(15)
and

$$\sigma =\mathrm{\Omega }_f\sqrt{\alpha \beta }$$
(16)
that give for Mercury the two periods $`\mathrm{\Omega }_f`$ and $`\sigma `$ (where $`L`$ is the angular momentum in $`OXYZ`$, $`H`$ the energy). $`\alpha `$ and $`\beta `$ are the dynamical coefficients of the body related to $`A,B`$ and $`C`$ by the following relations :

$$\alpha =\frac{CB}{A}\text{ and }\beta =\frac{CA}{B}$$
(17)
With the previous initial conditions given in Table 1 ($`\dot{\psi }=\dot{\theta }=0`$), the free rotation of Mercury is reduced to an elementary rotation of the $`\phi `$ angle only ($`\mathrm{\Omega }_f=\omega _z`$). Starting with the physical mean values of $`\dot{\psi }`$, $`\dot{\theta }`$, and $`\eta _0`$ evaluated in the present paper, $`\mathrm{\Omega }_f=58.646`$ days and $`\sigma =964.88`$ years. The period of $`58.646`$ days corresponds to the Hermean polar rotation while the one of $`964.88`$ years means the global period of rotation of Mercury in $`R^3`$, which is equivalent to the Earth’s classical Euler period, namely 305 days (Lambeck 1980). In the case of the Moon, this period is equal to $`148.129`$ years.

Using the SONYR model reduced to the free rotation of Mercury, we obtain the components of the instantaneous rotation vector $`\omega `$ in the body-fixed system $`Oxyz`$, as presented in Figure 14. $`\omega _z`$ is well found constant while the Euler period is equal to 964.92 years.


format: text

text:
# Chapter 1 Probing TeV gravity with extensive air-showers

## Chapter 1 Probing TeV gravity with extensive air-showers

### Abstract

Particle collisions with center-of-mass energy larger than the fundamental gravitational scale can generate non perturbative gravitational objects such as black holes and branes. In models with large extra dimensions, the fundamental gravitational scale may be around a TeV, making it possible for next generation particle colliders and ultra-high energy cosmic rays to produce such non perturbative gravitational objects. The decay of TeV gravitational objects is significantly different from standard model processes such that probes of these new ideas are within reach. We study the differences between standard model and TeV gravity interactions in extensive air showers (EAS) generated by ultra-high energy cosmic neutrinos. We show that discriminating TeV gravity from standard model interactions is generally difficult, but not impossible given a few unique signatures.

### 1. Introduction

In models with large extra dimensions (LEDs) \[6-8\], the fundamental scale of gravity may be $``$ TeV. In these models, particle collisions with center-of-mass (CM) energy larger than $``$ TeV may create non perturbative gravitational objects such as black holes (BH)  and branes . Next generation particle colliders  and interactions of ultra-high energy cosmic rays (UHECRs) with the atmosphere  can reach TeV CM energies and, therefore, create these non perturbative gravitational objects (see  for a complete review).

The shower-to-shower fluctuations in EASs initiated by hadronic primaries combined with the small branching ratio to BH formation makes the study of TeV gravity with ultra-high energy protons hopeless . Ultra-high energy neutrinos produced by the photo-pion production of ultra-high energy protons in the cosmic microwave background provide a cleaner beam to test departures from SM interactions . Here, we discuss the characteristics of EASs initiated by cosmogenic neutrinos that produce BHs in TeV gravity theories and contrast these with standard model (SM) interactions (based on a more detailed study reported in ). We find that due to the long interaction length of high energy neutrinos in the atmosphere and the large uncertainties in the BH formation cross section discriminating between the two scenarios is quite difficult. However, we also find that a few unique signatures involving tau leptons can help detect TeV BHs.

### 2. Black hole production in TeV gravity

The 4-dimensional Planck mass in natural units can be written as $`M_{Pl}G_4^{1/2}`$, where $`G_4`$ is the 4-dimensional gravitational constant. In the presence of extra dimensions, the fundamental Planck mass is given by $`M_{}=G_{n+4}^{1/(n+2)}`$, where $`n`$ is the number of extra dimensions. The 4-dimensional and $`(n+4)`$-dimensional gravitational constants are related by $`G_4=G_{n+4}/V_n`$, thus $`M_{Pl}^2=M_{}^{n+2}V_n`$ , where $`V_n`$ is the volume of the extra dimensions. In models where $`V_n`$ is large, $`M_{}`$ can be $``$ TeV. Depending on details of LED models (such as $`n`$, $`M_{}`$, and minimum BH formation mass, $`M_{BH,min}`$), the cross section for BH formation in neutrino-nucleon interactions can be either enhanced or suppressed with respect to SM cross section by several orders of magnitudes. Fig. 1 shows an example of the range in BH formation and SM cross sections, which include uncertainties in $`M_{BH,min}`$ and the parton distribution functions (PDFs)  (see  for other examples). If BHs form in TeV CM collisions, they evolve by shedding their hair first followed by Hawking evaporation, leading to the emission of SM particles such as quarks and gluons that subsequently hadronize into jets.

### 3. Extensive air-shower simulations

We simulated both SM showers and BH showers to compare their detectable characteristics. The most relevant SM process for the comparison is the $`\nu _e`$ charged current (CC) interaction which produces an electron shower that we modeled with AIRES . For the BH production process, the secondary particles from BH evaporation are hadronized with PYTHIA  before producing a shower with AIRES.

In Fig. 2, we show showers generated by neutrino primaries with energy $`E_\nu =10^7`$ TeV and TeV gravity parameters $`n`$=6, $`M_{BH,min}=2M_{}=2`$ TeV. On the left panel, the first interaction point ($`X_0`$) is fixed to be 10 km above sea level with shower zenith angle of $`70^{}`$, corresponding to a slant depth of 780 $`gcm^2`$, for both BH and SM showers ($`X_0^{CC}=X_0^{BH}`$). $`\nu _e`$CC and BH showers are clearly different by $`200`$ g/cm2. However, unlike the case of UHE protons, the interaction length for neutrinos in the atmosphere is large, thus, $`X_0^{CC}`$ is not fixed. By shifting $`X_0`$ such that the shower maxima, $`X_m`$, for both cases match ($`X_m^{CC}=X_m^{BH}`$), the differences in shower development are much harder to distinguish as seen in the right panel of Fig. 2. Given a large number of neutrino horizontal showers, one can distinguish the SM and the BH cases by studying the rise of the shower, $`X_mX_{0.1}`$, where $`X_{0.1}`$ is the slant depth containing 10 % of particles of $`X_m`$. In addition, the muon content of BH showers is larger than the SM case, since BHs produce hadrons while CC showers do not. A deep horizontal shower accompanied by many muon secondaries is a sign of BH formation.

Detecting TeV BH formation with UHECR detectors may be possible through the decay of $`\tau `$-leptons generated by $`\nu _\tau `$’s that interact in the Earth or in mountain ranges close to the detectors. A secondary $`\tau `$ generated through the decay of a BH has much less energy than the SM $`\tau `$ secondary. In addition, BHs may produce multiple $`\tau `$-leptons in their evaporation, a unique signature of TeV gravity. SM processes that generate multiple $`\tau `$-leptons are highly unlikely, the detection of multiple $`\tau `$’s in earth-skimming and mountain crossing neutrinos will be a smoking gun for BH formation.

### 4. Conclusions

We showed that given the uncertainties in the $`\nu `$-nucleon cross section for TeV gravity and the flux of cosmogenic neutrinos, distinguishing TeV LED models from the SM via the rate of neutrino induced EASs is unattainable. Although BH showers develop faster that the SM ones leading to a difference of 200 g/cm2 in $`X_mX_0`$, the variation in $`X_0`$ for neutrino showers make the distinction quite subtle. A large number of neutrino EASs with measured muon content and $`X_mX_{0.1}`$ is necessary for a clear distinction of SM and BH formation. Finally, a few background free signatures such as multiple $`\tau `$’s and lower energy $`\tau `$ secondaries may more clearly signal the existence of TeV LED models. (We thank NSF and DOE for financial support.)

### 5. References

1. Ahn, E. J., Cavaglià, M., and Olinto, A. V. 2003, Phys.Lett.B., 551, 1.

2. Ahn, E. J. and Cavaglià, M. 2002, Gen.Rel.Grav., 34, 2037.

3. Ahn, E. J., Ave, M., Cavaglià, M., and Olinto, A. V. 2003, submitted.

4. Anchordoqui, et al., 2002, Phys.Rev.D., 65, 124027.

5. Anchordoqui, L. and Goldberg, H., 2002, Phys.Rev.D, 65, 047502.

6. Antoniadis, I. 1990, Phys.Lett.B., 246, 377.

7. Antoniadis, I., et al., 1998, Phys.Lett.B., 436, 257.

8. Arkani-Hamed, N., Dimopoulos, S., and Dvali, G. 1998, Phys.Lett.B., 429, 263.

9. Banks, T. and Fischler, W. 1999, hep-th/9906038.

10. Cavaglià, M. 2003, Int.Jou.Mod.Phys.A., in press, hep-ph/0210296.

11. Dimopoulos, S. and Landsberg, G. 2001, Phys.Rev.Lett, 87, 161602.

12. Feng, J. L. and Shapere, A. D. 2002, Phys.Rev.Lett, 88, 021303.

13. Giddings, S. and Thomas, T. 2002, Phys.Rev.D., 65, 056010.

14. J. Pumplin,et al. 2002, JHEP., 0207, 012

15. Sciutto, S. J. , 1999, astro-ph/9905185

16. Sjöstrand, T., et al., 2001, Computer Physics Commun., 135, 238.

17. Tyler, C., Olinto, A. V., and Sigl, G., 2001, Phys. Rev. D 63, 055001.


text:
# 1 Introduction

## 1 Introduction

Recently, intensive development of techniques and new methods in astrophysical observations had lead to discovery of more than 80 giant black holes in centers of galaxies , . The masses of the black holes are distributed in the range $`210^6310^9M_{}`$ ($`M_{}=210^{33}`$g is solar mass), as is shown at Fig.1. The analysis of observational data revealed a connection between the giant black hole mass in the center of a galaxy and mass of the galaxy’s bulge.

We remind that bulge is the most dense central spherical component of a galaxy, containing mainly old stars. The total rotational velocity of stars in bulge is usually significantly smaller than velocity dispersion $`\sigma `$, that is why the form of bulge is often close to spherically symmetric. The density of stars in bulge grows rapidly to the center. The size of bulge $`r_b130`$kpc, it grows with the galaxy mass. The masses of giant black holes $`M_{bh}`$ are approximately three orders of magnitude smaller than the masses of host galaxy bulges.

According to modern knowledge, the main part of matter in the Universe is made up by nondissipative nonbaryon dark matter. The initial density fluctuations of the dark matter in the early Universe had grown up forming large-scale nonlinear gravitationally bound objects - galactic haloes. The important feature of these haloes is singular distribution of density in their centers . Baryon matter trapped by gravitational field of a halo has settled gradually to its center and formed a galaxy. The presence of singularity in the center of dark matter distribution allows the compressing baryon gas to form a primordial black hole  with mass of the order of $`10^3M_{}`$. Note, that in this scenario primordial black hole is situated exactly in the galaxy dynamical center. Later the seed black hole grows intensively due to accretion both baryonic and dark matter from the bulge. The present paper is devoted to the analysis of evolution of dark matter and its absorption by the black hole in the center of a galaxy.

The paper is organized as follows.

In the second section, the distribution of dark matter in a galaxy and its evolution under the influence of the averaged self-consistent gravitational field of baryon matter (bulge) is considered. The problem is solved in the adiabatic approximation in assumption of spherical symmetry of baryon matter distribution. The role of deviations from spherical symmetry is analyzed. It is shown that even if the baryon matter distribution is noticeably asymmetric, the distribution of dark matter defined by joint action of gravitational fields of baryon and dark matter remains spherically symmetric.

In the third section, growth of a seed black hole in galactic center on account of direct capture of dark matter particles moving in the averaged self-consistent potential is investigated. It is shown that the mechanism is of small efficiency in real conditions and cannot lead to significant growth of the seed black hole mass.

In the forth and fifth sections, the evolution of dark matter distribution function due to the gravitational scattering of particles by stars is considered. The scattering lead to violation of adiabatic approximation. The important peculiarity of the process is that the mean free path length of particles is much more than the size of the bulge, hence the collisions are rare. Taking this into account, in the section 4 we derive the collision integral averaged over dark matter particles oscillations, and write down the kinetic equation describing the evolution of dark matter distribution function. It is shown that because of special features of initial distribution function, the main role is played by diffusion in angular momentum space.

In the fifth section, the solution of the diffusion equation is obtained, the flux of dark matter on black hole is derived, and the law of the black hole’s growth is ascertained.

In conclusion, we make brief analysis of observational data showing that the growth of giant black holes on account of dark matter absorption should be essential.

On the whole, the presented theory is in reasonable agreement with existing observational data concerning giant black holes. Further development both of the observations and the theory and their detailed comparison are undoubtedly of great interest.

## 2 Adiabatic dynamics of gravitating dark matter and baryonic matter.

Observation show that over distances of the order of the horizon radius our Universe is homogeneous, isotropic and expands uniformly. The expansion leads to a rapid cooling of matter. A cold gravitating gas is unstable because of the action of universal gravitational forces. Growth of the Jeans instability creates regions of strong compression with dimensions much smaller than the horizon radius $`R_H`$. This is of decisive importance for the formation (in the CDM scenario) of the large-scale structure of matters in the Universe.

The main role in this process is played by dark matter, which manifests itself only in the gravitational interaction. According to the nowadays knowledge the baryonic matter in the Universe contribute $`\mathrm{\Omega }_B<0.05`$ and thus the gravitational processes are determined mostly by nonbaryonic dark matter. It is supposed as usually that the particles of which dark matter is composed interact very weakly with one another and with baryonic matter.

After recombination initial linear density fluctuations $`\delta \rho `$ are growing relatively slowly with time $`t`$ (due to expansion of Universe) according to power law: $`\delta \rho /\rho t^{2/3}`$. As soon as the disturbances amplitude reaches unity $`\delta \rho /\rho 1`$ the nonlinear processes dominate. Well known Zeldovich-Harrison spectrum of the initial linear fluctuations is growing with the wave vector $`k`$: $`\delta \rho (k)k`$. Due to this nonlinear stage $`\delta \rho (k_n)/\rho 1`$ is reached first for the small $`k`$. Nowadays the nonlinear boundary $`l_nk_{n}^{}{}_{}{}^{1}10100Mps`$. Thus the following conditions are fulfilled for the dark matter large-scale structures $`l`$, which have passed the nonlinear stage of the Jeans instability:

$$\frac{l}{l_f}0,\frac{l}{R_H}<<1,\frac{u}{c}\frac{l}{R_H}<<1.$$
(1)
Here $`l_f`$ is the free path of dark matter particles which we consider as noninteracting, $`u`$ \- characteristic velocities in the inhomogeniety. We see that the nonlinear stage of Jeans instability is developing in noninteracting, nonrelativistic, cold gravitating gas.

The general system of equation, describing the nonlinear gravitational compression of the dark matter can be presented in a form :

$$\frac{\rho }{t}+\frac{}{𝐫}\left(\rho 𝐮\right)=0$$

$$\left(\frac{}{t}+𝐮\frac{}{𝐫}\right)𝐮+\frac{\psi _d}{𝐫}=0$$
(2)
$$\mathrm{\Delta }\psi _d=4\pi G\rho $$

Here $`\rho `$ is the density, $`u`$ \- velocity of the gas, and $`\psi _d`$ \- gravitational potential. The nonlinear process is most strongly developing near the density maximums. Let us consider the vicinity of the one maximum, presenting the initial density distribution $`\rho `$ in the general form:

$$\rho _{t=0}(𝐫)=\rho _0\left(1\frac{x^2}{a^2}\frac{y^2}{b^2}\frac{z^2}{c^2}\right),𝐮_{t=0}=0$$
(3)
The solution of hydrodynamical system of equations (2), lead to singularity at the moment of time $`t=t_k`$. After that moment in a dark matter gas the multifluid motion is developed. It is accompanied by density oscillation. With time $`t\mathrm{}`$ the number of fluids is infinitely growing. The density oscillation amplitudes $`a_i(t)`$ at the same time is going to zero: $`a_i(t)O(1/t)`$. As a result of this kinetic mixing of fluids the stationary spherically symmetric cluster of dark matter is established. The main peculiarity of dark matter distribution in this cluster is the existence of a singularity in its core

$$\rho =Kr^{12/7}$$
(4)
This result was obtained in analytic theory by Gurevich and Zybin .

It should be noted that numerical calculation could not receive this result for a long time and showed only the plane structures of the Zeldovich pancake type. Only the usage of special numerical methods allowed to make a significant step demonstrating the existence of the stationary spherical symmetric structures . But the singularity of density distribution was different from (4). Further development of numerical methods allowed to improve this result and move on very close to the law (4) . So nowadays one can say that the numerical calculation fully confirmed the existence of a stable spherical clusters with the singularity law (4). This cluster was called nondissipative gravitational singularity (NGS).

The main method of analytical solution  was based on the solution of kinetic equation for dark matter particles using canonical action variables $`I=(I_R,m)`$, where

$$I_R(E,m)=\frac{\sqrt{2}}{\pi }_r_{}^{r_+}\sqrt{E\psi _d(r)\frac{m^2}{2r^2}}𝑑r,$$
(5)
\- radial action, $`m`$ is the angular momentum, $`E`$ \- energy, $`r_\pm `$-turning points (roots of the equation $`E\psi _d(r)m^2/2r^2=0`$). According to  the established distribution function is

$$f(I_R,m)=f_0I_R^{1/8}\delta \left(m^2l_{0}^{}{}_{}{}^{2}I_R^2\right),$$
(6)
where $`f_0`$ and $`l_0`$ are constants, determined by the initial distribution (3). In particular the case of spherical symmetric initial density distribution (3): $`a=b=c`$, is degenerate. The spherical symmetry is conserved at all stage of evolution, including initial compression, hydrodynamical singularity, kinetic mixing. In this case the angular momentum are always equal to zero, consequently in (6) $`l_0=0`$. In the case of weak initial asymmetry:
$`(ab)^2+(ac)^2<<a^2`$ parameter $`l_0`$ is small, proportional to $`ϵ`$:

$$l_0=0.16ϵ,ϵ=\frac{1}{a}\sqrt{(ab)^2+(ac)^2(ab)(ac)}$$
(7)
Accordingly as follows from (6) angular momentum in this case is also small and proportional to parameter $`ϵ`$ and to radial action $`I_R`$.

We will consider now the effect of baryonic matter on the distribution of dark matter. Let us take into account, that the baryonic matter particles lose their energy due to inelastic collisions accompanied by optic emission. That is why the baryonic matter falls down to the bottom of the potential wells, formed by the cold dark matter. The luminous baryonic matter then acts as indicator of the dark matter structure identifying in particular the position of the centre of NGS. This produces a unique object: baryonic matter with a halo of dark matter. Examples of such objects are galaxies: in this case the evidence for the presence of a dark matter halo is provided by flat rotation curves . Other examples of such objects are clusters of galaxies, which contain a trapped hot gas, again directly confirming the existence of a dark matter halo .

It should be pointed out that as the baryonic matter drops to the bottom of potential wells, its density rises at the central regions of the wells and can become equal to the density of nonbaryonic matter. That will lead to the considerable distortion of the canonical distribution of the dark matter (4). The main feature of the problem is that the baryonic matter loses its energy and falls to the wells bottom slowly - at the timescales comparable with the lifetime of the Universe. The dark matter particles trapped in the well oscillate many times. It means that the radial action $`I_R`$ is integral of motion (adiabatic invariant) and distribution function $`f_d`$ is conserved and described as previously by formula (6). But the dark matter potential in (5) should be changed to the full potential $`\mathrm{\Psi }`$, determined both by baryonic and dark matter $`\mathrm{\Psi }=\psi _d+\psi _b`$. Consequently the density $`\rho _d`$ takes the form:

$$\rho _d(r,t)=\frac{\sqrt{2}\pi }{r^2}_0^{\mathrm{}}𝑑m^2_{\mathrm{\Omega }_E}𝑑E\frac{f(I_R,m)}{\sqrt{E\mathrm{\Psi }(r,t)m^2/2r^2}},$$
(8)
$`\mathrm{\Omega }_E`$ is the region, where $`E\mathrm{\Psi }(r,t)m^2/2r^2>0`$.

$$I_R=\frac{\sqrt{2}}{\pi }_r_{}^{r_+}\sqrt{E\mathrm{\Psi }(r,t)m^2/2r^2}𝑑r.$$
(9)
Let us consider central part of the object where $`M_b>M_d`$ and $`\mathrm{\Psi }\psi _b`$ . Potential $`\psi _b(r)`$ is determined by the mass $`M_b(r)`$ of a baryonic matter. $`M_b(r)`$ is slowly changing with time $`t`$ due to slowly falling down into the well of a barionic component. Let us suppose, for example, that $`M_b`$ has a power-like dependence on $`r`$:

$$M_b(r)r^n,\psi _br^{n1},n0.$$
(10)
From (6), (8), (9) and (10) it follows:

$$\rho _d=\rho _0r^{(9n39)/16}$$
(11)
Thus we see that in conditions (10) the power scaling for $`\rho _d`$ exists also and just as in the case of purely dark matter (4), the scaling law (11) does not depend on parameters of dark matter initial distribution $`l_0`$ and $`f_0`$. For example, if $`M_b(r)r^{9/7}`$ both $`\rho _d\rho _br^{12/7}`$. If $`\rho _br^2`$, $`M_br`$ (isothermal sphere) and $`\rho _dr^{15/8}`$. We see, that dark matter is compressed by the gravitational action of baryonic gas. For point-like baryonic mass distribution could be reached maximal compression of a dark matter

$$\rho _dr^{39/16}.$$
(12)
We note that the mass distribution $`\rho r^m`$ with $`m>2`$ must lead to the black hole creation and the problem should be considered in more details taking this fact into account.

We discussed above the spherically symmetric distribution of baryonic matter. In real conditions due to rotational motion baryonic component even in the central parts both of the spiral and elliptic galaxies can have a noticeable nonspherical components. Can they lead to a significant discrepancy of dark matter distribution in comparison with the case of purely spherical symmetric baryonic matter? There are some general theorems concerning the behaviour of systems close to integrable ones. According to these theorems, the average action variables $`I_R,m,m_z`$ remain constant during particles motion in weakly asymmetric, slowly varying potential . That is why the distribution function (6) does not change significantly during evolution of the slowly rotating baryonic matter .

## 3 Black hole growth in the center of the dark matter halo (self-consistent approximation).

The aim of the present section is the investigation of possibility of seed black holes growth within the bounds of self-consistent, collisionless approximation.

The growth of central black hole leads to decreasing of dark matter particles number in the loss cone. On the other hand, it leads to the growth of the loss cone. As the result of the two competing processes the flux of a dark matter to black hole can decrease. Let us assume that the growth of black hole stops at total mass

$$M_{bh}=M_b+M_d,$$
(13)
where $`M_b`$ is an initial baryonic mass of a black hole and $`M_d`$ is absorbed dark mass. The flux of dark matter on a black hole with mass $`M_{bh}`$ is composed of particles with angular momentum

$$m<m_g=2cr_g,$$
(14)
where $`r_g=\frac{2M_{bh}}{c^2}`$ is the gravitational radius of a black hole.

In spherically symmetric potential $`\mathrm{\Psi }(r,t)`$ (with arbitrary time dependence) the angular momentum of particles does not change. That is why total mass of particles captured by a black hole is described by initial distribution function $`f_i`$:

$$N(f_i;M_{bh})=d^3rd^3vf_i(𝐫,𝐯)\theta \left(𝐦_𝐠𝐦\right),$$

where

$`\theta (x)=1`$, if $`x>0`$, and $`\theta (x)=0`$, if $`x<0`$. The input of a dark matter particles in total mass of black hole is

$$M_d=N(f_i;M_b+M_d)$$
(15)
Thus the scenario of growth of black hole mass caused by dark matter particles flux depends on initial distribution function $`f_i`$ and black holes baryonic mass $`M_b`$ only. If the equation (15) has the solution $`M_d=M_d(f_i;M_b)>0`$ the mass of a black hole will not exceed $`M_b+M_d(f_i;M_b)`$. If the equation (15) has no solution, the seed black hole will grow beyond all bounds.

As an example let us consider initial isothermal distribution of a dark matter. In this case the distribution function, density and potential have the form:

$$f_i(E,m)=\frac{\rho _0}{\left(2\pi \sigma _d^2\right)^{3/2}}e^{\frac{E}{\sigma _d^2}},\rho _0=\frac{\sigma _d^2}{2\pi G},$$
(16)
$$\rho (r)=\frac{\rho _0}{r^2},$$
(17)
$$\psi _d(r)=4\pi G\rho _0\mathrm{ln}(r),$$
(18)
where $`\sigma _d`$ is velocity dispersion of dark matter particles. Density of particles with angular momentum smaller than $`m_g`$ is

$$\rho _g=\frac{\sqrt{2}\pi }{r^2}𝑑E𝑑m^2\frac{f(E,m^2)}{\sqrt{E\psi _d\frac{m^2}{2r^2}}}\theta \left(E\psi _d\frac{m^2}{2r^2}\right)\theta (m_gm)$$

or, after calculations

$$\rho _g(r)=\frac{\rho _0}{r^2}\left(1e^{\frac{m_g^2}{2\sigma _d^2r^2}}\right),$$

Hense total mass of particles with the angular momentum smaller than $`m_g`$ is

$$N_g=2\pi \left(\frac{\rho _0}{G}\right)^{1/2}m_g.$$
(19)
Taking into an account (19) one can rewrite the equation (15) in the form:

$$\frac{M_d}{M_b}=Q\left(1+\frac{M_d}{M_b}\right),$$
(20)
where

$$Q=\frac{8\pi }{c}\left(G\rho _0\right)^{1/2}$$

The positive solution of the equation (20) exists only at $`Q<1`$. Thus, scenario of growth of a black hole in the case of isothermal distribution of a dark matter does not depend on black hole baryonic mass but depends strongly on the dark matter density determined by factor $`\rho _0`$. If $`\rho _0`$ is anomaly big: $`\rho _0>G^1(c/8\pi )^2`$ the initial isothermal distribution is unstable: any small seed black hole will result in absorption of all dark matter halo. At $`\rho _0<G^1(c/8\pi )^2`$ distribution is steady. The absorbed dark matter mass appears proportional to baryonic mass and is determined by expression (20).

For our Galaxy $`Q0.01`$ and $`M_d`$ is less that 1%, what means, that the accretion of dark matter on the initial black hole $`M_b`$ lead to a very small growth of its mass only.

Actually, distribution function of dark matter (6) differs essentially from isothermal one. In this case total mass of particles with angular momentum smaller than $`m_g`$ is determined by expression

$$N_g=(2\pi )^3_0^{\mathrm{}}𝑑I_R_0^{\mathrm{}}𝑑m_m^m𝑑m_zf_0I_R^{1/8}\delta \left(m^2l_0^2I_R^2\right)\theta \left(m_gm\right),$$

hence

$$N_g=(2\pi )^3\frac{8}{9}f_0\left(\frac{m_g}{l_0}\right)^{9/8}.$$
(21)
The equation (15) takes the form:

$$\frac{M_d}{M_b}=Q\left(1+\frac{M_d}{M_b}\right)^{9/8},$$
(22)
where $`Q=Q^{}M_b^{1/8}`$, $`Q^{}=(2\pi )^3\frac{8}{9}f_0\left(\frac{4G}{l_0c}\right)^{9/8}`$ and $`f_0`$ can be calculated by expression

$$M_H=(2\pi )^3_0^{I_{max}}𝑑I_R_0^{\mathrm{}}𝑑m^2f_0I_R^{1/8}\delta \left(m^2l_0^2I_R^2\right),$$

where the radial action of particles near the halo’s boundary is related to the halo’s full mass $`M_H`$ and radius $`R_H`$ by expression $`I_{max}\frac{1}{\pi }G^{1/2}M_H^{1/2}R_H^{1/2}`$, hense

$$f_0=\frac{9}{8}\frac{\pi ^{9/8}}{(2\pi )^3}\frac{M_H^{7/16}}{(GR_H)^{9/16}},$$
(23)
Using this relation, one can show that in our Galaxy the value of parameter $`Q`$ is again small (of the order of $`10^2`$), and mass of dark matter captured by black hole is not significant.

We would like to pay attention that investigation of this problem in previous papers  results to opposite statement: the initial black hole was growing without any limits. The reason of this discrepancy is in the fact that authors did not take into account the changing of dark matter particles distribution function caused by black hole absorption. That is why the loss cone was filled constantly and the flux on black hole was not exhausted. Considering the transformation of distribution function one comes to our result.

Thus, dark matter particles cannot cause a noticeable growth of the black hole because of conservation of their angular momentum in the spherically symmetric potential. As it follows from the remark in the end of Section 2, weak asymmetry of the potential does not change significantly this conclusion. The simplest process that can fill the loss cone is collisions of dark matter particles with stars. We shall discuss this point in next sections.

## 4  Kinetic equation

Gravitational interaction of dark matter particles and their collisions with stars is described by distribution function $`f(𝐫,𝐯,t)`$, which satisfies the kinetic equation :

$$\frac{f}{t}+\{H_0,f\}=St[f],$$
(24)
where $`H_0`$ is Hamilton function corresponding to the motion of particles in the averaged self-consistent potential, $`St[f]`$ is the collision term. Gravitational interaction of dark matter particles with individual stars is described by Coulomb law, hence the collision term may be written in Landau form

$$St[f]=\frac{}{v_k}W_{kp}\frac{}{v_p}f.$$
(25)
Here $`W_{kp}=2\pi G^2M_{}\mathrm{\Lambda }d^3v^{}w_{kp}F(𝐯^{},r),w_{kp}=\left(u^2\delta _{kp}u_ku_p\right)/u^3`$, $`v_k`$ are the dark matter velocity $`𝐯`$ components, $`F(𝐯^{},r)`$ is the stars distribution function, $`𝐮=𝐯^{}𝐯`$ is the relative velocity of a star and dark matter particle, $`\mathrm{\Lambda }10`$ is the gravitational Coulomb logarithm.

In fact, as we have already noticed, the frequency of collisions of particles with stars is much less than the frequency of their orbital motion. Thus the kinetic equation (24) may be significantly simplified. For this purpose, let us rewrite it in the action-angle variables $`I,\varphi `$:

$$\frac{f}{t}+\omega _k\frac{f}{\varphi _k}=St_{I,\varphi }[f]$$
(26)
Taking into account that the initial distribution function (6) does not depend on angle variables $`\varphi `$, we note that its time variation is determined by the collision term (25) only. However, the collisions are rare and $`St[f]\nu f`$, where $`\nu `$ is the collision frequency. Hence we search a solution of the equation (26) in the form $`f=f_0+\nu f_1`$, where $`\nu f_1`$ is a small correction to $`f_0`$ , and both the terms $`f_0`$ $`f_1`$ depend on ”fast” $`t_0=t`$ and ”slow” $`t_1=\nu t`$ times. In zero approximation by the parameter $`\nu `$ we obtain

$$\frac{f_0}{t_0}+\omega _k\frac{f_0}{\varphi _k}=0$$

From this equation follows that the main part of distribution function $`f_0`$ does not depend on angle variables and ”fast” time. The first approximation gives

$$\nu \left(\frac{f_0}{t_1}+\frac{f_1}{t_0}+\omega _k\frac{f_1}{\varphi _k}\right)=St_{I,\varphi }[f_0].$$

Averaging this equation over the angle variables and ”fast” time, taking into account that the second and third terms become zero, we obtain equation

$$\frac{f(I,t)}{t}=\overline{St}[f(I,t)]$$
(27)
(Here we omit the index ”$`0`$” in distribution function and returned to usual time $`t`$.)

So, under assumptions on collision frequency $`\nu `$ made above, it is possible to find the distribution depending on action variables only, as a solution of kinetic equation (27) with time-averaged collision term:

$$\overline{St}[f]=\frac{1}{(2\pi )^3}d^3\varphi St_{I,\varphi }[f]$$
(28)
The time-independent kinetic equation in the averaged form (27) describing the distribution of stars in the vicinity of a giant black hole was first studied in ,  and . A purely Coulomb potential field was considered in these works.

Our aim is to examine time-dependent solutions of the kinetic equation (27) describing dynamics of particles in arbitrary central symmetric field. First, it is necessary to obtain the averaged collision integral (28). To do this we will follow the technique, developed in .

Let us consider tensor differential form (25) formally in six-dimensional phase space $`X=(𝐯,𝐱)`$, assuming coordinate components of the tensor $`W`$ being zero:

$$\frac{}{v_k}W_{kp}\frac{}{v_k}=\frac{}{X_\mu }W_{\mu \nu }\frac{}{X_\nu };k,p=1,2,3;\mu ,\nu =1,\mathrm{},6$$

In this six-dimensional space we make canonical transformation to the action-angle variables

$$X_\mu Y_\mu ^{};Y=\{I_R,m,m_z;\varphi _1,\varphi _2,\varphi _3\}.$$

The differential form then takes the form

$$\frac{}{X_\mu }W_{\mu \nu }\frac{}{X_\nu }=\frac{1}{\sqrt{g}}\frac{}{Y_\mu ^{}}\sqrt{g}R_{\mu ^{}\nu ^{}}\frac{}{Y_\nu ^{}},R_{\mu ^{}\nu ^{}}=\frac{Y_\mu ^{}}{X_\mu }\frac{Y_\nu ^{}}{X_\nu }W_{\mu ,\nu }$$

where $`\sqrt{g}`$ is the Jacobian. However, Jacobian of canonical transformation is unity, hence the differential form in new variables preserves in the form:

$$\frac{}{X_\mu }W_{\mu \nu }\frac{}{X_\nu }=\frac{}{Y_\mu ^{}}R_{\mu ^{}\nu ^{}}\frac{}{Y_\nu ^{}}.$$

We search for the solution of kinetic equation (27) which does not depends on angle variables $`\varphi `$, hence $`\frac{}{Y}R\frac{}{\varphi }=0`$. After averaging over $`\varphi `$ terms $`\frac{}{\varphi }R\frac{}{I}`$ is vanished, and the expression for averaged collision integral becomes again three-dimensional:

$$\overline{St}[f]=\frac{}{I_k^{}}\overline{R}_{k^{}p^{}}\frac{}{I_p^{}}f,$$
(29)
where

$$\overline{R}_{k^{}p^{}}=\frac{1}{(2\pi )^3}d^3\varphi \frac{I_k^{}}{v_k}\frac{I_p^{}}{v_p}W_{kp}.$$
(30)
To give concrete expression to the differential (29), one should know the distribution function of stars $`F`$. The observations of stellar dynamics in the bulge indicate that in the first approximation the distribution function of stars could be assumed isotropic, i.e. depending on energy $`E^{}=\frac{v^2}{2}+\mathrm{\Psi }(r)`$ only. For simplicity we also assume that the dependence of stellar distribution function on energy is power-like:

$$F(v^{},r)=F_0E^\beta $$
(31)
One can show that the potential produced by stars distributed in accordance with (31) is also power function of $`r`$:

$$\mathrm{\Psi }(r)=\mathrm{\Psi }_0r^\alpha ,\alpha =4/(2\beta 1)$$
(32)
We note that both distribution function (31) and potential (32) can be determined if the velocity dispersion of stars $`\sigma =\sqrt{v_k^2}=\sqrt{v^2/3}`$ is known. The dispersion could be choosen as a function of distance to the galaxy center:

$$\sigma (r)=\sigma _0r^{\alpha /2}.$$
(33)
The distribution function (31) and the potential (32) are uniquely dependent functions of the velocity dispersion, and hence the parameters $`\sigma _0,\alpha `$ only:

$$\mathrm{\Psi }_0=3\sigma _0^2\frac{_0^{\mathrm{}}x^2\left(1+x^2/2\right)^\beta 𝑑x}{_0^{\mathrm{}}x^4\left(1+x^2/2\right)^\beta 𝑑x},F_0=\frac{\alpha (1+\alpha )}{(4\pi )^2G_0^{\mathrm{}}x^2\left(1+x^2/2\right)^\beta 𝑑x}\mathrm{\Psi }_0^{2/\alpha }$$

Note that in the case $`\alpha =0`$ the relations formally lose their mathematical meaning. However, one can show that when $`\alpha 0`$ the stellar distribution function becomes isothermal:

$$F(E^{})=\frac{\rho _0}{\left(2\pi \sigma ^2\right)^{3/2}}e^{\frac{E^{}}{\sigma ^2}},\rho _0=\frac{\sigma ^2}{2\pi G},$$

$$\mathrm{\Psi }(r)=2\sigma ^2\mathrm{ln}(r),$$

and the velocity $`\sigma `$ dispersion does not depend on distance. As a rule, the parameter $`\alpha `$ in a bulge is rather small.

Calculations show that for isotropic stellar distribution function $`F(E^{})`$ tensor $`W_{kp}`$ in (25) takes the form

$$W_{kp}=A(E,r)\delta _{kp}B(E,r)\frac{v_kv_p}{v^2},$$
(34)
where

$$A=\frac{16\pi ^2}{3}G^2M_{}\mathrm{\Lambda }_{\mathrm{\Psi }(r)}^{\mathrm{}}𝑑E^{}F(E^{})\{\begin{array}{ccc}\hfill 1& ,& E<E^{}\hfill \\ \hfill \frac{3}{2}\frac{v^{}}{v}\left(1\frac{v^2}{3v^2}\right)& ,& E>E^{}\hfill \end{array}$$

$$AB=\frac{16\pi ^2}{3}G^2M_{}\mathrm{\Lambda }_{\mathrm{\Psi }(r)}^{\mathrm{}}𝑑E^{}F(E^{})\{\begin{array}{ccc}\hfill 1& ,& E<E^{}\hfill \\ \hfill \frac{v^3}{v^3}& ,& E>E^{}\hfill \end{array}$$

We notes, that the radial action $`I_R`$ of particles moving in the potential (32) could be with good accuracy approximated by the expression

$$I_R(E,m)J(E)b_\alpha m,$$
(35)
where $`b_\alpha 1`$ is some positive constant,

$$J(E)=\frac{\sqrt{2}_0^1\left(1x^\alpha \right)^{1/2}𝑑x}{\pi \mathrm{\Psi }_0^{1/\alpha }}E^{\frac{1}{\alpha }+\frac{1}{2}},\alpha >0,$$

$$J(E)=\frac{\sigma }{\sqrt{\pi }}e^{\frac{E}{2\sigma ^2}},\alpha =0.$$

(We notice that in the case of Coulomb or oscillator potentials the equality (35) becomes exact and the constant $`b=1`$. In the case of isothermal stellar distribution $`b0.6`$)

The initial function (6) does not depend on the variable $`m_z`$, hence one can search for solution of the equation (29) as a function of $`I_R`$ and $`m`$. Calculating the coefficients of quadratic form (30) and taking into account (34), we rewrite the collision term in a form

$$\overline{St}[f]=\frac{1}{m}\frac{}{m}m\left(\overline{R}_{22}\frac{f}{m}+\overline{R}_{12}\frac{f}{I_R}\right)+\frac{}{I_R}\left(\overline{R}_{12}\frac{f}{m}+\overline{R}_{11}\frac{f}{I_R}\right),$$
(36)
where

$$\overline{R}_{11}=p2b_\alpha q+b_\alpha ^2s,$$

$$\overline{R}_{12}=qb_\alpha s,\overline{R}_{22}=s,$$

$$p=\left(\frac{dJ}{dE}\right)^2(AB)v^2_\varphi ,q=\left(\frac{dJ}{dE}\right)(AB)m_\varphi ,$$

$$s=Ar^2\frac{B}{v^2}m^2_\varphi ,$$

where $`<\mathrm{}>_\varphi =\frac{2}{T(E,m)}_r_{}^{r_+}\frac{dr}{v_r}(\mathrm{}),`$ and $`T=2_r_{}^{r_+}\frac{dr}{v_r}`$ is the period of dark matters partical radial oscillations, $`v_r=\sqrt{2}\left(E\mathrm{\Psi }(r)\frac{m^2}{2r^2}\right)^{1/2}`$ is radial velocity. Taking into account that angular momentum of dark matter particles are small, in the first approximation we shall calculate the coefficients $`\overline{R}_{ab}(I_R,m)`$ in (36) at $`m=0`$. Further, from the initial distribution function (6) one can see that for sufficiently large time, as far as the distribution function has sharp maximum at $`m=l_0I_R,l_01`$, one can neglect the diffusion by $`I_R`$ compared to diffusion by $`m`$. Summarizing, we rewrite the kinetic equation (27) as the diffusion equation depending on $`I_R`$ only parametrically:

$$\frac{f(I_R,m,t)}{t}=R(I_R)\frac{1}{m}\frac{}{m}m\frac{}{m}f(I_R,m,t)$$
(37)
Here the diffusion coefficient $`R=\overline{R}_{22}|_{m=0}=Ar^2_\varphi `$ calculated with account of (33), takes the form

$$R(I_R)=0.46GM_{}\mathrm{\Lambda }\sigma _0^{\frac{2}{2+\alpha }}I_R^{\frac{\alpha }{2+\alpha }}$$
(38)
We note that in the case of isothermal distribution of stars in the bulge $`\alpha =0`$, and the diffusion coefficient $`R`$ does not depend on $`I_R`$. We emphasize that the fact that $`I_R`$ comes into the final equation only as parameter is conditioned by the particular form of initial distribution function (6) containing small parameter $`l_0`$.

## 5 Flow of dark matter on a black hole.

As a result of direct capture of particles in the loss cone of seed black hole, the distribution function of dark matter particles takes the form

$$f(I_R,m,0)=f_0I_R^{1/8}\delta \left(m^2l_0^2I_R^2\right)\theta \left(mm_g\right).$$
(39)
It differs from (6) by factor $`\theta \left(mm_g\right)`$.

The diffusion of dark matter particles in angular momentum space, their drift into the loss cone (14) and further capture by black hole can be described by the diffusion equation (37) with initial condition (39) and boundary condition

$$f|_{m=m_g}=0.$$
(40)
The solution of the diffusion equation can be presented in a form

$$f(I_R,m,t)=_{m_g}^{\mathrm{}}𝑑m_1G(I_R,m,m_1,t)f(I_R,m,0)$$
(41)
where $`G`$ is the Green function of the boundary problem (37)-(40)

$$G=_0^{\mathrm{}}𝑑\lambda m_1Z_\lambda (m_1,m_g)Z_\lambda (m,m_g)e^{\lambda R(I_R)t},$$

$`Z`$ is the orthogonal system of fundamental functions of the boundary problem (37)-(40)

$$Z_\lambda (m,m_g)=\frac{J_0\left(\sqrt{\lambda }m_g\right)N_0\left(\sqrt{\lambda }m\right)N_0\left(\sqrt{\lambda }m_g\right)J_0\left(\sqrt{\lambda }m\right)}{\left(J_0^2\left(\sqrt{\lambda }m_g\right)+N_0^2\left(\sqrt{\lambda }m_g\right)\right)^{1/2}},$$

$`J_0,N_0`$ are Bessel and Neumann functions of zeroth order.

Let us now calculate the flux of dark matter on a black hole. Let $`D`$ be the area in the phase space where $`I_R0,mm_g,mm_zm`$. The whole mass of dark matter is $`N(t)=_Dd^3Id^3\varphi f(I_R,m,t)`$. Now, from the condition of mass conservation it follows that the flux through the boundary is $`S=dN/dt`$. Taking into account (25) we obtain

$$S(t)=_Dd^3Id^3\varphi \frac{}{I_k}\overline{R}_{kp}\frac{}{I_p}f(I_R,m,t).$$

Let us transform now the integral over volume into an integral over surface according to Stokes formula. Taking into account that the only flow of matter is through the surface $`m=m_g`$ , we find that

$$S(t)=2(2\pi )^3_0^{\mathrm{}}𝑑I_RR(I_R)m_g\frac{}{m}f(I_R,m,t)|_{m=m_g}.$$
(42)
Substituting the solution of diffusion equation (41) into (42), we obtain

$$S(t)=2\frac{(2\pi )^3}{\pi }f_0R_\alpha _0^{\mathrm{}}𝑑I_RI_R^{\gamma +\frac{1}{8}}_0^{\mathrm{}}𝑑\lambda e^{\lambda R_\alpha I_R^\gamma t}\frac{Z_\lambda (l_0I_R,m_g)}{\left(J_0^2\left(\sqrt{\lambda }m_g\right)+N_0^2\left(\sqrt{\lambda }m_g\right)\right)^{1/2}},$$

where

$$\gamma =\alpha /(2+\alpha ),R_\alpha =0.46\mathrm{\Lambda }GM_{}\sigma _0^{\frac{2}{2+\alpha }}.$$
(43)
It is convenient to introduce a normalized time $`T=l_0^\gamma R_\alpha t`$ and make a change of variables

$$\lambda \eta =\lambda T^{\frac{2}{2\gamma }},I_Ry=l_0I_RT^{\frac{1}{2\gamma }},m_gx=m_gT^{\frac{1}{2\gamma }}$$

The expression for the flux takes then a form

$$S(t)=2\frac{(2\pi )^3}{\pi }\frac{f_0}{l_0^{\frac{9}{8}+\gamma }}\frac{R_\alpha }{T^\zeta }\mathrm{\Phi }_\gamma (x),$$
(44)
$$\zeta =\frac{78\gamma }{168\gamma },\mathrm{\Phi }_\gamma (x)=_x^{\mathrm{}}𝑑yH_\gamma (x,y),$$

$$H_\gamma (x,y)=y^{\frac{1}{8}+\gamma }_0^{\mathrm{}}𝑑\eta e^{\eta y^\gamma }\frac{J_0\left(\sqrt{\eta }x\right)N_0\left(\sqrt{\eta }y\right)N_0\left(\sqrt{\eta }x\right)J_0\left(\sqrt{\eta }y\right)}{J_0^2\left(\sqrt{\eta }x\right)+N_0^2\left(\sqrt{\eta }x\right)}.$$

One can show that $`H_\gamma (x,y)`$ is finite positive function, it is equal to zero at $`y=x`$ and exponentially small at $`y>x+4`$. From this follows that at the moment $`t`$ the flow of dark matter to the black hole is made up by particles with radial action range in an interval

$$\frac{m_g}{l_0}<I_R<\frac{m_g}{l_0}+\frac{4}{l_0}T^{1/(2\gamma )}.$$
(45)
Thus, the area from where the dark matter flows into black hole, grows with time as $`t^{1/(2\gamma )}`$. For function $`\mathrm{\Phi }`$ in expression (45), there is an estimate $`\mathrm{\Phi }_\gamma (x)x^{\frac{1}{8}+\frac{\gamma }{2}}.`$ Finally, we find the expression for the flux of dark matter onto a black hole through the boundary of the loss cone $`m=m_g`$:

$$S(t)=2\frac{(2\pi )^3}{\pi }\frac{f_0}{l_0^{\frac{9}{8}+\gamma }}R_\alpha ^{1/2}\frac{m_g^{\frac{1}{8}+\frac{\gamma }{2}}}{t^{1/2}}.$$
(46)
The expression (46) has been obtained in assumption that the boundary of the loss cone $`m_g`$ does not depend on time. In reality, the boundary moves in accordance with the whole black hole mass. However, from (46) follows that the flux depends only weakly on the value $`m_g`$. It means, that the growing of $`m_g`$ is a quasistationary process and one can use the expression (46), assuming $`t`$ as a parameter: $`m_g(t)=4GM_{bh}(t)/c`$. Hence, the growth of the black hole is described by

$$\frac{dM_{bh}}{dt}=2\frac{(2\pi )^3}{\pi }\left(\frac{4G}{c}\right)^{\frac{1}{8}+\frac{\gamma }{2}}\frac{f_0}{l_0^{\frac{9}{8}+\gamma }}R_\alpha ^{1/2}\frac{M_{bh}^{\frac{1}{8}+\frac{\gamma }{2}}}{t^{1/2}}.$$

Assuming the seed black hole mass small, we obtain the solution of this equation

$$M_{bh}=Ct^{\frac{4}{74\gamma }},$$
(47)
$$C=\left(\left(\frac{74\gamma }{2}\right)\frac{(2\pi )^3}{\pi }\left(\frac{4G}{c}\right)^{\frac{1}{8}+\frac{\gamma }{2}}\frac{f_0}{l_0^{\frac{9}{8}+\gamma }}R_\alpha ^{1/2}\right)^{\frac{8}{74\gamma }}.$$

To estimate the parameter $`f_0`$ one can use the relation (23). Thus, the growth of black hole at the expense of absorption of dark matter scattered by stars follows a power law $`M_{bh}t^a,a4/7`$.

## 6 Conclusion.

In conclusion, we compare roughly the results of the theory presented above with the observational data. Fig.2 represents the dependence of stellar velocity dispersion on distance to the galactic center for our Galaxy and galaxies M31, NGC4258 . One can see that the distribution of stars in bulges of M31 and NGC4258 is close to isothermal, i.e. up to the area of action of black hole ($`r<(47)`$pc) the velocity dispersion is almost constant, it does not depend on $`r`$: $`\sigma =\sigma _0200`$km/s. In the case of isothermal distribution of stars in the bulge we can rewrite the expression (47):

$$M_{bh}=810^7M_{}\left(\frac{l_0}{0.1}\right)^{\frac{9}{7}}\left(\frac{M_H}{10^{12}M_{}}\right)^{\frac{1}{2}}\left(\frac{R_H}{100kps}\right)^{\frac{9}{14}}\left(\frac{\sigma }{200km/s}\right)^{\frac{4}{7}}\left(\frac{t}{310^{17}s}\right)^{\frac{4}{7}}$$
(48)
Let us suppose that parameter $`l_0`$ of the dark matter haloes around these galaxies is close to average value $`<l_0>=0.1`$. Implying, that masses $`M_H`$ and sizes $`R_H`$ of haloes are roughly equal and constitute $`10^{12}M_{}`$ and $`100`$kpc respectively we obtain from (48) the mass of the black hole $`M_{bh}810^7M_{}`$ (We assume, that the age $`t`$ of black hole is comparable with the age of the Universe $`t310^{17}`$s,). Thus, theoretical predictions of masses of black holes in these galaxies are close enough to the observed values: $`M_{bh}=(2.08.5)10^7M_{}`$ for M31 and $`M_{bh}=(3.84.0)10^7M_{}`$, for NGC4258.

The central part of our Galaxy’s bulge does not satisfy the assumption of isothermal distribution of stars at $`r100`$pc. At the same time, as it is clear from (45), the dark matter accumulates mainly from the central region. The dependence of stellar velocity dispersion on distance is shown at Fig. 3 . One can see that at the scales from $`10`$pc to $`100`$pc the dispersion can be roughly represented as $`\sigma =60`$km/s$`\left(\frac{r}{10nc}\right)^{0.3}`$. From (33) then follows $`\alpha 0.6`$. In the case, the formula (47) gives for the black hole mass a value $`210^7M_{}`$, which is much more than observed value $`M_{bh}2.610^6M_{}`$. To obtain more precise estimation of the black hole mass, more information on the stellar and dark matter distribution functions and their evolution is needed. Also, the contribution of the black hole into the joint gravitational potential may be important. However, we see that even the rough calculation of absorption of dark matter only may give a reasonable estimate of black hole mass. This estimate is consistent with observed masses of most giant black holes (see Fig.1). The fact may be considered as an additional prove of the general theory of the role of CDM in the large-scale structuring of the Universe .

Note also that the kinetic theory developed in Section 5 allows to describe the absorption of luminous (baryon) matter as well as dark matter. More precisely, it is possible to take into account the capture of stars by black hole because of their gravitational scattering. This process is especially important for giant black holes having mass of the order $`M_{bh}10^9M_{}`$ and active galactic nuclei.

Another important process described by the kinetic equation (27) is the possibility of diminishing of amount of dark matter in bulge. The change of dark matter particles’ energy as a result of gravitational scattering by stars leads to their ”heating” and pushing them out of bulge – the phenomenon analogous to Fermi acceleration. Therefore, the density of dark matter compared to that of luminous matter decreases. Qualitatively this corresponds to the results of recent observations discussed in the literature , .

Further development of the consistent kinetic theory discussed above and its detailed comparison with observational data would be undoubtedly useful for understanding of main physical processes in the galactic nuclei.

The authors are grateful to V.A. Sirota and M.I. Zelnikov for numerous and useful discussions.