Dynamics to the universal structure of one-dimensional self-gravitating systems in the quasi-equilibrium state

TL;DR

This paper explores the universal structure of one-dimensional self-gravitating systems in quasi-equilibrium, revealing conditions for universality, proposing a phenomenological model, and explaining the system's size determination.

Contribution

It demonstrates that the null virial condition leads to universal density profiles and introduces a Langevin-based model to describe these structures.

Findings

Universal density profile similar to higher-dimensional systems

Null virial condition ensures universality

A mechanism for system size determination

Abstract

We investigate the quasi-equilibrium state of one-dimensional self-gravitating systems. If the null virial condition is satisfied at initial time, it is found that the number density around the center of the system at the quasi-equilibrium state has the universality similar to two- and three-dimensional self-gravitating systems reported in \cite{Tashiro16,Tashiro10}. The reason why the null virial condition is sufficient for the universality is unveiled by the envelope equation. We present a phenomenological model to describe the universal structure by using a special Langevin equation with a distinctive random noise to self-gravitating systems. Additionally, we unveil a mechanism which decides the radius of the system.

Equations

$$\mbox{$\mathcal{N}$}(r)\simeq\frac{\mbox{$\mathcal{N}$}(0)}{(1+r^{2}/a^{2})^{\kappa}}$$

$$\Delta (t)=\frac{3}{5E}\sqrt{\frac{1}{N}\sum _{i=1}^N{\left[\overline{{\mathcal{E}}_i(t)}-\frac{5E}{3}\right]}^2}$$

$$\frac{{\mathrm{d}}^2\varphi (x)}{\mathrm{d}{x}^2}=4\pi G\rho (x)$$

$$\frac{{\mathrm{d}}^2{\varphi }_{\mathrm{sp}}(x)}{\mathrm{d}{x}^2}=4\pi Gm\delta (x),$$

$${\varphi }_{\mathrm{sp}}(x)=2\pi Gm|x|,$$

$$H=\sum _{i=1}^N\frac{m{u}_i^2}{2}+2\pi G{m}^2\sum _{i>j}|x_i-x_j|,$$

$$f_n(r,v)\propto {\left[{\varphi }_n({\mathcal{R}}_n)-\left\{\frac{1}{2}{v}^2+{\varphi }_n(r)\right\}\right]}^{n-1/2}\Theta ({\mathcal{E}}_n)$$

$$\frac{{\mathrm{d}}^2{\varphi }_n(r)}{\mathrm{d}{r}^2}=4\pi Gm{\mathcal{N}}_n(r),$$

$${\ddot{r}}_e(t)=\frac{M^2\epsilon (t{)}^2}{4r_e(t{)}^3}-\frac{\Omega (t)}{r_e(t)},$$

$$\sum _{i=1}^N\Delta {\mathcal{E}}_i(t)=\frac{1}{m}\left\{\Omega (t)-\Omega (0)\right\},$$

$$\sum _{i=1}^N\Delta {\mathcal{E}}_i^{\mathrm{qe}}=-\frac{1}{3m}\Omega (0),$$

$$m\ddot{r}(t)+m\gamma \dot{r}(t)=-F(r)\left\{1+\sqrt{2ϵ}\eta (t)\right\}+\sqrt{2D}\xi (t),$$

$$m\gamma \dot{r}(t)=-F(r)\left\{1+\sqrt{2ϵ}\eta (t)\right\}-\frac{\partial }{\partial r}\frac{ϵ}{2m\gamma }F(r{)}^2+\sqrt{2D}\xi (t).$$

$$\frac{\partial }{\partial t}P(r,t)$$

$$\left\{\frac{D}{(m\gamma {)}^2}+\frac{ϵF(r{)}^2}{(m\gamma {)}^2}\right\}{P^{\mathrm{qe}}}^{\prime }(r)+\left[\frac{2ϵF(r)F^{\prime }(r)}{(m\gamma {)}^2}+\frac{F(r)}{m\gamma }\right]P^{\mathrm{qe}}(r)=\text{mbox}const..$$

$${{\mathcal{N}}^{\mathrm{qe}}}^{\prime }(r)=-\frac{F(r)\left\{m\gamma +2ϵF^{\prime }(r)\right\}}{\left\{D+ϵF(r{)}^2\right\}}{\mathcal{N}}^{\mathrm{qe}}(r).$$

$$F^{\prime }(r)=4\pi G{m}^2{\mathcal{N}}^{\mathrm{qe}}(r),$$

$$[\mathrm{length}]=\sqrt{3\mathcal{T}/2\pi G{m}^3{\mathcal{N}}_0^{\mathrm{qe}}}$$

$$[\mathrm{force}]=2\sqrt{6\pi Gm{\mathcal{N}}_0^{\mathrm{qe}}\mathcal{T}}$$

$${\bar{\mathcal{N}}}^{\mathrm{qe}}{}^{\prime }(\bar{r})=-\frac{6F(\bar{r})\left\{1+2q{F}^{\prime }(\bar{r})\right\}}{\left\{1+14q+6qF(\bar{r}{)}^2\right\}}{\bar{\mathcal{N}}}^{\mathrm{qe}}(\bar{r}),$$

