{'fraction_non_alphanumeric': 0.038171394085697044, 'fraction_numerical': 0.010070911285455642, 'mean_word_length': 4.500622406639004, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stochastic Dynamics Toward the Steady State of Self-Gravitating Systems

We will construct a theory which can explain the dynamics toward the steady state self-gravitating systems (SGSs) where many particles interact via the gravitational force. Real examples of SGS in the universe are globular clusters and galaxies. The idea is to represent an interaction by which a particle of the system is affected from the others by a special random force. That is, we will use a special Langevin equation, just as the normal Langevin equation can unveil the dynamics toward the steady state described by the Maxwell-Boltzmann distribution. However, we cannot introduce the randomness into the system without any evidence. Then, we must confirm that each orbit is random indeed. Of course, it is impossible to understand orbits of stars in globular clusters from observations. Thus we use numerical simulations. From the numerical simulations of SGS, grounds that we use the random noise become clear. The special Langevin equation includes the additive and the multiplicative noise. By using the random process, we derive the non-Maxwellian distribution of SGS especially around the core. The number density can be obtained through the steady state solution of the Fokker-Planck equation corresponding to the random process. We exhibit that the number density becomes equal to the density profiles around the core by adjusting the friction coefficient and the intensity of the multiplicative noise. Moreover, we also show that our model can be applied in the system which has a heavier particle, corresponding to the black hole in a globular cluster.

A self-gravitating system (SGS) is a system where many particles interact via the gravitational force. When we shall explain a distribution of SGS in phase space, the Boltzmann-Gibbs statistical mechanics is not useful. This is because the statistical mechanics is constructed under the condition of the additivity of the energy: as well known, the total energy of several SGSs is not equal to the sum of the energy of each system. In fact, SGS does not have a tendency to become a state characterized with a temperature.
If we use the statistical mechanics assuming that the state of the SGS with an equal mass m becomes isothermal with the temperature T and that the particles of the system are distributed spherically symmetrically, what kind of distribution can be obtained? Then, the structure in phase space can be determined using the Maxwell-Boltzmann distribution.
For example, the number density at a radial distance r in real space is given by where Φ(r) is the mean gravitational potential per mass generated by this whole system at r and k B is the Boltzmann constant. This potential should satisfy a relation with number density by the Poisson equation, where G is the gravitational constant. A special solution to Eq. (1) and this Poisson equation is n MB (r) = k B T /2πGm 2 r 2 , known as the singular isothermal sphere [1]. However, this solution has two problems: infinite density at r = 0 and infinite total mass. Even though we solve the equations with a finite density at r = 0, the solutions behave ∝ r −2 at a large r, and so we cannot avoid the infinite total mass problem. In either case, the solutions are unrealistic.
Of course, real examples of SGS in the universe, e.g., globular clusters and galaxies, have various structures with a finite radius. As for most globular clusters, it is known that their number densities in real space have a flat core and behave as a power law outside this core.
King interpreted these profiles by introducing the new distribution function in phase space, known as the lowered Maxwellian; in which E is the total energy of a particle belonging to a globular cluster. This distribution becomes zero when the total energy is greater than mΦ t , and so Φ t can be understood as a potential energy per mass at the surface of the globular cluster. Because the velocity of the particle must be in the range of 0 ≤ v ≤ 2 {Φ t − Φ(r)}, the number density n KM (r) can be obtained integrating f (r, v) as Moreover, using a dimensionless potential W (r) ≡ −mβ {Φ t − Φ(r)} and integrating by parts, the number density becomes n KM (r) ∝ e W (r) As mentioned before, the potential energy and the number density has a relation through the Poisson equation. Thus, W (r) must satisfy the following equation: where a ≡ 9/{4πGm 2 βn KM (0)} corresponds to the core radius. The number density satisfying Eqs. (5) and (6) can be calculated numerically as shown in Fig. 1. This is called the King model [2]. When W (0) is larger than about 5, the number density around the origin can be represented by the following approximation: which is shown as the red curve in Fig. 1 Since King put forward this model, this number density has been applied to fitting for the surface brightness of many globular clusters, for example as in Ref. [3][4][5][6][7], that is, most exponents of power law outside the core of globular clusters are similar to −3 which cannot be explained by the model with the isothermal assumption. But, it is not easy to see what kind of dynamics occurred in the system, because his procedure was done to the distribution function in the steady state.
So, we will construct a theory which can explain the dynamics toward such a special steady state described by the King model especially around the origin. The idea is to represent an interaction by which a particle of the system is affected from the others by a special random force described by a position-dependent intensity noise, in other words the multiplicative noise, that originates from a fluctuation only in SGS. That is, we will use a special Langevin equation, just as the normal Langevin equation with a constant-intensity noise, in other words the additive noise, can unveil the dynamics toward the steady state described by the Maxwell-Boltzmann distribution. However, we cannot introduce the randomness into the system without any evidence. Then, we must confirm that each orbit is random indeed. Of course, it is impossible to understand orbits of stars in globular clusters from observations. Thus, we must use numerical simulations.
From the numerical simulations of SGS, the ground that we can use the random noise becomes clear. The special Langevin equation includes additive and multiplicative noises.
By using this stochastic process, we derive that non-Maxwell-Boltzmann distribution of SGS especially around the origin. The number density can be obtained through the steady state solution of the Fokker-Planck equation corresponding to the stochastic process. We exhibit that the number density becomes equal to the density profiles around the origin, Eq. (7), by adjusting the friction coefficient and the intensity of the multiplicative noise.
Moreover, we also show that our model can be applied in the system which has a heavier particle (5-10 times as heavy as the surrounding particle). The effect of the heavier particle in SGS, corresponding to a black hole in a globular cluster, has been studied for long time. If the black hole is much heavier than other stars, a cusp of the density distribution appears at the center of a cluster [8][9][10]. The observations which suggest that intermediate mass black hole (IMBH, ∼ 10 2 − 10 3 M ⊙ ) is in the globular cluster in recent years are accomplished one after another [12][13][14][15][16][17]. Although these studies are very interesting, our model does not treat these situations: in our model, the heavier particle is too lighter than IMBH. Our model corresponds small globular cluster (10 4 stars) with only a stellar black hole (≃ 1 − 10M ⊙ ).
Here, note that we have reported similar results in our previous letter [11]. In this paper, however, we demonstrate how we executed our numerical simulations. Moreover, a treatment for stochastic differential equations becomes precise, and so the analytical result derived by a different method changes a little. is exerted on particles in SGS and the gravitational potential is asymptotically flat, we must compute interaction of all particle pairs. When we treat N particles, the computation of interaction becomes O(N 2 ) by direct approach. By these reasons, we require huge computation for numerical simulation of the evolution of SGS.
For time evolution of SGS, many improvements of algorithm and hardware have been carried out. First, we consider integrator for simulation. For long-time evolution, both the local truncation error and the global truncation error are noticed. These error occur deviation of the conservation physical quantities such as total energy. For compression of the global truncation error, symplectic integrator has been developed. The symplectic integrator conserves the total energy for long-time evolution. We apply 6th-order symplectic integrator for the time evolution of SGS. Secondly, we apply special-purpose processor for the computation of the interaction. Most of the computation of the time evolution in SGS is 2-body interaction. As special-purpose processors, GRAPE system has been developed [18].
GRAPE system can compute 2-body interaction from position and mass of particles quickly.
In our study, we apply GRAPE-7 chip, which consists of Field-Programmable Gate Array (FPGA) for computation of the interaction [19]. GRAPE-7 chip implements GRAPE-5 compatible pipelines [35]. The performance of GRAPE-7 chip is approximately 100 GFLOPS and is almost equal to a processor of present supercomputers, but the energy consumption of the chip is only 3 Watts. Using sophisticated integrator and special-purpose processor, we have analyzed time evolution of SGS.

B. Symplectic integrator
For time evolution, we must choose reasonable integrator for simulation. For long-time evolution, not only the local truncation error but also the global truncation error is noticed.
For example, 4th-order Runge-Kutta method has been applied for time evolution of physical systems [21]. Although its local truncation error is O ((∆t) 5 ), because its error accumulates, the global truncation error increases during time evolution. For example, we apply 4th-order Runge-Kutta method for the harmonic oscillation.
Using exact solutions, the orbit of the harmonic oscillation in the phase space draws a circle.
Of course, the total energy is conserved. On the other hand, when we apply 4th-order Runge-Kutta method for time evolution, the total energy is decreased monotonically.
The orbit of the harmonic oscillation in the phase space draws a spiral and it converges to origin (p = q = 0). When we consider long-time evolution, 4th-order Runge-Kutta method is not reasonable integrator. If Hamiltonian of physical system is given, we can apply symplectic integrator which based on canonical transformation [22][23][24]26]. This method suppresses increase of the global truncation error. In generic case, the symplectic integrator is implicit method. If Hamiltonian is divided to coordinate parts and momentum parts, the integrator becomes explicit method. The procedure of low-order integration becomes easy more than Runge-Kutta method. The simplest integrator is called "leap-flog method" (2nd-order integrator).
Using leap-flog method for the harmonic oscillator, the following equation is satisfied.
Therefore the orbit in the phase space draws an oval. The deviation from the exact solution is suppressed. To suppress the local truncation error, higher-order symplectic integrators have been developed. We apply 6th-order symplectic integrator for time evolution of SGS [25].  The symplectic integrator conserves the total energy and the symplectic structure in generic cases. When we use n-th order symplectic integrator, the local truncation error of the total energy becomes O ((∆t) n+1 ). Furthermore, the global truncation error is not accumulated [27].
In SGS, because the interaction at zero distance diverges, the local truncation error would diverge in long-time evolution. For avoidance of this divergence, some kind of softening parameter has been introduced to gravitational interaction. When the nature of the pure gravity is analyzed, the regularization procedure of interaction is required [28,29].