$${\bar{F}}^{\prime }(\bar{r})={\bar{\mathcal{N}}}^{\mathrm{qe}}(\bar{r}),$$

$${\bar{\mathcal{N}}}^{\mathrm{qe}}(\bar{r})=\frac{1+14q-3f_q^{-1}(\bar{r}{)}^2}{1+14q+6q{f}_q^{-1}(\bar{r}{)}^2}$$

$$\bar{F}(\bar{r})=f_q^{-1}(\bar{r})$$

$$f_q(X)\equiv -2qX+\frac{(1+2q)(1+14q)}{\sqrt{3(1+14q)}}{\tanh }^{-1}\left(\sqrt{\frac{3}{1+14q}}X\right).$$

$$\kappa =\frac{3(1+2q)}{1+14q}.$$

$$\bar{r}=\frac{1}{\sqrt{3}}{\tanh }^{-1}\left[\sqrt{3}\bar{F}(\bar{r})\right].$$

$$\bar{F}(\bar{r})=\frac{1}{\sqrt{3}}\tanh \left[\sqrt{3}\bar{r}\right],$$

$${\bar{\mathcal{N}}}^{\mathrm{qe}}(\bar{r})=1-3{\bar{F}(\bar{r})}^2={\mathrm{sech}}^2\left[\sqrt{3}\bar{r}\right],$$

$$\Delta {\mathcal{E}}_{\mathrm{out}}^{\mathrm{qe}}=2\pi Gm\sum _{i=1}^N|x_i(0)|.$$

$${\mathcal{E}}_{\mathrm{out}}^{\mathrm{qe}}={\mathcal{E}}_{\mathrm{out}}(0)+\Delta {\mathcal{E}}_{\mathrm{out}}^{\mathrm{qe}}=2\pi Gm\sum _{i=1}^N|x_i(0)-x_{\mathrm{out}}(0)|+2\pi Gm\sum _{i=1}^N|x_i(0)|.$$

Full text

I Introduction

Self-gravitating systems (SGSs) exhibit various interesting properties which systems with short-range forces do not have Levin14 ; Campa09 . One of the examples is the presence of a stable state different from the equilibrium state. In this Letter, we refer to the state as a quasi-equilibrium state (qES). In general, the qES depends strongly on initial distributions. If a null virial condition is fulfilled at initial time (i.e., all the velocities are zero), however, it has been found that the number densities of the two-dimensional SGS (2DSGS) or the three-dimensional SGS (3DSGS) at qES around the center of the system can be universally approximated by

[TABLE]

with $\kappa\sim 1$ in 2DSGS Tashiro16 and $\kappa\sim 3/2$ in 3DSGS Tashiro10 regardless of the initial distribution in real space, where $r$ represents the distance from the center of the system. The results of 2DSGS and 3DSGS coincide with the observations of the number density of molecular clouds in IC 5146 or the Taurus molecular complex Arzoumanian11 ; Stepnik03 and the number density of stars in globular clusters Binney08 ; Peterson75 ; Trager95 , respectively. Thus, we can expect that if the null virial condition is satisfied, the value of $\kappa$ at qES is roughly the half of the dimension regardless of the initial distribution. It is natural that the question arises as to whether the same universality holds in one-dimensional self-gravitating systems (1DSGSs).

1DSGS with $N$ particles is the simplest SGS whose two-body interaction force is constant independent of the distance between the particles, which will be demonstrated later. The majority of studies about 1DSGS were aimed at investigating a relaxation to the thermal equilibrium state whose distribution function (DF) in phase space was obtained analytically in Ref. Rybicki71 . The relaxation time toward the state was originally estimated to be $N^{2}t_{d}$ Hohl67 ; Hohl67b ; Hohl68 where $t_{d}$ is the dynamical time which is the characteristic time for a particle to pass through the system. After the studies, however, the relaxation time was found to depend on initial conditions. If the initial distribution is a water-bag form, the relaxation time was revealed as much more than $N^{2}t_{d}$ Wright82 ; Luwel85 . On the other hand, for initial counterstreamed configuration with virial ratio 0.3, the system relaxes on the time scale $N^{2}t_{d}$ Luwel84 ; Severne84 ; Reidl91 . The relaxation process from the initial water-bag distribution was precisely described by the deviation from the equipartition law defined as

[TABLE]

where $E$ is the total energy of the system and $\overline{\mathcal{E}_{i}(t)}$ is the time averaged one particle energy per unit mass of $i$th particle until time $t$. It was unveiled that there are two different relaxation process, the microscopic relaxation occurring roughly after $Nt_{d}$ and the macroscopic relaxation occurring roughly after $4\times 10^{4}Nt_{d}$ Tsuchiya94 ; Tsuchiya96 .

The other studies about 1DSGS dealt with the phase space distribution at qES. When the initial water-bag distribution is in virial equilibrium (i.e., the virial relation is one), the distribution at qES can be described well by the Lynden-Bell statistics which assumes ergodicity in phase space based on the Vlasov equation Lecar71 ; Yamaguchi08 ; Joyce11 (Note that the Lynden-Bell statistics is not a unique theory for long-range interacting systems initially in the virial equilibrium. It has been reported that non-ergodic theories can describe the systems more accurately than the Lynden-Bell statistics deBuyl11 ; Ribeiro-Teixeira14 ; Benetti14 ). In contrast, if the initial distributions are far from virial equilibrium, the distribution at qES cannot agree with the statistics Lecar71 ; Yamashiro92 ; Yamaguchi08 ; Joyce11 ; Teles11 . The characteristic of the distributions from initial conditions near virial equilibrium is a core-halo structure. The halo is composed of particles acquired high energy by a parametric resonance of the oscillating potential Teles11 . Also, the distribution of particles losing energy approaches that of a degenerate Fermi gas, resulting in forming the core Teles11 .

Actually, in the previous researches about 1DSGS, the null virial condition was not focused on and the water-bag distributions were mainly adopted as the initial condition. In this Letter, therefore, we shall investigate qES of 1DSGS corresponding to the various initial distributions satisfying the null virial condition. Although 1DSGS is simpler than 2DSGS and 3DSGS, it possesses the essential properties of SGSs. Thus, if we could unveil the dynamics to qES of 1DSGS, by applying it to the other dimensional SGSs, we might clarify the origin of the universality of number density profiles of globular clusters or molecular clouds.

II One-dimensional self-gravitating system

The one-dimensional gravitational potential per unit mass $\phi(x)$ generated by mass source $\rho(x)$ satisfies the following Poisson equation in the one-dimensional space.

[TABLE]

When the mass source is a mass point with mass $m$ at the origin, $\rho(x)=m\delta(x)$. Thus, a single-particle potential $\phi_{\rm sp}$ generated by the mass point satisfies the following Poisson equation,

[TABLE]

and the potential can be solved as

[TABLE]

where we have used $\phi_{\rm sp}(x)=-\phi_{\rm sp}(-x)$ and $\phi_{\rm sp}(0)=0$. Therefore, an interaction force between two mass points is constant, and then, the Hamiltonian of 1DSGS of $N$ particles with equal mass can be expressed as

[TABLE]

where $x_{i}$ and $u_{i}$ mean a position and a velocity of the $i$th particle, respectively.

III $N$-body simulations

III.1 Initial conditions

As the initial distribution for $N$ body simulations, we shall adopt the following polytrope solution with a polytrope index $n$:

[TABLE]

where $r\equiv|x|$, $v\equiv|u|$, $\Theta$ is the Heaviside step function and $\mathcal{R}_{n}$ means a radius of the system. The potential energy per unit mass $\phi_{n}(r)$ satisfies the Poisson equation,

[TABLE]

where $\mathcal{N}_{n}(r)=\int{\rm d}vf_{n}(r,v)\propto\left\{\phi_{n}(R)-\phi_{n}(r)\right\}^{n}$. Between the local pressure $p_{n}(r)=\int{\rm d}vv^{2}f_{n}(r,v)$ and the number density $\mathcal{N}_{n}(r)$, the polytrope relation, $p_{n}(r)\propto\mathcal{N}_{n}(r)^{1+1/n}$, is established like the three-dimensional polytrope solution Binney08 . Moreover, the total kinetic energy $K_{n}=\int{{\rm d}}r\int{{\rm d}v}\frac{1}{2}mv^{2}f_{n}(r,v)$ and the total potential energy $\Omega_{n}=\frac{1}{2}\int{{\rm d}}r\phi_{n}(r)m\mathcal{N}_{n}(r)$ satisfy the virial relation for 1DSGS, $2K_{n}/\Omega_{n}=1$.

The polytrope solution with $n=\infty$ corresponds to the thermal equilibrium state distribution. Moreover, the number density of the the polytrope solution with $n=0$ in real space is constant, i.e., the water-bag distribution. Hence, by changing the value of $n$, we can make various distributions which are physically characterized by the polytrope relation.

DF with the virial ration $V$ can be derived from $f_{n}(r,v)$ as $f_{n,V}(r,v)\equiv f_{n}(r,v/\sqrt{V})$ where the velocity interval of integration in phase space is changed to be $0\leq v\leq\sqrt{2V\left\{\phi_{n}(R)-\phi_{n}(r)\right\}}$. Note that the number density of $f_{n,V}(r,v)$ in real space, $\mathcal{N}_{n,V}(r)$, is equivalent to that of $f_{n}(r,v)$: $\mathcal{N}_{n,V}(r)=\mathcal{N}_{n}(r)\propto\left\{\phi_{n}(R)-\phi_{n}(r)\right\}^{n}$.

When we run the $N$-body simulations, we adopt a unit system where $4\pi G=M=1$ in which $M(=mN)$ is the total mass and the radius of the case with $n=0$ (water-bag distribution) is unity.

III.2 Number densities at qES