C. Steady number density in numerical simulation
Now, we investigate the steady number density (SND) of the SGS with a mass m including a particle with a mass M by numerical simulation. In particular, we show that SNDs have a core and behave as a power law outside the core.
The system is composed of N = 10000 particles. At t = 0, all velocities of the particles are zero and they are distributed by n 0 (r) ∝ (1 + r 2 /a p 2 ) −5/2 (0 ≤ r ≤ 4a p ), which is the density in real space of Plummer's solution [1]. In this SGS, we put another particle with a mass M in the origin at t = 0. We shall change the mass as M/m = 1, 5, and 10.
Throughout this paper, we adopt a unit system where the core radius of Plummer's solution a p , the initial free fall time t f f , and the total mass N · m are unity.
We started the numerical simulation under these conditions. For dynamical evolution, we used GRAPE-7 at Center for Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan. For the computation of gravitational force, we applied Plummer's softening: the potential energy between the ith and jth particles separated by a distance where ε s is the softening parameter. We set ε s = 10 −3 . For time evolution, we used 6th-order symplectic integrator [25]. The time step for the simulation is defined as ∆t = 10 −5 . We carried out simulations until t = 100 t f f . During simulations, the error in total energy was less than 0.1%.
First, most particles collapse into the origin within several t f f . After approximately 20 t f f , the distribution becomes stable and the system reaches the steady state. Of course, we can confirm whether the system becomes steady or not from the profile of the number density. However, furthermore we also focus on the number of particles inside a sphere. at the origin C increases as M increases, which is simply understood to be a result of many particles being attracted by the heavier particle.