When the null virial condition is satisfied at initial time (i.e., the initial distribution in real space is $\mathcal{N}_{n}(r)$ and all the velocities are zero), the numerical simulations about 1DSGS make it clear that there are qESs and DF in phase space is symmetric about the origin. Moreover, it is found that the number density in qES has the universality regardless of the initial polytrope index and the density around the center of the system can be fitted well by Eq. (1).

We can make a decision whether the system is at qES by paying attention to the change of physical quantities charactering the system in time, e.g. the total kinetic energy and the potential energy. As one of the examples of such physical quantities, we show the time evolution of the deviation from the equipartition law, Eq. (2), for $n=0$, 1, 5, 10 and 50 (100) in Fig. 1. The curve for $n=50$ overlay that for $n=100$ in the figure. Although the behavior of the saturation of the deviation for $n=0$ (corresponding to the water-bag distribution) is exceptional, all the deviations does not change after a time and we regard this circumstance as qES. Note that the behavior of $\Delta$ is completely different from that in Refs. Tsuchiya94 ; Tsuchiya96 : If the initial conditions is the water-bag form in virial equilibrium, $\Delta$ temporarily decreases and then increases.

Figure 2(a) shows time averaged number densities at qES of $N=10^{4}$ particles simulations where the null virial condition is satisfied, for several initial polytrope indices ($n=0$, 1, 5, 10, 50 and 100). So as to be seen easily, each density is shifted by two digits to the below. Thus, we can find that the number densities around the center of the system universally have the same form, which can be approximated by Eq. (1) shown by black curves in the same figure, independent of the initial polytrope index. The optimal values of $\kappa$ obtained by fitting the number densities with Eq. (1) are plotted in Fig. 2(b) as a function of the initial polytrope index $n$. Since the values are almost 0.5, the expectation stated in Introduction has been verified numerically. Contrary to the universality, the radius of the system increases as $n$ gets larger, which can be seen in Figs. 2(a) and (c).

IV envelope equation

In the previous section, it is found that the number density around the center of the system has the universality depicted by Eq. (1) with $\kappa\sim 0.5$ when the initial distribution satisfies the null virial condition. In order to investigate why this condition yields the universality, we utilize the envelope equation (EEq) of 1DSGS. If DF, $f(x,u,t)$, of 1DSGS satisfies the Vlasov equation (i.e., in dynamical equilibrium) and it is origin-symmetric and isotropic in velocity, EEq can be represented as follows111Note that, although EEq has already been reported for the initial water-bag distribution in Refs. Teles11 ; Levin14 , the derivation in the Letter requires only the symmetry about the origin in phase space. (see Appendix A):

[TABLE]

where $r_{e}(t)\left(\equiv\sqrt{M\left\langle x^{2}\right\rangle_{t}}\right)$ is the envelope of the system, $\Omega$ is the total gravitational potential and $\varepsilon(t)^{2}\equiv 4(\left\langle x^{2}\right\rangle_{t}\left\langle u^{2}\right\rangle_{t}-{\left\langle xu\right\rangle_{t}}^{2})$. The average $\left\langle\bullet\right\rangle_{t}$ means the integration over phase space with DF: $\left\langle\bullet\right\rangle_{t}\equiv\int{\rm d}x\int{\rm d}u\bullet f$. $\varepsilon(t)^{2}$ cannot be less than 0 because of the Cauchy-Schwarz inequality.

From Eq. (9), it is found that the dynamics of the envelope is determined by the opposite effects, collapse and expansion. If and only if the initial condition of a system satisfies null virial condition (i.e., $\varepsilon(0)=0$), the envelope begins to collapse. When the system is origin-symmetric in phase space, $\varepsilon(t)^{2}$ does not change in time: By using the Vlasov equation, the time derivative of $\varepsilon(t)^{2}$ is calculated as $\partial_{t}{\varepsilon(t)^{2}}=8\left\langle x^{2}\right\rangle_{t}\left\langle u\partial_{x}\phi\right\rangle_{t}-8\left\langle xu\right\rangle_{t}\left\langle x\partial_{x}\phi\right\rangle_{t}$. Because the system has the symmetry about the origin in phase space, $\int{\rm d}uuf=0$ holds resulting in $\left\langle u{\partial_{x}\phi}\right\rangle_{t}=\left\langle xu\right\rangle_{t}=0$. That is, if and only if the null virial condition is satisfied at initial time and the system can be described by the Vlasov equation, the system keeps on collapsing toward the origin. In other words, all the particles are attracted toward the origin regardless of the initial distribution.

From Eq. (9) with $\varepsilon^{2}=0$, $r_{e}$ seems to become 0 within a finite time. However, the actual envelope (calculated by the product of $\sqrt{M}$ and the root mean square position which derived from $N$-body simulations) does not reach 0. The time evolution of the envelope derived from $N=10^{4}$ particles simulation with $n=1$ is shown in Fig. 3 by (blue) open circles. The solution of Eq. (9) with setting $\varepsilon^{2}=0$ and replacing $\Omega(t)$ by $\Omega(0)$ is also shown in the same figure by the (red) curve. Although the solution can approximate the behavior of the actual envelope derived from the simulation around $t=0$, it cannot describe the minimum of the actual envelope.

The deviation in behavior of the actual envelope from the solution of Eq. (9) with $\varepsilon^{2}=0$ is due to the increment of $\varepsilon^{2}$. The reason why $\varepsilon^{2}$ becomes non-zero is that it becomes impossible to describe the system by the Vlasov equation, as the system is collapsing toward the origin. This fact can be understood from the sum of the difference in the one particle energy per unit mass from the initial time which can be calculated as (see Appendix B)

[TABLE]

where $\Delta\mathcal{E}_{i}(t)\equiv\mathcal{E}_{i}(t)-\mathcal{E}_{i}(0)$ in which $\mathcal{E}_{i}(t)$ is one particle energy per unit mass of the $i$th particle at time $t$. Since the total gravitational potential decreases during the collapse, the sum of the difference in the one particle energy per unit mass becomes negative. This means that exchanges of energy occur among the particles, and then, the Vlasov equation cannot describe this situation of the system since the equation is valid for collisionless SGSs.

V Core-Halo structure

To investigate the evolution of the system after the collapse caused by the initial null virial condition, let us focus on energy of each particle.

In Fig. 4, we show the time evolution of the kinetic, potential and total energy of the outermost particles at initial time. The total energies of particles initially at the outside of the system (outer particles) gradually increase like this example.

One can simply explain the mechanism of the increment of energy of the outer particles. Since the system is origin-symmetric, the gravity acting on a particle at $r$ is proportional to the number of particles within a smaller distance than $r$. While all the particles move toward the origin owing to the null virial condition, they cannot pass each other. Note that a difference occurs in the time of arrival at the origin: The particles initially at the inside of the system (inner particles) can reach the origin faster than the outer particles (As discussed in the previous section, the actual envelope does not become 0, which means that all the particles do not reach the origin at the same time). Therefore, each energy of the outer particles does not change until when the inner particles go through the origin and pass the outer particles. Then, the outer particles gain the potential energy because of the spread of the inner particles. After that, the outer particles also spread out from the origin. Because the inner particles return to the origin again, the second passing between the inner and outer particles occurs. After these movements are repeated for a several times, the system reaches qES.

While the outer particles acquire energy until qES, there are particles losing energy. This fact can be clear from the following equation,

[TABLE]

where $\Delta\mathcal{E}_{i}^{\rm qe}$ represents the difference of energy of the $i$th particle between the initial time and qES (see Appendix B). Since the total gravitational potential energy for 1DSGS is positive, the particles reducing energy must exist.

The movements of the particles gaining and losing energy at qES are as follows. A region where the lower energy particles can move in real space is restricted to the neighborhood of the origin and we refer to the region as a core. Because the particles at the core are influenced by many two-body interactions each other, the trajectories of them change randomly. Contrary to this, the high energy particles can move in a much wider region in real space which we refer to as a halo. Owing to the high energy, the halo particles are not influenced by two-body interaction when passing the core, which means that the movements of them are governed by the mean potential of the system and the trajectories are smooth. These are a scenario for constructing the core-halo structure of the system satisfying the initial null virial condition.

VI Fokker-Planck model for the core distribution

Here, we shall derive the core distribution by using a special Fokker-Planck equation. The Fokker-Planck model approach is appropriate rather than the kinetic equations, i.e., the Boltzmann equation and the Landau equation, because the collision terms of these equations are zero in 1DSGS with equal mass Eldridge63 ; Kadomtsev70 ; Chavanis06 ; Sano11 .

Before constructing the Fokker-Planck equation, we shall model a force influencing an element of the core particles. That is, we shall begin by constructing a Langevin equation. Since the system is origin-symmetric, we use $r(t)\equiv|x(t)|$ where $x$ is the position of the element.

Let us introduce the frictional force $-m\gamma\dot{r}(t)$ and the random noise with constant intensity $\sqrt{2D}\xi(t)$, which are essential for a many-body system to reach the thermal equilibrium state, where $\xi(t)$ means a Gaussian-white noise. Also, the element must be influenced by a mean gravitational force $-F(r)$, which is derived by differentiating $m\phi(r)$ as $-F(r)=-m\partial_{r}\phi(r)$ where $\phi(r)$ is the mean gravitational potential at qES. However, this is just a mean gravity. It is natural to consider that the element actually is influenced by a fluctuating gravity around the mean value: The number of particles within a smaller distance than $r$ producing $-F(r)$ through the Poisson equation is the mean value, and the actual number must fluctuate around the value. This means that another noise which prevents the system from reaching the thermal equilibrium state simply is added to the ordinary Langevin equation, and so this system goes to another stable state, that is, qES. Therefore, we can regard the noise as distinctive to 1DSGS.

If assuming the intensity of the noise to be a constant denoted by $\sqrt{2\epsilon}$, we can obtain the following Langevin equation:

[TABLE]

where $\eta$ also means a Gaussian-white noise and has no correlation with $\xi$. In the over-damped limit, the equation (12) becomes

[TABLE]