III. SIMPLE MODEL
A. Forces acting on each particle of SGS As shown in the last section, SND is the King-like profile even though the system includes the heavier particle. In this section, in order to explain these results and derive this non-Maxwell-Boltzmann distribution around the origin, we demonstrate a simple model based on stochastic process, which is quite different from the King model.
The reason why stochastic process appears in the SGS is as follows. After the collapse, the density around the origin becomes high. Thus, the particles around the region disturb the orbits of other particles repeatedly, so that their movements become random [36]. As the time at which this disturbance occurs, we introduce the local 2-body relaxation time t rel [30]: where σ(r) is the standard deviation of the velocity at r; we adopted ln(1/ε s ) as the Coulomb logarithm. Figure 4 shows the logarithm of t rel , which is calculated using the σ(r) and n(r) during the steady state obtained from our numerical simulation, as a function of log r. As expected, t rel around the origin is short. Our simulation continues after the collapse at about 80 t f f , which is sufficiently longer than the t rel around the core. As radius increases, however, t rel becomes longer than the rest of our simulation time, which means that no stochastic motion occurs at a large r. Therefore, note that our model is valid only in the neighborhood of the core. When constructing our model, the following points are premised: the model describes the stochastic dynamics near the steady state and the mean distribution is spherically symmetric.
As is well known, the gravitational force at r arising from such a spherically symmetric system depends only on the particles existing inside a sphere with a radius r, and this attractive force acts along the radial direction. In other words, this mean force −F (r) is the gradient of the mean potential: −F (r) = −m∂Φ(r)/∂r (< 0). Indeed, lim r→0 F (r) = 0.
Hence, F (r) can be expanded around the origin as It will be clarified later that the lowest exponent must be 1 and the coefficient c 0 is related to the number density at the origin C as For M/m = 1, we can identify another particle together with the other particles. Contrary to this, we must consider the effect of the particle in the case M/m = 1. Now, we suppose that the heavier particle exists at the origin. Then, the attractive force by this particle at r is −F (r) = −GmM/r 2 . We can estimate F (r) around this region as F (r) ∼ c 0 r = 4πGm 2 Cr/3, where we used Eq. (18). Thus, F (r)/F (r) ∼ 3Mr −3 /4πCm. This ratio becomes significant when r < ∼ 10 −2 , since C ∼ 10 6 as shown in Tab. II. Therefore, if r is smaller than the radius, particles are influenced by not only F (r) but also F (r), so that the core disappears. In fact, we have performed a numerical simulation with the heavier particle fixed at the origin, where this result is confirmed. On the other hand, Miocchi improved the King model in order to describe the steady state of a globular cluster including an IMBH and reported that the density becomes cuspy as the mass of the black hole increases [31]. The nature of a globular cluster when a massive black hole is much heavier than the surrounding star, have been studied as mentioned in Introduction. In this case, the massive black hole stays at the center mostly. Then, a cusp of the density distribution at the center appears.
Because the heavier particles in our numerical simulation are not very heavy, the particles are not trapped at the origin. Therefore, we do not consider the effect of heavier particles explicitly and we suppose that the particles influence SGS through the density at the origin C: as M becomes larger, it attracts more particles and C increases, as shown in Tab. II.
Thus, c 0 is an increasing function of M.
It is natural to consider that the distribution fluctuates around the mean because of the many disturbances. In fact, as shown in Fig. 2, the number of particles existing inside the sphere with a radius 1 fluctuates around the mean value. The fluctuating part of the distribution should not be spherically symmetric, so that this produces forces along not only the radial direction, but also other directions. We assume that they are random forces and set their intensity at r 2H(r). In addition to such random forces resulting from the fluctuating distribution, a particle at r is expected to be influenced by random forces generated from neighboring particles. We set the intensity 2D, which is independent of position.