The second term on the right hand side of the above equation is a correction term in order to regard products as the Storatonovich product Sekimoto99 .

The corresponding Fokker-Planck equation is given by

[TABLE]

When the system reaches qES, $\partial_{t}P^{\rm qe}(r)=0$. Here, by integrating the Fokker-Planck equation by $r$, we can obtain

[TABLE]

Let us determine the constant of the right-hand side of the above equation by using the boundary condition at $r=0$. Because of the symmetry, the mean field force can be represented as $F(r)\propto\int_{0}^{r}{\rm d}r^{\prime}NP^{\rm qe}(r^{\prime})$. Thus, $F(0)=0$ where we have used the fact that $P^{\rm qe}(0)$ is bounded. Moreover, $P^{\rm qe}{}^{\prime}(0)=0$, because of the symmetry of the system222We shall introduce $\tilde{P}^{\rm qe}(x){\rm d}x$ which represents a probability that we can find an element in a range from $x$ to $x+{\rm d}x$ at qES. Because of this definition and the symmetry of the system, $\tilde{P}^{\rm qe}(-x)=\tilde{P}^{\rm qe}(x)=P^{\rm qe}(r)/2$. By differentiating the equation with $x(>0)=r$, we can obtain $-\tilde{P}^{\rm qe}{}^{\prime}(-x)=\tilde{P}^{\rm qe}{}^{\prime}(x)={P^{\rm qe}}^{\prime}(r)/2$. Therefore, $\tilde{P}^{\rm qe}{}^{\prime}(0)={P^{\rm qe}}^{\prime}(0)=0$.. Therefore, the constant is determined to be equal to 0. In addition, if we use the number density at qES, $\mathcal{N}^{\rm qe}(=NP^{\rm qe})$, the equation (15) becomes

[TABLE]

The mean gravitational force $F(r)$ satisfies the following equation

[TABLE]

which can be derived by substituting the gravitational potential per unit mass $\phi(r)(=\frac{1}{m}\int{\rm d}rF(r))$ into the Poisson equation \mbox{\bigtriangleup}\phi=4\pi G\rho^{\rm qe}=4\pi Gm\mathcal{N}^{\rm qe}(r).

Now, we nondimensionalize these equations by using the following units of length and force:

[TABLE]

and

[TABLE]

where $\mathcal{T}=D/\{\gamma+56\pi\epsilon Gm\mathcal{N}^{\rm qe}_{0}\}$ and $\mathcal{N}^{\rm qe}_{0}\equiv\mathcal{N}^{\rm qe}(0)$. Then, equations (16) and (17) are altered to

[TABLE]

and

[TABLE]

where $q\equiv 4\pi\epsilon Gm\mathcal{N}^{\rm qe}_{0}/\gamma$ and variables with overbar denote dimensionless. We should solve these equations with boundary condition $\bar{\mathcal{N}}^{\rm qe}(0)=1$.

The solution of the equations can be represented as

[TABLE]

and

[TABLE]

where $f_{q}^{-1}$ is the inverse function of $f_{q}$ defined as

[TABLE]

In Fig. 5, we show the non-dimensional number density $\bar{\mathcal{N}}^{\rm qe}$, obtained by calculating the inverse function of $f_{q}$ numerically, as a function of $\bar{r}$. Each curve corresponds to the number density with $q=0$, 1, 3, 4 and 6 respectively from left to right. The dashed curve, representing $1/(1+\bar{r}^{2})^{1/2}$, means the typical number density around the core at qES found numerically. Thus, we can describe the numerical results by our model with appropriate $q$.

We can examine the dependence of $\kappa$ in Eq. (1) on $q$ from the analytical formula of the number density at qES. By approximating the number density around the origin as $\bar{\mathcal{N}}^{\rm qe}(\bar{r})\simeq 1/(1+\bar{r}^{2})^{\kappa}$, we can obtain

[TABLE]

Since $q\geq 0$, the range of $\kappa$ is $3/7<\kappa\leq 3$, which includes the results of numerical simulations. However, the numerical results show that $\kappa\simeq 0.5$ which means that the value of $q$ is limited. Because $q$ is proportional to a ratio of the intensity of mean gravity fluctuation $\epsilon$ to the friction coefficient $\gamma$, this limited value of $q$ can be regarded as a kind of fluctuation-dissipation relation Kubo91 .

Finally, we consider that the number density with $q=0$. Note that $q$ is proportional to $\epsilon$ which prevent the system from reaching the thermal equilibrium state. Thus, the number density with $q=0$ must agree with the thermal equilibrium state of 1DSGS.

By substituting $q=0$ and Eq. (23) on Eq. (24), we can obtain

[TABLE]

Therefore, $\bar{F}$ can be solved as

[TABLE]

and then,

[TABLE]

which corresponds to the number density of 1DSGS in the thermal equilibrium state Rybicki71 .

VII estimation of the radius of the system