B. Langevin equations
Stochastic dynamics under the above assumptions is described by the following Langevin equations in spherical coordinates: the radial direction the elevation direction and the azimuth direction where a r , a θ , and a φ are accelerations along those directions; γ is the coefficient of dynamical friction in the low velocity limit, independent of velocity [1]. In the Chandrasekhar dynamical friction formula, the coefficient is more complicated [1,32]. However, we use the coefficient in such a limit, because the density around the core is so high that particles around there move slowly. Now, we focus on the overdamped limit of these equations, because we have interests in the stochastic dynamics near the steady state. In the case of the normal Langevin equation with a constant-intensity noise, we only neglect the inertial term. But, as for special Langevin equations with noises whose intensity depends on a position, the new force −∇H(r)/2mγ should be considered additionally [37]. Thus, the Langevin equations in the overdamped limit becomes as follows: azimuth direction : mγr sin θφ = 2H(r)η φ (t) where the prime indicates a derivative with respect to r. The noises in each Langevin equation, ξ i (t) and η j (t) (i, j = r, θ, and φ), are zero-mean white Gaussian and correlate only with themselves. Indeed, the correlation function is the Dirac delta function [38].
Here, revisit the position-dependent intensity noise. We have introduced such a noise in order to represent a random force originating from the fluctuation of the distribution around the mean value which yields the mean force −F (r). Thus, the first and the second terms on the right-hand side of Eq. (22) must denote that the mean force acting along radial direction is fluctuating. As a minimal formulation describing this situation, we propose the following one: in which ǫ is a positive constant. This can be realized by setting Note that this fluctuating mean force is the essential feature for SGS. Since the gravitational force is a long-range one, each particle is influenced from the whole system. The mean force is produced by the mean potential which is decided by the number density in the steady in which we have replaced H(r) by F (r) using Eq. (26). Then, the PDF with the Jacobian ρ(r, t) ≡ 4πr 2 P (r, t) satisfies the following Fokker-Planck equation.
This equation is useful when integrating with respect to r.
The steady state solution ρ st (r) satisfies Eq. (27) with the left-hand side zero. By integrating the equation with respect to r, we have Now, we impose the binary condition that P st (r) ≡ ρ st (r)/(4πr 2 ) and the derivative do not diverge at the origin. Then, when r → 0, and lim r→0 ρ ′ st (r) = lim r→0 4π(2rP st (r) + r 2 P ′ st (r)) = 0 , by which the constant on the right-hand side of Eq. (28) is decided and we obtain Thus, if F (r) is obtained, ρ st (r) can also be obtained. Here, F (r) relates with SND, n st (r), through the following relation, since −F (r) = −mΦ ′ (r).
Incidentally, the SND can be obtained by multiplying PDF in the steady state by total number N: Therefore, equation (32) can be represented as In short, ρ st (r) and F (r) are closely connected with each other through Eqs. (31) and (34). Here, we focus on the asymptotic behaviors of them around the origin, since our model is valid around there as mentioned before. Furthermore, due to this approach, we can treat them analytically.
Assume that F (r) can be expanded around the origin with the lowest exponent k as follows.
Substituting this expression into Eq. (34), we find that ρ st (r) can also be expanded like After substituting both Eqs. (35) and (36) into Eq.(31), we can obtain Firstly, we compare the lowest order terms on the both hand sides of Eq. (37), so that the following relation can be seen: Therefore, we can conclude that k = 1. Secondly, compare the next lowest order terms proportional to r 3 and we get and so c 1 = 0. Lastly, selecting only terms proportional to r 4 from Eq. (37), we can find from which the following relation can be obtained: Without going into detail, we can see that c 3 = 0 by comparison with terms containing r 5 .
So, ρ st (r) becomes as follows: Here, if we set [39] ρ st (r) can be expressed around the origin like which yields Thus, we can derive the number density around the origin of SGS from the model using stochastic dynamics.