In this section, we shall estimate the radius of the system at qES. As discussed in Sec. V, the particles gaining energy at the initial time construct the halo. The size of the system can be understood as the range of the movement of the highest-energy particles among them. We shall approximate the mechanism about gaining the highest energy as follows: When the outermost particles at initial time reach the origin, the other particles go back to the initial positions instantaneously, and then, the outermost particles gain the highest energy. Therefore, the increment of energy of the outermost particles from the initial time to qES can be estimated as

[TABLE]

Since the outermost particles initially have the energy $\mathcal{E}_{\rm out}(0)=2\pi Gm\sum_{i}|x_{i}(0)-x_{\rm out}(0)|$ where $x_{\rm out}$ represents the position of the outermost particles, the energy of them at qES can be approximated as

[TABLE]

The system is assumed to be origin-symmetric. Then, the energy of the highest-energy particles at rest for a moment becomes $2\pi GM\mathcal{R}_{\rm est}^{\rm qe}$ where $\mathcal{R}_{\rm est}^{\rm qe}$ is the estimated radius of the system. Hence, the radius of the system can be estimated as

[TABLE]

In Fig. 6, the estimated radiuses of $10^{4}$ particle systems at qES from Eq. (31) are plotted by triangles with the actual radiuses of the same systems shown by circles with varying the initial polytrope index $n$. We can see that the estimated radiuses almost corresponds to the actual ones although we roughly approximate the mechanism that the outermost particles gain energy.

VIII description of the whole of density profile at qES by our model

Now, we shall describe the whole of density profile of 1DSGS at qES by the Fokker-Planck model for the core and the estimated radius as follows. First, the radius of the system $\mathcal{R}_{\rm est}^{\rm qe}$ is estimated from Eq. (31) with the initial distribution $x_{i}(0)$. Next, we use constrained conditions about the total number and the total energy. The total number must satisfy

[TABLE]

which can be computed as

[TABLE]

where $a$ represents the unit of length. With this relation and Eq. (24), the following relation can be obtained:

[TABLE]

The total energy $E$ can be expressed only by the total gravitational potential owing to the virial relation as $E=3\Omega^{\rm qe}/2$. Since $\Omega^{\rm qe}$ can be computed by integrating $\bar{\mathcal{N}}^{\rm qe}(\bar{r})$, the constrained condition for $E$ is (see Appendix C)

[TABLE]

Finally, the parameter $q$ the constrained conditions contain must be decided. Because there are no conditions for deriving the parameter, we shall obtain the value by fitting the power law behavior of the numerical number density at qES with $r^{-2\kappa}$ and utilizing Eq. (25).

From the initial distribution of particles, $N$, $E$ and $q$ calculated in this manner, we compute $\mathcal{R}_{\rm est}^{\rm qe}$, $\mathcal{N}^{\rm qe}_{0}$ and $a$ by Eqs. (34) and (35). Figure 7(a) shows the number density of $10^{4}$ particle system with the initial polytope index $n=1$ at qES by open-circles. The curve means the number density derived by the above method. We can see that our model can describe $\mathcal{R}_{\rm est}^{\rm qe}$, $\mathcal{N}^{\rm qe}_{0}$ and $a$ well. Figure 7(b) shows the degree of agreement between our model ($\mathcal{N}^{\rm qe}_{0}\bar{\mathcal{N}}^{\rm qe}$) and the numerical result ($\mathcal{N}^{\rm num}$) defined by

[TABLE]

with varying the initial polytrope index $n$. The same degree of agreement is found for the other initial polytrope indices.

IX Concluding remarks

In the Letter, we have made it clear numerically that 1DSGS has the qES when the initial conditions satisfy the null virial condition and the density profile around the origin is found to have the universality independent of initial profile in real space. The universal density profile around the center is approximated by Eq. (1) with $\kappa\simeq 0.5$.

The reason why the null virial condition is sufficient for the universality is unveiled by the envelope equation. In general, this equation contains the opposite effects, collapse and expansion. If and only if the initial condition satisfies the null virial condition and the system can be described by the Vlasov equation, the envelope keeps on collapsing toward the origin. In other words, for the other virial conditions, the envelope is attracted to the other points depending on the initial condition. After the collapse, the particles exchange energy each other resulting in the core-halo structure: The core and the halo are composed by the particles gaining and losing energy, respectively. When the initial distribution is the water-bag form and near virial equilibrium, the exchange of energy also causes the core-halo structure Teles11 . Note that, however, a mechanism that the halo particles acquire high energy is a parametric resonance of the oscillating potential and is quite different from what is discussed in the Letter.

Since the core particles are influenced by the two-body interactions many times over, the trajectories of them change randomly. Thus, we have modeled the movement by the special Fokker-Planck equation, Eq. (14), reflecting a fluctuation of the mean gravitational force. The number density profile derived as a steady solution of Eq. (14) has no contradiction with the numerical results.

Contrary to this, since the particles of the halo have high energy, the movements are governed by the mean gravitational potential and the trajectories are smooth. We have approximated the mechanism of gaining the highest energy and estimated the radius of the system by the range of the movement of the particles. The estimated radiuses almost correspond to the real ones as shown in Fig. 6.