IV. DISCUSSION
In this section, we investigate the results derived in the preceding section and understand the roles of two noises and the heavier particle in Eq. (45). Additionally, we discuss the difference between the King model and our model.
As in Eq. (43), the exponent κ must be larger than 1, which does not contradict our numerical simulation shown in Tab. II. In order for Eq. (45) to correspond completely to the King model, κ = 3/2 or γ = ǫc 0 /m = 4πGmǫC/3 must hold. We can regard this relation between the friction coefficient γ and the intensity of the multiplicative noise ǫ as a kind of fluctuation-dissipation relation [34], which usually plays an important role when a stochastic process with a constant-intensity noise goes to the equilibrium state described by the Maxwell-Boltzmann distribution.
The core radius a is proportional to a square root of the intensity of the additive noise D owing to Eq. (43). Then, this intensity spreads the region where the density is almost constant. This is recognized as the effect of this noise, which makes a system homogeneous and isothermal. The existence of the core at globular clusters shows that such a diffusive effect does not disappear even for the system with long-range force. In other words, all statistical mechanical features observed in a system with short-range force, that is, normal Therefore, our result goes to a singular isothermal sphere, as discussed at the beginning of this paper, by which the number density of globular clusters cannot be explained.
Here, we examine the relationship between the King model and our model. King transformed the distribution function in the phase space in order to avoid a singular isothermal sphere. In our model, we introduce the multiplicative noise into the system influenced by the mean force and the additive noise whose PDF becomes Maxwellian in the steady state, as shown in Eq. (47), so that the non-Maxwell-Boltzmann distribution (45) is derived. In short, although these procedures seem to be different, they may have the same meaning at least around the origin. However, we emphasize that the stochastic dynamics around there near the steady state becomes clear owing to Eqs. (19), (20), and (21).

V. CONCLUSION
In conclusion, the non-Maxwell-Boltzmann distribution (45) has been obtained using the stochastic dynamics with the fluctuating mean force and the additive white noise. This number density can be the same as that of the King model around the origin by controlling friction coefficient and the intensity of multiplicative noise. Furthermore, our model can describe the SGS having a heavier particle. Of course, these results are consistent with our numerical simulation. We can say that such a stochastic dynamics occurs behind the background of the King model. In short, the diffusive effect, which is represented by the additive noise, is universal even in SGS, and it is particular to SGS that the fluctuation of the distribution around the mean value producing the mean force makes influence on each particle of this system, which our simple model can describe.
Finally, note that our result is available only in the neighborhood of the origin. Therefore, we must derive the density globally by further extended model and investigate the difference between the model and the King model.