We also have compared the number density profiles of our model with the estimated radius, with the numerical simulation results. The number density of our model can be obtained by the initial positions of the particles, the total energy, the total number of particles and the index of the power law behavior of the number density. From Figs. 7(a) and (b), we can see that our model can describe the numerical results well for the case of $n=1$ and the same degree of agreement is found for the other initial polytrope indices.

Although a factor causing $\kappa\sim 0.5$ remains unclear, the dynamics to the core-halo stricture is unveiled and the movements at the core and the halo can be explained in the Letter. Thus, it is expected that the scenario to the universal density profile of 1DSGS will help to clarify the origin of the universality of globular clusters or molecular clouds. Moreover, from the fact obtained numerically in the Letter and the universal density profile of 2DSGS and 3DSGS Tashiro10 ; Tashiro16 , we can make a conjecture that the density profile of $n$ dimensional SGS at qES from the null virial condition is described by Eq. (1) with $\kappa\simeq n/2$. We can consider that the conjecture gives a new direction to the research about SGS and must be examined.

X Acknowledgments

The author would like to thank members of astrophysics laboratory at Ochanomizu University for extensive discussions. He also appreciates the advice about parallel computing of $N$-body simulations from T. Takayuki.

Appendix A derivation of the envelope equation

Since $f(x,u,t)$ satisfies the Vlasov equation, $\partial_{t}f=-u\partial_{x}f+\partial_{x}\phi\partial_{u}f$, the second derivative of $r_{e}(t)$ with respect to time can be calculated as

[TABLE]

where the average $\left\langle\bullet\right\rangle_{t}$ means the integration over phase space with DF: $\left\langle\bullet\right\rangle_{t}\equiv\int{\rm d}x\int{\rm d}u\bullet f$. By multiplying $M\left\langle x^{2}\right\rangle_{t}/r_{e}(t)^{2}\ (=1)$ by the first term on R.H.S., the equation becomes

[TABLE]

where $\varepsilon(t)^{2}\equiv 4(\left\langle x^{2}\right\rangle_{t}\left\langle u^{2}\right\rangle_{t}-{\left\langle xu\right\rangle_{t}}^{2})$.

Let us assume that the system is origin-symmetric in phase space. Then, with $r\equiv|x|$, $v\equiv|u|$, $\left\langle x\partial_{x}\phi\right\rangle_{t}$ is calculated as

[TABLE]

where $\mathcal{R}$ and $\mathcal{V}$ are respectively the radius and the maximum velocity of the system, and we have used the Poisson equation, $\partial_{r}^{2}\phi=4\pi G\rho$ in which $\rho=M\int 2{\rm d}vf$.

The total gravitational potential, $\Omega\ (=\int 2{\rm d}r\rho 2\pi Gr\int 2{\rm d}r^{\prime}\rho)$, also can be calculated as

[TABLE]

where we have used $\left.\partial_{r}\phi\right|_{r=0}=0$, and then,

[TABLE]

Therefore, we can obtain EEq,

[TABLE]

Appendix B sum of the difference in the one particle energy per unit mass from the initial time

The total energy of 1DSGS, $E$, can be described as

[TABLE]

Let us introduce the gravitational potential per unit mass at $x$, $\phi(x)\equiv\sum_{j}2\pi Gm|x-x_{j}|$, and then, the above equation becomes

[TABLE]

where $\mathcal{E}_{i}\equiv\frac{1}{2}{u_{i}}^{2}+\phi(x_{i})$. Thus, the sum of the one particle energy per unit mass can be depicted as

[TABLE]

Then, the sum of the difference in the one particle energy per unit mass from the initial stag can be obtained as

[TABLE]

where $\Delta\mathcal{E}_{i}(t)\equiv\mathcal{E}_{i}(t)-\mathcal{E}_{i}(0)$.

By calculating the difference of this equation between the initial time and qES, we can obtain

[TABLE]

where $\Delta\mathcal{E}_{i}^{\rm qe}\equiv\mathcal{E}_{i}^{\rm qe}-\mathcal{E}_{i}(0)$ and the superscript qe indicates the quantity at qES.

If the system satisfies the null virial condition at the initial time, the energy conservation law can be written as

[TABLE]

At qES, the virial theorem must hold: $2K^{\rm qe}=\Omega^{\rm qe}$. Therefore,

[TABLE]

and the equation (47) becomes

[TABLE]

Appendix C Constrained condition for the total energy

The total gravitational potential at qES, $\Omega^{\rm qe},$ can be computed by integrating $\bar{\mathcal{N}}^{\rm qe}$ as

[TABLE]

where $\bar{\mathcal{R}}^{\rm qe}_{\rm est}\equiv\mathcal{R}^{\rm qe}_{\rm est}/a$ and we have used Eq. (21) and $\bar{F}(0)=0$.

We can perform the integral of the second term on the last R.H.S of the above equation by substituting $\bar{F}(\bar{r})=f_{q}^{-1}(\bar{r})=X$:

[TABLE]

where we have used Eqs. (33) and (34).

Thus, by using Eq. (33) and $E=3\Omega^{\rm qe}/2$, we can obtain

[TABLE]

and Eq. (35).