text:
# Origin of scaling structure and non-gaussian velocity distribution in self-gravitating ring model

## I Introduction

Many astrophysical objects in our universe consist of mutually interacting elements through gravity. If they are almost isolated systems, they are called self-gravitating system (SGS). For example, galaxies, clusters of galaxies, globular clusters, and molecular clouds are thought to be typical SGSs. Their statistical properties are often characterized by non-Gaussian velocity (or pair-wise velocity) distributions, fractal structures, and the scaling relation between the mass density and the system size. Most of these objects are thought to be gravitationally virialized. Therefore pure gravitational force seems to play an essential role in characterizing the above statistical properties of SGS independently of initial conditions.

There have been some theoretical approaches to explain the fractal structures in SGS from the viewpoint of criticality and phase transition in gravo-thermo dynamics. Strictly speaking, the ordinary SGSs in three-dimension (3-D) cannot attain genuine stable equilibrium, because the gravitational force does not vanish at long distances (IR-divergence) and diverges at short distances (UV-divergence). These properties of gravity cause gravothermal catastrophe in a self-gravitating gas system enclosed by a solid adiabatic wall. In fact, the isothermal sphere is not always stable, since the entropy does not necessarily take the local maximum for this configuration. Therefore the introduction of a small-scale cutoff as well as a large-scale cutoff is inevitable for such unstable systems to attain the final equilibrium. Though the introduction of the cutoff prevents the gravothermal catastrophe, 3-D gravitational system has a phase with negative specific heat when treated in the microcanonical ensemble. Even in the canonical ensemble, this system is highly unstable. When the temperature decreases, the system experiences the violent first-order phase transition from the gas phase into the cluster phase. According to these arguments, no stable equilibrium states are theoretically expected in the 3-D gravitational system with or without cut-off in any ensemble.

On the other hand, in the real world, we do not need the stable equilibrium states for describing SGS; meta-stable equilibrium states do actually appear in the dynamical description such as collisionless Boltzmann equation. The finite lifetime of such structures is sufficient to explain the present structures of SGS, even if they are expected to evolve further into different quasi-stable states through two-body relaxation.

The relaxation process of SGS has been mainly discussed in the one-dimensional gravitational sheet model (1DS) . In this model, many parallel sheets interact with each other through constant force which never decays at distant places. Though the interaction is long-ranged, no phase transition occurs in 1DS. Thermodynamics of 1DS is exactly solved and actually, in numerical calculations, the system reaches thermal equilibrium long after it attains virial equilibrium. In this model, the virial condition gives the relation $`2K=V`$ between the time averaged kinetic energy $`K`$ and the potential energy $`V`$. Therefore, contrary to the 3-D SGS, specific heat of 1DS is always positive. Thus the relaxation process in 1DS would be quite different from that in 3-D SGS. Another well-known one-dimensional model which has long-range force is the Hamilton Mean Field (HMF) model, in which phase transition does occur. There have been a lot of studies on the relaxation process of HMF. Actually, Lévy type jumping motion of constituent particles and super-diffusion process have been revealed in HMF. In this HMF, though the specific heat becomes negative in the quasi-stable state , it remains positive in the thermal equilibrium state and the phase transition turns out to be the second-order. Despite some common properties of the models, the interaction form in HMF is also quite different from that of real 3-D gravitational systems.

Several other simple models have been proposed in order to characterize the 3-D SGS much faithfully. For example, the cell model and its extended version are the simplest models which show phase transition. In these models, the pair interaction potential changes its value only at the cell boundary. At low temperatures, most particles are trapped within several clusters, and at higher temperatures, these clusters melt and the particles can move freely. However, it seems difficult to examine the relaxation process of this system numerically, since the interaction is not efficient to cause sufficient relaxation; particles have no interaction within each cell. There are also several numerical analysis with the different type interactions including the lower and higher cutoffs. However, the value of the lower cutoff used in realistic numerical calculations of these models seems to be too large to extract the intrinsic properties of realistic 3-D SGS.

Thus it seems indispensable for us to have a model which faithfully reflects the characteristic properties of 3-D SGSs; the model has negative specific heat and shows the phase transition similar to 3-D SGSs, and is capable for us to analyze quasi-equilibrium states which would be realized before the system reaches the complete equilibrium.

In this paper, therefore, we first propose a new model (Self Gravitating Ring: SGR) which a) has negative specific heat, b) shows phase transition representing 3-D SGS, and c) is numerically tractable. This model consists of $`N`$-particles, which can freely move on a ring with a fixed radius, mutually interacting through 3-D gravity. An excellent point of this model is that the force is genuine 3-D gravity while the calculation is essentially one-dimensional. Moreover, the Hamiltonian permits the accurate symplectic integration method by which we can analyze the nature of quasi-equilibrium states and phase transitions within very long time scales in this model.

In section II, we introduce the SGR model, and in section III identify three quasi-equilibrium states including the state with negative specific heat. In section IV, we analyze the particle motions from a statistical point of view, and then, we study the relaxation process toward the thermodynamical equilibrium in section V, and the scaling structure in section VI. Finally, we discuss our results in section VII.

## II Self-gravitating ring model

In this section, we introduce the self-gravitating ring (SGR) model, in which particles interact with each other through genuine 3-D gravity while the particle motion is constrained on a 1-D ring. By utilizing this model, we study the gravitational phase transition, which was not feasible in the 1DS model. The lack of phase transition in 1DS is related with the fact that the potential of 1DS increases linearly without bound, and therefore there is no characteristic energy scale required for a particle to escape from the cluster. On the other hand, we have the characteristic energy scale to bind particles in our SGR model.

We consider a system of self-gravitating $`N`$ particles with mass $`m`$, whose motions are smoothly constrained on a circular ring with a fixed radius $`R`$ without friction (Fig.1). Each pair of particles interact with each other through 3-D gravitational force. The distance between a pair is measured by the length of the straight line combining the pair and not by the minor arc of the ring.

The Hamiltonian of this system becomes

$`H_\mathrm{p}={\displaystyle \frac{1}{2mR^2}}{\displaystyle \underset{i=1}{\overset{N}{}}}P_i^2{\displaystyle \underset{i<j}{\overset{N}{}}}{\displaystyle \frac{Gm^2}{\sqrt{2}R\sqrt{1\mathrm{cos}(\theta _i\theta _j)+ϵ}}}.`$ (1)
The position of $`i`$-th particle is fully described by the angular variable $`\theta _i`$ as $`\text{ }𝒓_i=(R\mathrm{sin}\theta _i,R\mathrm{cos}\theta _i)`$. The momentum conjugate to $`\theta _i`$ is given by $`P_i=mR^2d\theta _i/dt`$. The UV-cutoff parameter $`ϵ`$ truncates the diverging gravitational force at around the distance $`\xi \sqrt{2ϵ}R`$.

We first introduce three dynamical time-scales which are apparent in the above Hamiltonian. They are parameterized by the ring radius $`R`$, the cutoff scale $`\xi `$, and the total energy $`E`$ of the system.

When the system is almost uniformly filled by moving particles, we have the longest dynamical time $`t_R`$ defined by

$$t_R\sqrt{\frac{R^3}{GNm}}.$$
(2)
During this time $`t_R`$, a typical particle goes around the ring once. Therefore, the time-scale for the whole system to attain thermodynamical equilibrium, if any, is at least larger than this time scale.

When all the particles collapse completely into a core, we find the shortest dynamical time-scale $`t_\xi `$ defined as

$$t_\xi \sqrt{\frac{\xi ^3}{GNm}}.$$
(3)
During this time $`t_\xi `$, a typical particle bound to the core oscillates once. *** In the next section, we classify all the particles into three species, core , halo , and gas . The typical time scales $`t_R`$ and $`t_\xi `$ respectively characterize the gas and the core species.

There is an another time-scale in between the above two extreme time scales. Suppose the system is stably confined in a region with a scale $`r`$ which satisfies

$`\xi <<r<<R.`$ (4)
For this scale $`r`$, the ring can be approximated as an infinite straight line and the cutoff $`\xi `$ can be neglected. Then the leading term of the denominator of the potential term in the Hamiltonian (1) becomes $`|\theta _i\theta _j|`$ and

$`H_p{\displaystyle \frac{1}{2mR^2}}{\displaystyle \underset{i=1}{\overset{N}{}}}P_i^2{\displaystyle \underset{i<j}{\overset{N}{}}}{\displaystyle \frac{Gm^2}{R|\theta _i\theta _j|}}.`$ (5)
In this case, the potential term $`V_\mathrm{p}`$ satisfies Euler’s theorem for homogeneous function,

$$\underset{i=1}{\overset{N}{}}\theta _i\frac{V_\mathrm{p}}{\theta _i}=V_\mathrm{p},$$
(6)
and the ordinary virial condition holds

$`2K_\mathrm{p}=V_\mathrm{p},`$ (7)
where angle brackets represent the long-time average.  Outside of the above region in SGR model, this form of virial relation would be modified even in the limit of $`\xi 0`$,as is shown in Appendix A. According to this virial relation, the typical size of a system $`r`$ is related with the total energy $`E`$ as

$$E=\frac{V_\mathrm{p}}{2}\frac{Gm^2N^2}{4r}.$$
(8)
Hence the above condition for $`r`$ (4) reads

$$\frac{Gm^2N^2}{4\xi }|E|\frac{Gm^2N^2}{4R}.$$
(9)
Moreover from the virial relation, the velocity dispersion is given as

$$2K_\mathrm{p}=V_\mathrm{p}=mNv^2,$$
(10)
where $`vdr/dt`$ and therefore

$$\sqrt{v^2}=\sqrt{\frac{2|E|}{mN}}.$$
(11)
From this equation (11), the crossing time is defined as

$$t_E=\frac{r}{\sqrt{v^2}}=\frac{Gm^{5/2}N^{5/2}}{4\sqrt{2}|E|^{3/2}}.$$
(12)
This is the intrinsic time scale associated with genuine gravity independent of the cutoff $`\xi `$ and of the system size $`R`$.  As we will show in Sec.IV, this is also the time scale for the particles in halo species in SGR model.

The above introduced cutoff parameter $`\xi `$ connects the two limiting cases in the following sense. In the limit $`\xi 0`$, the SGR model becomes genuine 3-D gravity at small scales, while in the limit $`\xi R`$, the model almost becomes the HMF model. The latter is because, in this limit, the shape of the potential around $`\theta =k\pi (k=0,\pm 1,\mathrm{})`$ becomes almost identical as does in the HMF model §§§ The exact HMF model is reproduced in the limit $`\xi \mathrm{}`$. .

For numerical simulations, we need to make all physical variables non-dimensional; we use $`m`$, $`R`$, and $`t_R`$ for the unit of mass, distance, and time, respectively. In these units, the physical Hamiltonian (1) reads

$$H_\mathrm{p}=Gm^2\frac{N}{R}H,$$
(13)
where

$$H=\frac{1}{2}\underset{i=1}{\overset{N}{}}p_i^2\underset{i<j}{\overset{N}{}}\frac{1}{\sqrt{2}N\sqrt{1\mathrm{cos}\theta _{ij}+ϵ}}.$$
(14)
The dimensionless momentum $`p_i`$ is given by $`p_i=d\theta _i/d\tau `$ and the dimensionless time $`\tau `$ is introduced as $`\tau t/t_R`$. The above form of Hamiltonian permits us to use a powerful symplectic integrator , with which the total energy is conserved with extremely high accuracy, even beyond thousands of dynamical time. Typical magnitudes of errors for the total Hamiltonian $`H(\tau )`$ and the total momentum $`P(\tau )`$ in our simulations up to $`\tau 10^4`$, are $`(H(\tau )H(0))/H(0)𝒪(10^5)`$ and $`P(\tau )/P_{\text{r.m.s}}(\tau )𝒪(10^8)`$, respectively, where $`P_{\text{r.m.s}}(\tau )\left(_{i=1}^Np_i^2\right)^{1/2}`$.

## III Classification of phases and particles

We now study the quasi-equilibrium state that appears in a transient stationary stage in our SGR system. Though this state is not absolutely stable, it generally appears in SGSs during sufficiently long time before the system finally approaches the equilibrium state characterized by the equipartition of particle energy.

We would like to extract universal properties observed in this transient state; only the transient description is possible and necessary to explain observations.

### A Negative specific heat of SGS

First, we study the phase diagram of SGR model in the temperature-energy plane. The $`T`$-$`U`$ relation is shown in Fig.2 and we observe that the region with negative specific heat, i.e. negative slope region, apparently exists. The temperature $`T`$ and the internal energy per particle $`U`$ of the present system are defined respectively by

$`T`$ $``$ $`{\displaystyle \frac{2K(\tau )}{N}},`$ (15)
$`U`$ $``$ $`{\displaystyle \frac{H}{N}}.`$ (16)
The internal energy $`U`$ is conserved in the present microcanonical system, while the temperature $`T`$ is defined as the twice of the time-averaged (represented by angle brackets) kinetic energy per particle. The phase diagram in Fig.2 is time averaged until $`\tau =2`$ where the system archives virial equilibrium (see Fig.3). Throughout this paper, except explicitly mentioned, we fix the total number of particles $`N=100`$ for simplicity and therefore the remaining relevant parameters which characterize the quasi-equilibrium state of the system would be $`ϵ`$ and $`U`$.

From Fig.2, we find two characteristic energy scales where $`T/U=0`$; $`U_{\text{c-low}}(ϵ)`$ at a low energy side and $`U_{\text{c-high}}(ϵ)`$ at a high energy side. The energy scale $`U_{\text{c-high}}(ϵ)`$ corresponds to a mean gravitational binding energy per particle, which is estimated as $`1/(\sqrt{2}N\sqrt{1\mathrm{cos}\theta +ϵ})\times (N(N1)/2)/N𝒪(1)`$, while $`U_{\text{c-low}}(ϵ)`$ strongly depends on the cutoff $`ϵ`$. This cutoff dependence can be estimated from the condition (9) which, in our normalization, becomes

$$\frac{1}{4\sqrt{2ϵ}}U\frac{1}{4}.$$
(17)
Under this condition, the negative specific heat condition in virialized state would be justified. Actually, substitution of three different cutoffs $`ϵ=2.5\times 10^7`$, $`1.0\times 10^6`$, and $`2.5\times 10^5`$ into (17) yields the lower limits of (17) $`353`$, $`177`$, and $`35.3`$, respectively. Thus the condition (17) correctly describes the region of negative specific heat in Fig.2 with sufficient accuracy. However, the slope of each T-U line in Fig.2 is less steep than the value $`2`$ which is expected for the virial condition of 3-D gravity. This discrepancy is probably because the energy range satisfying the condition (17) is too narrow for the ideal T-U curves according to our choice of the cutoff parameters; even the smallest cutoff we took may not be enough to make the relation converge.

In the range between these two energy scales, i.e. $`U_{\text{c-low}}(ϵ)<U<U_{\text{c-high}}(ϵ)`$, the system has negative specific heat, which suggests the existence of phase transition . Actually in the system with negative specific heat, a slight energy injection from outside decreases the system’s temperature and induces further energy flow from outside. Then this catastrophic temperature reduction induces rapid cluster formations in the system. As we will see soon below, such phase transition from the gaseous state to the cluster state actually appears and characteristic structures are realized in this intermediate energy range.

### B Three phases in SGR model

As is seen in Fig.2, there are apparently three phases according to the energy per particle $`U`$; a) low energy collapsed phase ($`𝒞`$ -phase) $`U<U_{\text{c-low}}(ϵ)`$, b) intermediate energy phase ($``$ -phase) $`U_{\text{c-low}}(ϵ)<U<U_{\text{c-high}}(ϵ)`$, and c) high energy gaseous phase ($`𝒢`$ -phase) $`U_{\text{c-high}}(ϵ)<U`$. The $`𝒢`$ -phase (c) is stabilized by the infrared cutoff ($`\theta 2\pi `$, or a largest physical scale $`R`$), without which the particles would escape into spatial infinity. The $`𝒞`$ -phase (a) is stabilized by the ultraviolet cutoff $`ϵ`$, without which the particles would fall into a singularity. The specific heat for these particles (a and c) is positive, in accordance with the stability of these phases. On the other hand, the $``$ -phase (b) (the most specific to gravity) has negative specific heat and therefore is unstable. The nature of this phase is independent of any artificial cutoffs and therefore is thought to represent intrinsic properties of gravity.

The existence of three different phases, as is explained above, distinguishes the SGR model from other models with long-range force. For example, HMF model has only high and low energy phases, and shows a second-order phase transition between these phases. On the other hand in SGR model, the intermediate phase with negative specific heat exists and is strongly unstable when the system is in contact with a heat bath. As for 1DS model, only a single phase exists because the system has no characteristic energy scale; there is no phase transition in this model.

The limit $`ϵ0`$ in SGR model represents genuine gravity without collision and the range of the intermediate phase increases without bound.

### C Three species: gas , halo , and core particles

In each phase, the particles of the system prevail in various energy ranges. For example, in the intermediate phase, some particles evaporate from a cluster and move along the ring almost freely with the time scale $`t_R`$, while some of the others fall into the center of a cluster and oscillate with the time scale $`t_\xi `$. Thus, the overall phase information does not precisely specify the nature of individual particles. In order to obtain much fine information, we define three species in each phase by using the energy of the particles: The energy of the $`i`$’th particle is given by

$$E_i\frac{1}{2}p_i^2\underset{ji}{\overset{N}{}}\frac{1}{\sqrt{2}N\sqrt{1\mathrm{cos}\theta _{ij}+ϵ}}.$$
(18)
The classification is based on particle energy, that is, a) core particles for $`E_i<U_{\text{c-low}}(ϵ)`$, b) halo particles for $`U_{\text{c-low}}(ϵ)<E_i<U_{\text{c-high}}(ϵ)`$, and c) gas particles for $`U_{\text{c-high}}(ϵ)<E_i`$.  Note that even in the $`𝒞`$ -phase, there exist a few gas particles. We term the core particles and halo particles as a cluster in this paper; they form an apparent single bound state.

In the low-temperature phase of SGR, the quasi-equilibrium state at very low temperature is highly inhomogeneous and a single cluster is formed. Most particles are condensed in this cluster and the total potential is deep. In the intermediate phase, many particles spread around a cluster. As the temperature increases within this phase, halo particles gradually dominate core particles and eventually there appear gas particles, which evaporate from a cluster and go along the ring. In the high-temperature phase, all particles move almost freely without forming any cluster.

We describe the ratio of particle number in each states  We term a state to represent the situation that a particle belongs to one of the three species.  by $`\lambda _{\text{core}},`$ $`\lambda _{\text{halo}},`$ and $`\lambda _{\text{gas}}`$, with $`\lambda _{\text{core}}+\lambda _{\text{halo}}+\lambda _{\text{gas}}=1`$. Their evolution is shown in Fig.4 for the case $`ϵ=2.5\times 10^7`$. From this, we observe those ratios seem to have approached asymptotic values beyond $`\tau 1`$: $`\lambda _{\text{core}}0`$, $`\lambda _{\text{halo}}0.2`$, and $`\lambda _{\text{gas}}0.8`$ for $`U=100`$ ($`𝒢`$ -phase), $`\lambda _{\text{core}}0.41`$, $`\lambda _{\text{halo}}0.23`$, and $`\lambda _{\text{gas}}0.36`$ for $`U=100`$ ($``$ -phase), $`\lambda _{\text{core}}0.96`$, $`\lambda _{\text{halo}}0.02`$, and $`\lambda _{\text{gas}}0.02`$ for $`U=500`$ ($`𝒞`$ -phase). Gas particles dominate in $`𝒢`$ -phase and core particles dominate in $`𝒞`$ -phase in number. While in $``$ -phase, all of the three species of particles coexist almost equally. According to our various calculations changing the particle number and initial conditions, this coexistence seems to represent the prominent property of SGS and not the finiteness of particles.

In Fig.5, the energy dependence of the ratio of each species in $`\tau =5`$ for $`ϵ=2.5\times 10^7`$ is shown at the center. The energy scale where halo particles exit corresponds to the one where negative specific heat appears ($``$ -phase). Together with it, the relative ratio of time intervals of three states in which each particle stays until $`\tau =5`$ is shown. As we expect, in both $`𝒢`$ \- and $`𝒞`$ -phase, most particles stay just in one state for quite a long time. For the $``$ -phase, however, many particles experience at least two states. Some particles wander from one state to another in three states.

## IV Particle motions

In this section, we examine individual particle motion and velocity distribution function in each phase.

### A Recurrent motion of halo particles

In the HMF model at the state near the critical energy, Lévy-type flight and anomalous diffusion of particles have been reported . These peculiar behaviors of particles are apparently caused by the transitions from core particles to gas particles and vice versa. We observe that such peculiar behaviors in HMF model are caused by the periodicity of the configuration space and not by the long range force itself. In fact, for all Lévy-type flights numerically shown in HMF, the flight distances turn out to be longer than the period $`2\pi `$. Therefore the artificial periodicity in the potential, and not the long-rage nature of the force itself, is thought to have caused the Lévy-type flight in HMF model. In our SGR model, we do have the possibility to observe the same Lévy-type flight motions since particles in SGR also move along a closed ring. Since we would like to extract the intrinsic property of the long-range nature of the gravity itself, we pay attention to the recurrent motion of halo particles and disregard the round-trip motion along the entire ring.

First, we choose the parameters as $`\xi 7.1\times 10^4R`$ and $`U=100`$ for which the system is in the intermediate energy phase. The particle motion is shown in Fig.6(a), in which core particles form a firm cluster and they oscillate around the center of the cluster with the time scale $`t_\xi `$ and halo particles go in and out of the core region without any typical time scale and amplitude. Zoomed in ten-times (Fig.6(b)), and even in hundred-times (Fig.6(c)), the particle motion is always similar recurrent movements. These repetition of the similar recurrent pattern suggests a self-similar structure of the system.

Note that the recurrent motion of halo particles is confined within the range $`[0,2\pi ]`$ and particles never experience a round-trip along the ring, quite contrary to the HMF model. Moreover, this recurrent motion of particles is quite robust and is observed in any region of the intermediate energy phase. This robustness is a remarkable contrast to the HMF model in which such motion is observed only at the critical point in the phase diagram. \**\**\** The smaller the cutoff $`ϵ`$, the larger the range of the total energy per particle $`U`$ where the recurrent motion appears.

In order to analyze this behavior more quantitatively, we examine the frequency distribution of the recurrent time $`\tau _{\mathrm{rec}}\tau _{\mathrm{in}}\tau _{\mathrm{out}}`$, that is, the time period from the moment $`\tau _{\mathrm{out}}`$ when a particle leaves the barycenter of the core to the moment $`\tau _{\mathrm{in}}`$ when it first returns to the barycenter again. Here we have defined the location of the barycenter of the core as

$$\theta _{bc}\underset{i=1}{\overset{N_c}{}}\theta _i/N_c,$$
(19)
where $`\theta _i`$ is the location of $`i`$’th particle in the core and $`N_c`$ is the total number of core particles at each moment.

In Fig.7, we depicted the frequency distribution of the recurrent time $`\tau _{\text{rec}}`$ for core particles and halo particles separately. For core particles, as we expected, we find almost Gaussian distribution around the center $`\tau _\xi `$ (Fig.7(a)). On the other hand for halo particles, although a peak is found around $`\tau _E(t_E/t_R)`$, the distribution shows that a long tail spreads widely characterized by the power law $`f(\tau _{\text{rec}})\tau _{\text{rec}}^p`$ with $`p2.0`$ (Fig.7(b)). We have checked that this power law range $`3<\mathrm{log}\tau _{\text{rec}}<1.5`$ is consistent with the eye-fitted region of self-similar motion in Fig.6. We will see that the above difference of frequency distribution for core and halo particles leads to the difference in relaxation time of them in Sec.V.

We have also examined the energy dependence of the above power $`p`$. The result is given in Fig.7(b) showing that the value of $`p`$ ($`2.0`$) is almost independent of the choice of the energy $`U`$ throughout the intermediate energy phase.

### B Velocity distribution of particles

As is shown in Appendix A, the time averaged kinetic energy $`K`$ is, in general, expressed as the sum of $`N^2`$ independent stochastic variables. However, in the case of core particles, which are located within a few cutoff scale, the quantity $`K`$ turns out to be expressed as the sum of $`N`$ independent stochastic variables (see Appendix B). In this section, we show that the $`N`$-dependence of $`K`$ characterizes the velocity distribution of particles.

In Fig.8, we show the velocity distributions of core and halo particles in three cases; (a) $`\xi 5\times 10^4R`$, $`U=100`$ (b) $`\xi 5\times 10^4R`$, $`U=500`$ and (c) $`\xi 5\times 10^2R`$, $`U=0.65`$. We have superposed the velocity data at $`\tau =1`$, $`2`$, $`3`$, $`4`$, and $`5`$, with ten different random initial conditions fixing the total energy. Thus the size of the whole data we used is $`5\times 10^3`$ for each velocity distribution function. The distribution of core particles is well fitted by the Gaussian distribution

$$P(v)=\frac{1}{\sqrt{2\pi }\sigma }e^{\frac{v^2}{2\sigma ^2}},$$
(20)
with the dispersion $`v_{\mathrm{𝑐𝑜𝑟𝑒}}^2^{1/2}\sigma =7.29(a),10.9(b)`$, and $`0.51(c)`$.

On the other hand, the distribution of halo particles in (a) is manifestly non-Gaussian, although that of the same particles in the case (c) is Gaussian. This distribution is at least unchanged until $`\tau =200`$ .

What is the origin of the above non-Gaussian distribution for halo particles in the intermediate energy phase in case (a)? Here we analyze this issue from the viewpoint of the generalized central-limit theorem. It is well known that the limit distribution for sums of independent random variables is the Gaussian distribution provided that the dispersion is finite (the central-limit theorem). However, it is less well-known that the limit distribution for sums of independent random variables is the stable distribution in general cases, including those where the dispersion is divergent. This stable distribution is defined to satisfy the relation

$$\underset{i=1}{\overset{N}{}}x_i\stackrel{d}{=}N^{1/\alpha }x,$$
(21)
where $`\stackrel{d}{=}`$ means that the distributions of both sides are equal to each other. The parameter $`\alpha `$ classifies the stable distributions and must satisfy $`0<\alpha 2`$ for the normalizability and the positivity of the probability distribution function. In the definition above, $`x_i`$ and $`x`$ are the probabilistic variables obeying the same distribution. Let us apply this to the velocity distribution of our model. We square both sides of the above equation, and obtain the non-extensive property,

$$v^2_{N\mathrm{particle}}=N^{2/\alpha 1}v^2_{\mathrm{one}\mathrm{particle}}$$
(22)
and the method of characteristic function yields the explicit form of the distribution function

$$P(v)=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}𝑑y\mathrm{exp}[ivys|y|^\alpha ].$$
(23)
Note that this stable distribution includes the Gaussian distribution (20) as a special case $`\alpha =2`$, where the kinetic energy becomes extensive, and the dispersion $`\sigma =\sqrt{2s}`$ is finite. The parameter $`s`$ is thought to be a generalized temperature.

Let us first consider the velocity distributions of core particles. We observe, from our numerical calculations, that all core particles oscillate within the narrow region of the cutoff size $`\xi `$, and the gravitational two-body interaction is dominated by the artificial potential force. Therefore the averaged kinetic energy is described as the sum of $`N`$ independent statistical elements (Appendix B). This leads to the normal extensivity for the velocity distributions of core particles and therefore we expect the Gaussian distribution. This is consistent with the results in Fig.8.

On the other hand, for halo particles which interact through genuine gravity, the kinetic energy behaves as $`K_\mathrm{p}N^2`$ for fixed $`R`$ as is shown in Appendix A. Thus from (22) we find the index $`\alpha =1`$ for the physical velocity distribution of SGR. In this case with $`\alpha =1`$, the distribution (23) becomes the Lorentz form

$$P(v)=\frac{1}{\pi }\frac{s}{v^2+s^2}.$$
(24)
This is also consistent with the results in Fig.8. The essence of the appearance of this non-Gaussian distribution is the non-extensivity of the energy for SGS; this is the intrinsic property of gravity.

Then what is the origin of the Gaussian distribution of halo particles in case (c)? As we show in the next section, this is mainly because the relaxation times of core and halo particles are close with each other. In the next section, we will discuss the relation between the profile of velocity distribution and the relaxation process.

## V Relaxation

One of the most important issues in statistical physics of the $`N`$-particle system is the relaxation process. In order to study the relaxation process of our model, we choose several different initial conditions for the same values of $`ϵ`$ and $`U`$. Here we have examined two parameter sets, i.e. (i) $`ϵ=2.5\times 10^3,U=0.65`$ and (ii) $`ϵ=2.5\times 10^7,U=100`$, for both of which negative specific heat appears.

As we discussed in Sec.II, there are three dynamical time scales: $`t_R`$, $`t_\xi `$, and $`t_E`$. In the units of our normalization, the ratios of these time scales are

$`t_E/t_R`$ $`=`$ $`{\displaystyle \frac{1}{4\sqrt{2}}}|U|^{3/2}`$ (25)
$`t_\xi /t_R`$ $`=`$ $`(2ϵ)^{3/4}.`$ (26)
For example, if we choose $`U=100`$ and $`ϵ=2.5\times 10^7`$, then $`t_E/t_R1.8\times 10^4`$ and therefore the three time scales are separated as

$$t_\xi <t_Et_R.$$
(27)
Then the contribution of gas particles to the dynamics of halo or core particles is negligible during a few $`t_R`$.

The parameter set (i) does not meet the condition(17) and $`t_E/t_R0.34`$; $`t_E`$ is the same order as $`t_R`$. On the other hand, the parameter set (ii) meets (17) and $`t_E/t_R1.8\times 10^4`$; $`t_E`$ is much shorter than $`t_R`$. This means that gas species play an important role in relaxation in case (i), while they do not in case (ii). In order to examine the relaxation process for these two cases, it is convenient to normalize the time with the unit $`t_E`$, since the degree of relaxation of 3-D gravitational system is conventionally measured with the dynamical time. So here we introduce the dimensionless time $`\tau _Et/t_E`$.

As for initial conditions, we locate several clusters of the same size with the same interval. The number of clusters is chosen as 5, 10, 20, 25, and 50 for the model (i) and 5, 10, 20, and 25 for the model (ii). The initial velocities of particles are set randomly but the total energy is fixed: $`U=0.65`$ for (i) and $`U=100`$ for (ii).

In order to study the Ergode property of the system, we examine the mixing property of the particles in the three states( core , halo , and gas ). In other words, we examine the degree of isolation of each states. We directly measure how extent each particle experiences these three states. First, we pick up discrete times $`(\tau _{E1},\tau _{E2},\mathrm{},\tau _{Ei},\mathrm{})`$ with the equal interval $`\mathrm{\Delta }\tau _E`$; at each time, each particle stays at one of the three states. Then we count how many times the $`k`$-’th particle stays at each states before the time $`\tau _{Ei}`$. Finally, after normalization, we obtain the relative frequencies of three states $`\lambda _{\text{core}}^k(\tau _{Ei}),\lambda _{\text{halo}}^k(\tau _{Ei})`$ and $`\lambda _{\text{gas}}^k(\tau _{Ei})`$ (with $`\lambda _{\text{core}}^k(\tau _{Ei})+\lambda _{\text{halo}}^k(\tau _{Ei})`$+$`\lambda _{\text{gas}}^k(\tau _{Ei})=1`$ ) for the $`k`$’th particle.

At the first time $`\tau _{E1}`$, each particle definitely stays at one of the three states; relative frequencies of stay $`\lambda _{\text{core}}^k(\tau _{E1}),\lambda _{\text{halo}}^k(\tau _{E1})`$ and $`\lambda _{\text{gas}}^k(\tau _{E1})`$ are either 1 or 0. In the next time $`\tau _{E2}`$, some particles may change the state; then $`\lambda _{\text{core}}^k(\tau _{E2}),\lambda _{\text{halo}}^k(\tau _{E2})`$ and $`\lambda _{\text{gas}}^k(\tau _{E2})`$ for these particles are $`1/2`$. In this way, $`\lambda _{\text{core}}^k(\tau _{Ei}),\lambda _{\text{halo}}^k(\tau _{Ei})`$ and $`\lambda _{\text{gas}}^k(\tau _{Ei})`$ will evolve in time $`\tau _{Ei}`$. Various distribution of $`\{\lambda _{\text{core}}^k(\tau _{Ei})\}_{k=1}^N`$ for all the particles $`k=1,2,\mathrm{},N`$ will define the distribution function $`N(\lambda _{\text{gas}},\tau _E)`$ for values $`\lambda _{gas}`$ at time $`\tau _{Ei}`$ by simply counting the number of the particles which take the value $`\lambda _{\text{core}}`$ at time $`\tau _E`$. Similarly distribution functions for other states can also be defined: $`N(\lambda _{\text{halo}},\tau _E)`$ and $`N(\lambda _{\text{core}},\tau _E)`$. As the system is thermally relaxed and the equipartition of energy is attained, $`\lambda _{\text{core}}^k(\tau _E),\lambda _{\text{halo}}^k(\tau _E)`$ and $`\lambda _{\text{gas}}^k(\tau _E)`$ will converge to the same value $`\lambda _{\text{core}}^{},\lambda _{\text{halo}}^{}`$ and $`\lambda _{\text{gas}}^{}`$ independent of the particle label $`k`$. Then $`N(\lambda _{\text{gas}},\tau _E)`$, $`N(\lambda _{\text{halo}},\tau _E)`$ and $`N(\lambda _{\text{core}},\tau _E)`$ would have sharp peaks at $`\lambda _{\text{core}}^{},\lambda _{\text{halo}}^{}`$ and $`\lambda _{\text{gas}}^{}`$ in the thermal equilibrium state. Thus the development of the peaks in the distributions functions $`N(\lambda _{\text{gas}},\tau _E)`$, $`N(\lambda _{\text{halo}},\tau _E)`$ and $`N(\lambda _{\text{core}},\tau _E)`$ can be a good measure of the degree of thermal equilibrium.

In the large UV-cutoff case (i), each distribution $`N(\lambda _{\text{gas}},\tau _E)`$, $`N(\lambda _{\text{halo}},\tau _E)`$ and $`N(\lambda _{\text{core}},\tau _E)`$ starting from $`\tau _{E1}=6.5\times 10^4`$ seems to develop a single peak in-between 0 and 1 at least by the time $`\tau _E10^5`$ (Fig.9 (b)), and the variance around the peak seems to reduce in time (Fig.9 (c)). This means that most particles experience all of the three states at least by $`\tau _E2.0\times 10^5`$; thermal relaxation proceeds in the time scale $`\tau _E10^5`$. Contrary in the small UV-cutoff case (ii), each distribution $`N(\lambda _{\text{gas}},\tau _E)`$, $`N(\lambda _{\text{halo}},\tau _E)`$ and $`N(\lambda _{\text{core}},\tau _E)`$ starting from $`\tau _{E1}=1.1\times 10^5`$ seems to develop peaks at the edge of the domain 1 and/or 0, even after the time duration, $`\tau _E3.8\times 10^5`$(Fig.9 (e)) . This means that most particles stay in each single state even in the time scale $`\tau _E10^5`$; the system does not approach thermal relaxation at all in this time scale.

By combining these results, the Gaussian velocity distribution realized in the case (i) in Fig.8 seems to reflect the achieved thermal relaxation of the system. On the other hand, the Lorentzian velocity distribution realized in the case (ii) in Fig.8 seems to reflect some quasi-equilibrium intermediate stage before the thermal relaxation of the system. Latter peculiar properties in velocity distributions may make 3-D gravitational system quite different from those in the system with short-range interactions or the system with positive specific heat. If we set $`ϵ`$ larger, then the encounter of halo particles with normal particles (gas and core particles) becomes more frequent, and the halo particles cannot be isolated from the normal particles, resulting trivial thermodynamical equilibrium of the whole system as in the case (i).

## VI Scaling properties of the HALO

In Sec.IV, we have found that the velocity distribution of halo particles is non-Gaussian with anomalous $`v^2`$-tail and each halo particle shows intermittent recursive motion around the core without a definite time scale. These scale-free properties of the halo particles turn out to appear in the spatial distribution profile of them itself. In this section, we study the scaling property of the spatial distribution with the box-counting method .

For this purpose, we use the box counting method: We divide the entire configuration space $`2\pi R`$ into segments with equal size $`\mathrm{}`$ and count the number $`N(\mathrm{})`$ of the segment which contains at least one halo particle. Then we define the quantity $`(\mathrm{log}N(\mathrm{}))/(\mathrm{log}\mathrm{})`$ which turns out to be the scaling exponent, provided that the quantity is almost independent of the scale $`\mathrm{}`$. In order to increase the statistical significance, we used multiple particle-distribution data at different times for each run. Moreover, we made many runs of calculations with different initial conditions for the fixed total energy $`U`$. We have applied the box-counting method to each data and then all the individual results are superposed. When we extract the distribution data, we have chosen the time interval $`[2000t_\xi ,5000t_\xi ]`$, in which the particle ratios $`\lambda _{\text{core}}`$, $`\lambda _{\text{halo}}`$, and $`\lambda _{\text{gas}}`$ almost reach the relaxed constant values. Thus we calculate the following averaged quantity

$$D(\mathrm{})\frac{(\mathrm{log}N(\mathrm{}))}{(\mathrm{log}\mathrm{})}_{\mathrm{ens}},$$
(28)
over all the data we thus prepared. If any scaling property exists in the particle distribution, $`D(\mathrm{})`$ would become constant for a finite range of $`\mathrm{}_1\mathrm{}\mathrm{}_2`$. As for the bounds $`\mathrm{}_1`$ and $`\mathrm{}_2`$, we have technical restrictions originating from our numerical calculation method. The averaged distance between gas particles is estimated as

$$l_{\text{gas}}(ϵ)\frac{2\pi R}{N\lambda _{\text{gas}}}.$$
(29)
If the box size $`\mathrm{}`$ is larger than $`l_{\text{gas}}(ϵ)`$, we expect $`D1.0`$, because almost all the boxes contain at least one particle. On the other hand, if $`\mathrm{}`$ is smaller than the cutoff distance $`\xi `$, the genuine property of gravity is lost. Therefore scaling property in particle distribution is relevant only within the rage $`\xi \mathrm{}l_{\text{gas}}`$. Actually in our calculation, for $`ϵ=2.5\times 10^7`$ and $`U=100`$, we can see the typical behavior of $`N(\mathrm{})`$ as a function of $`\mathrm{}`$ in Fig.10. In this plot, we can see scaling behavior $`N(\mathrm{})\mathrm{}^d`$ with the small exponent in the range $`10^3R<\mathrm{}<10^1R`$, which is well inside the relevant region. This scaling seems to originate from halo particles, because they not only dominate in the above scaling range but also show non-Gaussian velocity distributions and self-similar recursion jumps.

However, the scaling property of the box counting method does not always conclude the existence of the fractal structure in the particle distributions. This is because the box counting method itself cannot distinguish the two possibilities; a) genuine fractal structure and b) the power-law tail of the particle distributions around the core center.

In order to distinguish the above possibilities, we compare two different superposition methods as follows: a) We simply superpose all the data with the bare coordinate $`\theta `$ and b) we superpose all the data with the coordinate adjusting so that the mass-center of each data comes to the same position; i.e. we introduce a new coordinate $`\overline{\theta }`$ as $`\overline{\theta }\theta \theta _{bc}`$ for each data with the mass center $`\theta _{bc}`$.

As a result, the scaling property of the particle distributions, which was observed in the original method, disappears completely in the case a). On the other hand in the case b) we observe the same scaling property as was observed in the original method. This result suggests that the observed scaling property in particle distribution is not the genuine fractal property of the system, but simply due to the power-law distribution of particles around the mass-center. We have confirmed this with our further analysis on the direct observation of the power-law particle distributions around the mass-center.

The energy and cutoff dependence on the scaling exponent $`D`$ is shown in Fig.11. Here we show the scaling exponent of halo particles in the region of negative specific heat. We find that the scaling exponent $`D`$ is almost constant in this region. Moreover, smaller cutoff $`ϵ`$ corresponds to the smaller dimension $`D`$. On the other hand in the limit of $`ϵ1`$, the scaling property cannot be observed.

In the limit of real gravitational interaction ($`ϵ0`$), the region of negative specific heat extends in energy range and so does the scaling region. This constant property of the exponent $`D`$ is likely to be independent of the temperature or the energy of the system.

## VII Conclusion

In this paper we have studied the self-gravitating ring (SGR) model which is a one-dimensional (1-D) model of a self-gravitating system. Although the system is 1-D, the particles constrained on a circular ring mutually interact with the 3-D full gravitational force. We found that an interesting phase with the negative specific heat appears at the intermediate energy scale, which reflects the virial condition of the 3-D gravitating system. Classifying the particles in each phase into three species (core , halo , and gas specie), we examined dynamical properties of this phase. Using SGR model, we further examined characteristic properties of relaxation, velocity distributions, and density profiles in each phase.

The cutoff parameter $`\xi `$ introduced in SGR model controls the singular properties of the system at short-distances. By changing this parameter, we could systematically study the effect of the short-distance singularity of the force upon the statistical properties of the system.

When the cutoff parameter is of order one ($`\xi R`$ ), the system resembles the Hamilton Mean Field model (HMF). In fact, the velocity distribution becomes Gaussian at $`t𝒪(10^3)t_R`$. This is because the two time scales $`t_E`$ and $`t_R`$ become the same order and therefore halo particles frequently encounter and interact with gas particles. This strong interaction induces frequent exchange of particles between these two species. Thus the energy diffuses effectively and complete thermodynamical equilibrium of the entire system is established. Especially the temperature is the same for all the particle species. In other words, the negative specific heat region (halo particles) does exist, but the strong interaction with gas and core particles (positive specific heat) quickly dissolves the negative specific heat region.

In much more smaller cutoff case with $`\xi 7.1\times 10^2R(ϵ=2.5\times 10^3)`$, the velocity distribution becomes almost Gaussian after the virialization, although the temperature of each species is different from each other. This temperature difference, as well as the required time for establishing the equilibrium, is increased for smaller value of the cutoff $`\xi `$.

On the other hand, as $`\xi 0`$, the difference of the two time scales increases $`t_Rt_\xi `$, and the energy region of the negative specific heat extends wider. Therefore we can expect that halo particles become almost isolated from the gas particles. This isolation of unstable halo particles from other normal species enables the halo species to last sufficiently long. In the case of $`\xi 7.1\times 10^4R(ϵ=2.5\times 10^7)`$, the initial condition dependence of halo particles survives even after $`5\times 10^4t_R`$ in our simulation, which prevents the full system from reaching thermal equilibrium. These results resemble the case of core particles with $`\xi 7.1\times 10^2R`$. However, the remarkable difference between the two cases is that the velocity distribution of halo particles in $`\xi 7.1\times 10^4R`$ is not Gaussian but Lorentzian distribution, while that of core particles in $`\xi 7.1\times 10^4R`$ is Gaussian. This property of halo particles qualitatively reflects genuine 3-D gravitational interaction and is the main difference from the HMF models or 1DS models We have also found the scaling structure of halo particles in the intermediate energy phase. When we reduce the cutoff $`\xi 0`$, the scaling region as well as the energy range of negative specific heat is extended wider, although the value of exponent depends on the cutoff and reduces with $`\xi `$.

In 1DS model, the system reaches thermodynamical equilibrium at least at $`t𝒪(10^7)t_E`$, where $`t_E(1/4\pi GmN)\sqrt{4E/mN}`$ is the crossing time of 1DS model. It is also shown that the system reaches the collisionless mixing phase during much shorter time interval. Our simulation with $`\xi 7.1\times 10^4R`$ shows non-Gaussian velocity and scaling property in such a short time interval $`4\times 10^3t_E2\times 10^5t_E`$ for $`U=100`$. So the non-Gaussian velocity distribution might reflect the character in collisionless mixing phase which appears in 1DS simulation . However, it might also be true that the relaxation process in SGR is quite different from that in 1DS, where no time scale separation between core and halo particles exists. So the exotic character of halo particles might be intrinsic in SGR or the system with real 3-D gravitational interactions because of the time scale separation.

In summary, SGR model is characterized (i) by the phase separation and the particle state separation and (ii) by similar recurrent motion and non-Gaussian velocity distribution of halo particles . Both of these characters may play a key role in explaining the observational properties in real 3-D self-gravitating objects. For example, the recent observations with Hubble Space Telescope support the existence of supermassive black holes or cusps at the center of triaxial elliptical galaxies . This observation suggests that the above property (i) may appear in these galaxies, since they have high density region corresponding to core particles at their centers and low density region corresponding to halo particles surrounding them. Hence, the property (i) may give a hint for explaining the relaxation process and the origin of stationally configuration shape of these galaxies, since the stochastic halo orbits which interact with core particles at the center of a triaxial galaxy may affect the equilibrium shape through the continued mixing near the center . In addition, the propety (ii) might help to explain the observed fractal structures and non-Gaussian velocity distributions in the inter stellar medium, which are likely to be gravitationally virialized.

###### Acknowledgements.

We would like to thank A. Nakamichi, I. Joichi, K. Nakamura, and M. Hotta for useful discussions and comments. This work was supported partially by a Grant-in-Aid for Scientific Research Fund of the Ministry of Education, Science and Culture (Specially Promoted Research No. 08102010), and by the Waseda University Grant for Special Research Projects.

## A Virial condition for SGR model

We start with the non-dimensional Hamiltonian of the SGR model,

$$H=\frac{1}{2}\underset{i=1}{\overset{N}{}}p_i^2\underset{i<j}{\overset{N}{}}\frac{1}{\sqrt{2}N\sqrt{1\mathrm{cos}\theta _{ij}+ϵ}}.$$
(A1)
Canonical equations of motion for this system become

$`{\displaystyle \frac{d\theta _i}{d\tau }}`$ $`=`$ $`p_i,`$ (A2)
$`{\displaystyle \frac{dp_i}{d\tau }}`$ $`=`$ $`{\displaystyle \underset{ji}{\overset{N}{}}}{\displaystyle \frac{\mathrm{sin}\theta _{ij}}{2\sqrt{2}N(1\mathrm{cos}\theta _{ij}+ϵ)^{3/2}}}.`$ (A3)
Then the second time derivative of the inertia moment of this system becomes

$`{\displaystyle \frac{d^2}{d\tau ^2}}{\displaystyle \underset{i=1}{\overset{N}{}}}\theta _i^2`$ $`=`$ $`2\left({\displaystyle \underset{i=1}{\overset{N}{}}}p_i{\displaystyle \frac{d\theta _i}{d\tau }}+{\displaystyle \underset{i=1}{\overset{N}{}}}\theta _i{\displaystyle \frac{dp_i}{d\tau }}\right)`$ (A4)
$`=`$ $`2{\displaystyle \underset{i=1}{\overset{N}{}}}p_i^22{\displaystyle \underset{i<j}{\overset{N}{}}}{\displaystyle \frac{\theta _i\mathrm{sin}\theta _{ij}}{\sqrt{2}N(1\mathrm{cos}\theta _{ij}+ϵ)^{3/2}}}.`$ (A5)
Taking long-time average for this equation, we obtain

$$\frac{d^2}{d\tau ^2}\underset{i=1}{\overset{N}{}}\theta _i^2=4NT2\underset{i<j}{\overset{N}{}}\frac{\theta _i\mathrm{sin}\theta _{ij}}{\sqrt{2}N(1\mathrm{cos}\theta _{ij}+ϵ)^{3/2}}.$$
(A6)
Therefore if the condition

$$I(\tau )\frac{d^2}{d\tau ^2}\underset{i=1}{\overset{N}{}}\theta _i^2/(NT)1,$$
(A7)
is realized, then the virial relation

$$T\frac{1}{2N}\underset{i<j}{\overset{N}{}}\frac{\theta _i\mathrm{sin}\theta _{ij}}{\sqrt{2}N(1\mathrm{cos}\theta _{ij}+ϵ)^{3/2}},$$
(A8)
holds.

By using the total potential

$$V=\underset{i<j}{\overset{N}{}}\frac{1}{\sqrt{2}N\sqrt{1\mathrm{cos}\theta _{ij}+ϵ}},$$
(A9)
the virial relation (A8) can be expressed as

$$T\frac{1}{2N}\underset{i<j}{\overset{N}{}}\theta _i\frac{V}{\theta _i}.$$
(A10)
Note that even in the limit of $`ϵ0`$, the virial relation $`K=2NT=V/2`$ is not justified, because of the potential form (A9). However, as is shown in Sec.II, we can get the above ideal virial condition for 3-D gravity within the appropriate energy range.

## B Virial condition in a $`𝒞`$ -phase for SGR model

In $`𝒞`$ -phase, almost all of the particles are trapped inside the core , whose size is almost $`\xi `$. Therefore in this state canonical equations (A3) are approximated as

$`{\displaystyle \frac{d\theta _i}{d\tau }}`$ $`=`$ $`p_i,`$ (B1)
$`{\displaystyle \frac{dp_i}{d\tau }}`$ $`=`$ $`{\displaystyle \underset{ji}{\overset{N}{}}}{\displaystyle \frac{\theta _{ij}}{2\sqrt{2}Nϵ^{3/2}}}`$ (B2)
$`=`$ $`{\displaystyle \frac{\theta _i}{2\sqrt{2}ϵ^{3/2}}}\left(1{\displaystyle \frac{1}{N}}\right)+{\displaystyle \underset{ji}{\overset{N}{}}}{\displaystyle \frac{\theta _j}{2\sqrt{2}Nϵ^{3/2}}},`$ (B3)
since $`|\theta _{ij}|𝒪(ϵ)1`$.

In this approximation, the second time derivative of the inertia moment of this system becomes

$`{\displaystyle \frac{d^2}{d\tau ^2}}{\displaystyle \underset{i}{}}\theta _i^2`$ $`=`$ $`2{\displaystyle \underset{i}{}}p_i^2{\displaystyle \frac{1}{\sqrt{2}ϵ^{3/2}}}\left(1{\displaystyle \frac{1}{N}}\right){\displaystyle \underset{i=1}{\overset{N}{}}}\theta _i^2+{\displaystyle \frac{2}{\sqrt{2}Nϵ^{3/2}}}{\displaystyle \underset{i<j}{\overset{N}{}}}\theta _i\theta _j.`$ (B4)
Taking long-time average for this relation, we obtain

$$\frac{d^2}{d\tau ^2}\underset{i}{}\theta _i^2=4NT\frac{1}{\sqrt{2}ϵ^{3/2}}\left(1\frac{1}{N}\right)\underset{i=1}{\overset{N}{}}\theta _i^2+\frac{2}{\sqrt{2}Nϵ^{3/2}}\underset{i<j}{\overset{N}{}}\theta _i\theta _j.$$
(B5)
Since $`\theta _i`$ and $`\theta _j(ij)`$ are independent stochastic variables, the last term of (B5) vanishes. Hence from the virial relation,

$$T\frac{1}{4\sqrt{2}NNϵ^{3/2}}\left(1\frac{1}{N}\right)\underset{i=1}{\overset{N}{}}\theta _i^2,$$
(B6)
and

$$K_\mathrm{p}\frac{Gm^2}{2\sqrt{2}ϵ^{3/2}R}\left(1\frac{1}{N}\right)\underset{i=1}{\overset{N}{}}\theta _i^2,$$
(B7)
hold. This shows that $`K_\mathrm{p}`$ in $`𝒞`$ -phase is given by the sum of $`N`$ statistical elements.


format: text

text:
# P–vortices and Drama of Gribov Copies. ITEP-TH-50/00

## 1 Introduction

The idea that the center vortices are the objects responsible for confinement in the nonabelian gauge theories is rather old . Recently it has been argued that the center projection might provide a powerful tool to investigate this idea . It is suggested that projection dependent P-vortices defined on the lattice plaquettes are able to locate thick gauge invariant center vortices and thus provide the essential evidence for the center vortex picture of confinement. So far 3 different center gauges have been used in practical computations: the indirect maximal center (IMC) gauge , the direct maximal center (DMC) gauge  and the Laplacian center gauge  (see also  for a new proposal).

IMC and DMC gauges suffer from the Gribov copies problem  and this technical problem gave rise to several claims and counterclaims. At present we have a drama in five acts:

i) Initially there was the claim  that the projected string tension (the string tension which is due to P-vortices) reproduces the full $`SU(2)`$ string tension. Thus P-vortices are responsible for the confinement of color in gluodynamics.

ii) The problem of Gribov copies has been raised in Ref. , where it has been demonstrated that there are gauge copies which produce P-vortices evidently with no center vortex finding ability since projected Wilson loops have no area law. At the same time these gauge copies correspond to higher maxima of the gauge fixing functional, than those used in .

iii) In our previous publications  we confirmed the existence of two classes of gauge copies: with nonzero and zero projected string tension $`\sigma _{Z(2)}`$. We resolved the problem raised in  since we found the copies with the highest maxima of the gauge fixing functional which correspond to nonzero $`\sigma _{Z(2)}`$. At the same time the value of $`\sigma _{Z(2)}`$ was essentially lower than that obtained in  and in is disagreement with the physical string tension $`\sigma _{SU(2)}`$.

iv) In Refs.  it has been argued that the disagreement between $`\sigma _{Z(2)}`$ and $`\sigma _{SU(2)}`$ is due to strong finite volume effects and it disappears on large enough lattices. In Refs.  the usual relaxation plus overrelaxation (RO) algorithm was used to fix the DMC gauge.

v) Below we present results111The preliminary results are presented in Ref. , obtained with the help of more powerful gauge fixing algorithm, simulated annealing (SA), and show that the problem of low value of $`\sigma _{Z(2)}`$ persists even on large lattices with physical extension up to $`3fm`$. We discuss here DMC gauge. Our results for IMC gauge will be presented elsewhere.

## 2 Gauge fixing procedure

### 2.1 Gribov copies

Direct maximal center gauge  in $`SU(2)`$ lattice gauge theory is defined by the maximization of the following functional:

$$F(U)=\frac{1}{4V}\underset{n,\mu }{}\left(\frac{1}{2}\text{Tr}U_{n,\mu }\right)^2=\frac{1}{4V}\underset{n,\mu }{}\frac{1}{4}\left(\text{Tr}_{adj}U_{n,\mu }+1\right),$$
(1)
with respect to local gauge transformations, $`U_{n,\mu }`$ is the lattice gauge field, $`V`$ is the lattice volume. Condition (1) fixes the gauge up to $`Z(2)`$ gauge transformation, and can be considered as the Landau gauge for adjoint representation. Any fixed configuration can be decomposed into $`Z(2)`$ and coset parts: $`U_{n,\mu }=Z_{n,\mu }V_{n,\mu }`$, where $`Z_{n,\mu }=\text{sign}\text{Tr}U_{n,\mu }`$. The plaquettes $`Z_{n,\mu \nu }`$ constructed from the links $`Z_{n,\mu }`$ have values $`\pm 1`$. The P-vortices (which form closed surfaces in 4D space) are made from the plaquettes, dual to plaquettes with $`Z_{n,\mu \nu }=1`$.

Some evidence has been collected, that P-vortices in DMC gauge can serve to locate gauge invariant center vortices. It has been reported  that the expectation value of the projected Wilson loops $`W_{Z(2)}`$, computed via linking number of the static quarks trajectories and P-vortices, have area law. The projected string tension $`\sigma _{Z(2)}`$ is very close to the string tension of the nonabelian theory $`\sigma _{SU(2)}`$. This fact has been called center dominance. Another important observation was that the density of P-vortices scales as a physical quantity .

The main problem of DMC gauge fixing is that the functional $`F(U)`$ (1) has many local maxima. This is the analogue of the Gribov problem in continuum gauge theories  and we call configurations corresponding to these local maxima Gribov copies. For gauge conditions with Gribov copies the usual expression for gauge dependent quantities  does not provide unambiguous definition (see recent discussion of this subject in ). It was shown that for some gauge conditions the gauge dependent quantities depend strongly on the local maxima picked up . Thus to remove this ambiguity it is necessary to find the global maximum or, if this is impossible, to approach the global maximum as close as possible222For another proposal to solve this problem see e.g. ..

The usual algorithm for the gauge of the type (1) is the relaxation algorithm which performs maximization iteratively site by site. The relaxation becomes more effective if one uses the overrelaxation. It is also known that the SA algorithm is more effective than RO algorithm. It is very useful when ambiguities induced by Gribov copies become severe . For this reason we use this algorithm for DMC gauge.

We follow procedure proposed and checked in : for given configuration we generate $`N_{cop}`$ gauge equivalent copies applying random gauge transformations to the initial configuration, after that we fix the gauge for each gauge copy thus producing $`N_{cop}`$ Gribov copies for each configuration. Then we compute a gauge dependent quantity $`X`$ on the configuration corresponding to the Gribov copy with the highest local maximum of (1), $`F_{max}(N_{cop})`$. Averaging over statistically independent gauge field configurations and varying $`N_{cop}`$ we obtain the function $`<X(N_{cop})>`$ and extrapolate it to $`N_{cop}\mathrm{}`$ limit. This should provide a good estimation for $`<X>`$ computed on the global maximum unless the algorithm in use produces local maximum far from the global ones.

Note that there exists another proposal  for the gauge fixing procedure which is free of gauge copies problem. In some particular limit this procedure corresponds to the search of the global maximum. There is also a class of gauge conditions  which do not suffer from the gauge copies problem.

### 2.2 Simulated annealing algorithm

The functional $`F(U)`$ can be considered  as the spin action:

$$S(s)=F(U^g)=\frac{1}{16V}\underset{n,\mu }{}\text{Tr}_{adj}\left(s_nU_{n,\mu }s_{n+\widehat{\mu }}+1\right)$$
(2)
where the “spin variable” $`s_n`$, is the gauge transformation matrix in the adjoint representation: $`s_{n,ab}=\frac{1}{2}tr(\sigma _ag_n\sigma _bg_n^{})`$ .

The lattice field $`U_{n,\mu }`$ plays the rôle of (almost) random local couplings. Maximization of the functional $`F(U^g)`$ is equivalent to decreasing the auxiliary temperature $`T`$ of the statistical system with the partition function,

$$Z=\underset{\{s_n\}}{}\mathrm{exp}\left(\frac{1}{T}S(s)\right).$$
(3)
One starts with equilibrating this spin glass at high temperature. Subsequently the temperature, $`T`$, is decreased adiabatically. In the limit $`T0`$ the system approaches the ground state, i.e. the maximal value of $`S(s)`$. Our SA implementation consists of the following three steps:

1. Thermalization at $`T=1.16`$.
2. Gradual decreasing of $`T`$ down to $`T=0.02`$.
3. Final maximization by means of the RO algorithm.

In steps 1 and 2 an updating of the spin configuration was done with Metropolis algorithm. The temperature was lowered after every sweep by the quantity $`\delta T(T)`$. $`\delta T(T)`$ has been tuned such that the spin action increased about linearly with the number of iteration sweeps. This has been realized by subdividing full range \[0.02,1.16\] into 38 intervals of width $`\mathrm{\Delta }T=0.03`$. The corresponding differences of the average action $`\mathrm{\Delta }\overline{S}(T)>=\overline{S}(T)\overline{S}(T\mathrm{\Delta }T)`$ have been computed on equilibrated configurations and were found to be very stable against statistical fluctuations among different Monte Carlo configurations. We used the same $`\delta T(T)`$ for all volumes at given value of gauge coupling $`\beta `$ and slightly modified it for different $`\beta ^{}s`$. The number of sweeps to be performed within each interval $`(T\mathrm{\Delta }T,T]`$ was chosen to be proportional to $`\mathrm{\Delta }\overline{S}(T)`$ and, subsequently, the corresponding value of $`\delta T(T)`$ has been determined. To approach close enough to the equilibrium we have to perform about 2000 sweeps at the step 2 (for step 1 only 20 sweeps was enough). Finally, at step 3 it is necessary to make $`O(100)`$ sweeps of the RO algorithm to satisfy convergence criterion. In total this procedure was essentially more costly than the RO algorithm, also it gives the higher maximum of the functional (1).

## 3 Results

### 3.1 Details of simulation

Our computations have been performed on lattices $`12^4,16^4,20^4`$ and $`24^4`$ for $`\beta =2.4`$ and $`16^4,20^4`$ for $`\beta =2.5`$. The number of configurations was $`150`$ for $`12^4,`$ $`50`$ for $`16^4,`$ $`30`$ for $`20^4`$, $`20`$ for $`24^4`$.

Using the described above gauge fixing procedure we calculate various observables as functions of the number of randomly generated gauge copies $`N_{cop}`$. To make computations feasible we reduce the maximal value of $`N_{cop}`$ on large lattices, see Table 1.

The projected Wilson loops $`W_{Z(2)}(C)`$ are defined as

$$W_{Z(2)}(C)=<\mathrm{exp}\{i\pi (\mathrm{\Sigma }_P,C)\}>$$

where $`(\mathrm{\Sigma }_P,C)`$ is the 4D linking number of the closed surface, $`\mathrm{\Sigma }_P`$, formed by P-vortex and closed loop $`C`$. We used two ways to estimate projected string tension $`\sigma _{Z(2)}`$: from the Creutz ratios

$$\chi _{Z(2)}(I,I)=\mathrm{log}\frac{W_{Z(2)}(I,I)W_{Z(2)}(I+1,I+1)}{W_{Z(2)}(I,I+1)W_{Z(2)}(I+1,I)};$$
(4)
and fitting static potential

$$V_{Z(2)}(R)=\underset{T\mathrm{}}{lim}\mathrm{log}\frac{W_{Z(2)}(R,T+1)}{W_{Z(2)}(R,T)}$$

to the usual form

$$V_0+c/R+\sigma _{Z(2)}a^2R.$$
(5)
The P-vortex density $`\rho `$ is another quantity of interest. It is defined as follows

$$\rho =\frac{1}{12V}\underset{n;\mu >\nu }{}(1Z_{n,\mu \nu }).$$

### 3.2 Maximized functional and extrapolation to $`N_{cop}\mathrm{}`$

In DMC gauge the average value $`<F_{max}>`$ is an important indicator of the quality of the gauge fixing since gauge noninvariant observables we are measuring are strongly correlated with $`<F_{max}>`$. In Fig.1 we plot $`<F_{max}>`$ as a function of $`N_{cop}`$. For comparison we show the results for RO algorithm and the algorithm suggested in  : before applying RO algorithm the random gauge copy was fixed to Landau gauge. We call it LRO algorithm. One can see that SA algorithm indeed produces highest maxima. We checked that this is true for other lattices as well.

One can draw two conclusions from Fig.1. First, it is clear that the local maxima generated with LRO procedure of  can be safely ignored. Second, RO algorithm does not permit to reach the global maxima and results obtained with this algorithm must be revised with the help of SA algorithm. As we mentioned the gauge dependent quantities are strongly correlated with the value of $`F_{max}`$.

The important question is how close we approach the global maximum. Theoretically SA algorithm permits to find the global maximum. In practice this is complicated task because of large autocorrelation length at small temperature $`T`$. We checked for a few configurations that increasing the number of sweeps in our SA algorithm by factor 5 does not lead to any essential change in $`F_{max}`$. This indicates that we are close to the equilibrium during SA cooling and thus our results are close to the global maximum when we make the extrapolation to $`N_{cop}\mathrm{}`$ limit.

As we show in the Appendix there are 2 possibilities for dependence of $`<F_{max}>`$ on $`N_{cop}`$. Our data show that at small $`N_{cop}`$ eq.(A.2) provides better fit while at large $`N_{cop}`$ eq.(A.1) is better. For this reason we use eq.(A.1) to extrapolate the data to $`N_{cop}\mathrm{}`$ limit. In Fig.2 we depict $`\chi _{Z(2)}(5,5)`$ and $`<\rho >`$ as functions of $`<F_{max}>`$. The linear dependence seen in this figure proves that the use of the same fitting function (eq. (A.1)) for $`<F_{max}>`$ and for $`\chi _{Z(2)}(I,I)`$ and $`<\rho >`$ is indeed well grounded.

### 3.3 P-vortex density

The physical density of P-vortices $`<\rho /a^2>`$ is the important indicator of their relevance for the infrared physics. If, e.g., it goes to zero in the continuum limit, then P-vortices cannot play any important role. It is known that for the completely uncorrelated vortices $`<\rho /a^2>=\frac{1}{2}\sigma _{Z(2)}`$ .

To compare the scaling of density with that of the unprojected string tension we show in Fig.3 the ratio333$`\sigma _{SU(2)}a^2=0.0728(6)`$ for $`\beta =2.4\text{[23]},0.0350(4)`$ for $`\beta =2.5`$  $`<\rho >/\sigma _{SU(2)}a^2`$. One can see that for two $`\beta `$ values this ratio is almost the same, i.e. scaling properties of $`<\rho /a^2>`$ and $`\sigma _{SU(2)}`$ are similar. The volume dependence is rather weak. Our value for $`<\rho >`$ is essentially lower than that obtained in  and this difference increases with $`\beta `$. It will be clear from our results in the next section that $`<\rho /a^2>\sigma _{Z(2)}`$, i.e. P-vortices are far from being uncorrelated. On the other hand (at least on the lattices up to $`24^4`$) the density of P-vortices does not reproduce the full string tension, $`\sigma _{SU(2)}`$, as it was suggested in Ref. .

### 3.4 Projected string tension $`\sigma _{Z(2)}`$

#### 3.4.1 Creutz ratio $`\chi _{Z(2)}`$

At first we discuss the final volume effects. In Fig.4 we show the dependence of $`\chi _{Z(2)}(4,4)`$ on $`N_{cop}`$ for $`\beta =2.4`$. From the figure one can see that finite volume effects are essential for $`L=12`$ only. For the other values of $`I`$ we obtained the similar behavior. At $`\beta =2.5`$ we found a sizable finite volume effects only on $`L=16`$ lattice.

Usually one uses the Creutz ratio $`\chi _{Z(2)}(I,I)`$ (4) to estimate the projected string tension because the dependence on $`I`$ was found to be rather weak. In our study we found that it was not quite true. In Fig.5 we plot the ratio $`\chi _{Z(2)}(I,I)/\sigma _{SU(2)}a^2`$ as a function of the dimensional distance 444The physical distance scale is set by the value $`\sigma _{SU(2)}=(440Mev)^2`$ $`r=Ia`$. The data shown were obtained for maximal values of $`N_{cop}`$ given in Table 1. The values extrapolated to $`N_{cop}\mathrm{}`$ limit are only slightly lower. One can see that $`\chi _{Z(2)}(I,I)`$ is not constant but increases with $`I`$ at least in the range $`0.3fmIa0.7fm`$. The similar behavior can be seen in Fig.3 of . We found that this effect weakens, when $`N_{cop}`$ increases, but does not disappear completely. It is evident from this figure that even on our largest lattice ($`L=24`$ at $`\beta =2.4`$), which has extension in physical units of about $`3fm`$, $`\chi _{Z(2)}(I,I)`$ is much lower than the unprojected string tension. This invalidates the suggestion of  that $`\chi _{Z(2)}(I,I)`$ approaches $`\sigma _{SU(2)}a^2`$ on lattices larger than the size of $`Z(2)`$-vortex ($`1fm`$). Let us remind that the authors of Ref.  have used the RO algorithm to fix the gauge, while we are using the SA algorithm.

Our estimation for $`\sigma _{Z(2)}/\sigma _{SU(2)}`$ made from projected Creutz ratio measured at distances $`0.70.9fm`$ is 0.66(2). Comparing our data obtained at given $`\beta `$ on lattices of various size (see Fig.5), we can conclude that for lattices of physical extension larger than $`1.7fm`$ there are no sizable finite volume effects for the string tension extracted from the Wilson loops of the size $`1fm`$. Still we cannot exclude that the projected string tension, measured from projected Wilson loops with a very large extension, approaches the full string tension.

#### 3.4.2 Static potential $`V_{Z(2)}(R,T)`$

In view of the dependence of $`\chi _{Z(2)}(I,I)`$ on $`I`$ it is interesting to study the projected static potential $`V_{Z(2)}(R)=lim_T\mathrm{}V_{Z(2)}(R,T),V_{Z(2)}(R,T)=\mathrm{log}\frac{W_{Z(2)}(R,T+1)}{W_{Z(2)}(R,T)}`$. As it is seen from Fig.6 $`V_{Z(2)}(R,T)`$ rises as a function of $`T`$. If we assume that for $`W_{Z(2)}(R,T)`$ the usual expansion

$$W_{Z(2)}(R,T)=C_0(R)exp\left(V_{Z(2)}(R)T\right)+C_1(R)exp\left(V_1(R)T\right)+\mathrm{}$$
(6)
is valid, then such behavior implies that $`C_1(R)<0`$.

To compute $`V_{Z(2)}(R)`$ we first extrapolated $`V_{Z(2)}(R,T)`$ as a function of $`N_{cop}`$ to $`N_{cop}\mathrm{}`$ limit and then evaluated $`V_{Z(2)}(R)`$ using large $`T`$ values. The results are depicted in Fig.7.

The value of $`\sigma _{Z(2)}`$ obtained from the projected static potential is in agreement with that obtained from the projected Creutz ratio (Section 3.4.1).

## 4 Conclusions

We made the most careful up to date gauge fixing of DMC gauge employing SA algorithm. We have shown that this algorithm permits to get the higher local maximum of the functional $`F(U)`$ than that obtained before with other algorithms. Though increment in $`<F_{max}>`$ looks not impressive it gives rise to substantial change in the values of gauge noninvariant observables. This, in turn, brought us to conclusions different from those made in . We computed the projected string tension $`\sigma _{Z(2)}`$ and P-vortex density $`<\rho >`$ on volumes with physical extension from $`1.4fm`$ up to $`3fm`$. Our results imply that only on our smallest lattices ($`L=12`$ for $`\beta =2.4`$ and $`L=16`$ for $`\beta =2.5`$) there are moderate finite volume effects. Our data for lattice size larger than $`1.7fm`$ show the absence of the finite volume effects.

The ratio of P-vortex density to $`\sigma _{SU(2)}`$ weakly depends on $`\beta `$ (see Fig.3) and the early asymptotic scaling of this density (observed in ) does not exist.

The projected string tension, $`\sigma _{Z(2)}`$, does not reproduce the full string tension, $`\sigma _{SU(2)}`$. The obtained value of the ratio $`\sigma _{Z(2)}/\sigma _{SU(2)}0.66(2)`$ is rather low. Comparing results for two lattice spacings we do not see any convergence of this ratio to 1 in the continuum limit (further studies at larger values of $`\beta `$ on larger lattices are important). Let us remind that abelian monopoles are responsible for $`92\%94\%`$ of the full $`SU(2)`$ string tension .

Since the P-vortices in DMC gauge do not reproduce the full $`SU(2)`$ string tension they cannot be responsible for confinement. On the other hand the fact that $`\sigma _{Z(2)}`$ is of the same order as $`\sigma _{SU(2)}`$ probably means that the DMC gauge is close to some “ideal gauge” which detect perfectly  the thick gauge independent center vortices . Recent results  show that an example of such gauge is the Laplacian Center Gauge which is free of Gribov copies problem and $`\sigma _{Z(2)}\sigma _{SU(2)}`$ in this gauge.

## Acknowledgements

We are grateful to V. Anikeeva and A. Veselov for the help at some stages of this work. We thank T. Kovacs for the helpful remark. This work was partially supported by grants RFBR 99-01230a, INTAS 96-370, Monbushu grant and CRDF award RP1-2103.

## Appendix A

Let the probability distribution of the random variable $`x`$ is $`p(x)`$ and let $`p(x)`$ is nonzero on the interval $`[x_0\delta _1,x_0+\delta _2]`$. Let us introduce

$$P(x)=_{x_0\delta _1}^x𝑑yp(y)$$

with properties

$$P^{}(x)=p(x),P(x_0\delta _1)=0,P(x_0+\delta _2)=1.$$

Then probability distribution $`P^{(N)}(x)`$ to get value x as maximal out of N attempts is equal to

$$P^{(1)}(x)=p(x)$$

$$P^{(N)}(x)=C_NP(x)^{N1}p(x)$$

where normalization constant

$$C_N^1=_{x_0\delta _1}^{x_0+\delta _2}𝑑xP(x)^{N1}p(x)=_{P(x_0\delta _1)}^{P(x_0+\delta _2)}𝑑PP^{N1}=\frac{1}{N}P^N|_0^1=\frac{1}{N}$$

We have to evaluate

$$\overline{x}^{(N)}=_{x_0\delta _1}^{x_0+\delta _2}𝑑xxP^{(N)}(x)=N_{x_0\delta _1}^{x_0+\delta _2}𝑑xxp(x)P^{N1}(x)$$

$$=x_0+N_{\delta _1}^{\delta _2}𝑑yyp(y)P^{N1}(y)(x=y+x_0)$$

$$=x_0+_{\delta _1}^{\delta _2}𝑑yy\frac{dP^N(y)}{dy}=x_0+yP^N(y)|_{\delta _1}^{\delta _2}_{\delta _1}^{\delta _2}𝑑yP^N(y)$$

$$=x_0+\delta _2_{\delta _1}^{\delta _2}𝑑yP^N(y)=\overline{x}_{\delta _1}^{\delta _2}𝑑yP^N(y)$$

where $`\overline{x}=x_0+\delta _2`$ is evidently the value of $`\overline{x}^{(N)}`$ in the limit $`N\mathrm{}`$. The integral in the last equation can be evaluated for large $`N`$ in the following way:

$$I=_{\delta _1}^{\delta _2}𝑑yP^N(y)=_{\delta _1}^{\delta _2}𝑑ye^{Nlog(P(y))}$$

We expand the function $`f(y)=log(P(y))`$ at its maximum $`y=\delta _2`$:

$$f(y)=\frac{P^{}(\delta _2)}{P(\delta _2)}ϵ+\frac{1}{2}[\frac{P^{\prime \prime }(\delta _2)}{P(\delta _2)}(\frac{P^{}(\delta _2)}{P(\delta _2)})^2]ϵ^2+\mathrm{}=p(\delta _2)ϵ+\frac{1}{2}[p^{}(\delta _2)p^2(\delta _2)]ϵ^2+\mathrm{},ϵ=y\delta _2$$

Then

$$I=_{\delta _1\delta _2}^0𝑑ϵe^{N\left(p(\delta _2)ϵ+\frac{1}{2}[p^{}(\delta _2)p^2(\delta _2)]ϵ^2+\mathrm{}\right)}$$

Let us consider 2 cases:

* $`p(\delta _2)0`$.

$$I=\frac{1}{N}_{N(\delta _1+\delta _2)}^0𝑑ϵe^{p(\delta _2)ϵ+\frac{1}{2N}[p^{}(\delta _2)p^2(\delta _2)]ϵ^2+\mathrm{}}=\frac{1}{Np(\delta _2)}\left(1+O\left(\frac{1}{N}\right)\right)$$

Thus in this case the fitting function for numerical data should be:

$$f_1(N)=\overline{x}+\frac{C_1}{N}+\frac{C_2}{N^2}+\mathrm{}$$
(A.1)
* $`p(\delta _2)=0`$.

$$I=\frac{1}{\sqrt{N}}_{\sqrt{N}(\delta _1+\delta _2)}^0𝑑ϵe^{\frac{1}{2}[p^{}(\delta _2)p^2(\delta _2)]ϵ^2+\mathrm{}}$$

Then fitting function is

$$f_2(N)=\overline{x}+\frac{C_1}{\sqrt{N}}+\frac{C_2}{N}+\mathrm{}$$
(A.2)

If $`p(\delta _2)`$ is a very small number then there might exist interval of $`N`$ where eq. (A.2) provides better fit, while in the limit of large $`N`$ eq.(A.1) is valid. The data for our observables ($`<F_{max}>`$, $`<\rho >`$, $`\chi _{Z(2)}(I,I)`$) indicates just this situation. For this reason we made fits with eq.(A.1) discarding small values of $`N_{cop}`$.


text:
# Matter-wave entanglement and teleportation by molecular dissociation and collisions

## Abstract

We propose dissociation of cold diatomic molecules as a source of atom pairs with highly correlated (entangled) positions and momenta, approximating the original quantum state introduced by Einstein, Podolsky and Rosen (EPR) \[Phys. Rev. 47, 777 (1935)\]. Wavepacket teleportation is shown to be achievable by its collision with one of the EPR correlated atoms and manipulation of the other atom in the pair.

The fundamentally profound notion of quantum teleportation is the prescription of how to map, in a one to one fashion, any quantum state of system A onto that of a distant system B: one must measure the pertinent observables of A, then manipulate their counterparts in B according to the communicated results of the measurements on A . Teleportation has thus far been explored for photon polarizations , optical field quadratures  and multi-atom spin components . Here we pose the basic question: How to teleport the translational quantum states (wavepackets) of material bodies over sizeable distances? We propose dissociation of cold diatomic molecules as a source of atom pairs with highly correlated (entangled) positions and momenta, approximating the original quantum state introduced by Einstein, Podolsky and Rosen (EPR) . Wavepacket teleportation  is shown to be achievable by its collision with one of the EPR correlated atoms and manipulation of the other atom in the pair.

Consider a cold beam of ionized molecules that are moving with high, constant velocity $`v_z`$ along the $`z`$-axis (which thus plays the role of time) to a region where they are dissociated by means of an electromagnetic pulse. Let the size $`L`$ of the dissociation region along the perpendicular $`x`$-axis be defined by means of an aperture (Fig. 1). Here the molecule splits and its two constituents start receding.

In order to be coherent within the dissociation region $`L`$, the translational center-of-mass (COM) state of the molecule along the $`x`$-axis should be close to the minimum uncertainty state (MUS) of its position and momentum, described by a Gaussian exp$`[P_x^2/2\mathrm{\Delta }P_x^2]`$, where $`P_x`$ is the $`x`$-component of the COM momentum and $`\mathrm{\Delta }P_x`$ its spread. Such a state can be prepared by translationally cooling the molecules to the ground state of a trapping potential, then ionizing and accelerating them to the required speed $`v_z`$ (prior to dissociation). The required temperature is typically $`T`$ $``$ $`\mathrm{}^2/(Mk_BD^2)`$, where $`M`$ is the molecular mass, $`k_B`$ is the Boltzmann constant and $`D`$, the COM wavepacket size, is chosen to be $`D`$ $``$ $`L`$. A size $`D`$ $``$ 300 nm would require $`T`$ $``$$`\mu `$K for H$`{}_{}{}^{+}{}_{2}{}^{}`$ and $`T`$ $``$ 0.4 $`\mu `$K for Li$`{}_{2}{}^{}{}_{}{}^{}`$. Such temperatures are achievable by Raman photoassociation of atomic pairs in Bose condensates .

We would like the two-atom translational state obtained by molecular dissociation to closely resemble the original EPR state , whose realization has not been investigated thus far. The perfect EPR state of particles 0 and 1 is described by the wavefunction $`\mathrm{\Psi }(x_0,x_1)=\delta (x_0x_1)`$: the positions and momenta of the two particles along $`x`$ are completely uncertain, yet perfectly correlated (entangled): $`x_0=x_1`$, $`p_{x0}=p_{x1}`$. The resemblance of the dissociated state to the perfect EPR state depends on the extent to which the variances of the correlated atomic variables are below the Heisenberg uncertainty limit, so that

$`\mathrm{\Delta }x_{01}\mathrm{\Delta }P_x\mathrm{},`$ (1)
$`\widehat{x}_{01}`$ $`=`$ $`\widehat{x}_0`$ $``$ $`\widehat{x}_1`$ and $`\widehat{P}_x`$ $`=`$ $`\widehat{p}_{x_0}`$ $`+`$ $`\widehat{p}_{x_1}`$ being respectively the internuclear separation and COM momentum operators. Inequality (1) is possible since the two variances do not pertain to canonically conjugate variables. The diatomic molecule prior to dissociation is describable by an internuclear wavepacket whose spatial size is $`\mathrm{\Delta }d_v`$ (typically 0.1 nm), determined by its vibrational state. Dissociation by a short laser pulse can yield an almost exact replica of the initial internuclear wavepacket on a dissociative electronic potential surface of the molecule (Fig. 1-inset). The much broader, cold COM wavepacket keeps during dissociation the molecular MUS spread (Fig. 1) $`\mathrm{\Delta }P_x`$ $``$ $`\mathrm{}/D`$. As the internuclear wavepacket becomes nearly free shortly after dissociation, it may still resemble an EPR state. Its proximity to the perfect EPR state at $`t`$ $`=`$ $`0`$, chosen to signify the completion of dissociation, can be measured by the parameter $`s`$ that is inversely proportional to the left-hand side of (1)

$`sD/\mathrm{\Delta }d.`$ (2)
Here the width $`\mathrm{\Delta }d`$ $`=`$ $`\mathrm{\Delta }x_{01}(0)`$ $`>`$ $`\mathrm{\Delta }d_v`$ is determined by the dynamical spreading of the internuclear wavepacket, caused by the gradient of the dissociative potential surface. The ratio in (2) should be sought to satisfy $`s1`$, in accordance with inequality (1), the perfect EPR state having $`s`$ $``$ $`\mathrm{}`$. Even for the realistic values $`\mathrm{\Delta }d1`$ nm $`\mathrm{\Delta }d_v`$, and $`D`$ $``$ 300 nm, this parameter is remarkably large: $`s`$ $``$ $`300`$.

There is a noteworthy analogy between the two-particle EPR state and that of entangled two-mode light from parametric downconversion (PDC) : $`s`$ is analogous to the exponential of the “squeezing parameter” of such light, which represents the ratio of the standard deviation of the “stretched” field quadrature to that of its “squeezed” counterpart (one quadrature corresponding, e.g., to the sum and the other to the difference of the two fields, analogously to our $`\widehat{x}_{01}`$ and $`\mathrm{}/\widehat{P}_x`$, respectively). In current optical experiments  $`s`$ $``$4.

Outside the dissociation region (typically, for separations beyond a few nm) the receding particles evolve freely, with diminishing position correlation. Yet the two-particle wavefunction remains inseparable, having the form $`\mathrm{\Psi }(x_0,x_1)`$ $`=`$ $`\psi (x_{01})\psi _{\mathrm{COM}}(x_0+x_1)`$, i.e., a product of the internuclear wavefunction and its COM counterpart. The spread of the coordinate difference $`x_{01}`$ then grows with time as $`\mathrm{\Delta }x_{01}^2(t)`$ $`=`$ $`\mathrm{\Delta }d^2`$ $`+`$ $`\mathrm{\Delta }v_{01}^2t^2`$, $`\mathrm{\Delta }v_{01}`$ being the spread of their velocity difference acquired during dissociation (here we assume homonuclear diatomic molecules, for simplicity). Yet the growth of $`\mathrm{\Delta }x_{01}(t)`$ can be compensated by a suitable lens, which ideally “time-reverses” $`\mathrm{\Delta }x_{01}`$, i.e., projects it back to its initial “spot” of size $`\mathrm{\Delta }x_{01}(0)`$ $`=`$ $`\mathrm{\Delta }d`$. In practice, such compensation is limited by the resolution of the focussing lens (see below).

We suggest that the translationally-entangled particles 0 and 1 can be used to demonstrate the hitherto unobserved original prediction of EPR , concerning the nonlocality of quantum correlations in a free-propagating two-atom state satisfying inequality (1): detecting the momentum or position of particle 1 would make the corresponding variable of particle 0 assume a nearly well-defined value (with fluctuations $`\mathrm{}/D`$ or $`\mathrm{\Delta }d`$, respectively). By contrast, there will be large fluctuations if we measure uncorrelated variables (Fig. 2 (a), (b)): the momentum of particle 1 and the position of particle 0, or vice versa. We note that the preparation of EPR states of internal atomic observables (unlike the present external or translational observables) by diatomic dissociation was discussed in  (for atomic excitations or pseudospin states), and in  (for atomic spin states). An EPR correlation of trapped (rather than free-propagating) atoms is proposed in .

The major challenge is to teleport the quantum state $`|\psi _{in}`$ of the transversal motion of particle 2 (along the $`x`$ axis) to particle 0 (Fig. 1). To this end, particle 2 collides with one member of the EPR pair – particle 1 (in a synchronous fashion determined by laser pulses—see below), after which both particles 1 and 2 are detected. Essentially, the collision of particles 1 and 2 allows us to project their pre-collisional joint state onto the basis of EPR-correlated states (specified below) and detect the result of the projection. The results of the post-collision detection are used to control the evolution of the “readout” particle 0, the other member of the EPR pair. In the optimal case, the resulting translational state $`|\psi _{out}`$ of particle 0 would closely imitate the input $`|\psi _{in}`$ of particle 2. We may assume that $`|\psi _{in}`$ is prepared by diffraction on a double-slit screen or grating . If measured directly, the beam of particle 2 would exhibit a characteristic diffraction pattern. If the teleportation scheme is applied instead, then the input beam is destroyed but nearly the same diffraction pattern can be observed in the transformed output beam of particle 0. We note that the collision region in Fig. 1 plays an analogous role to the beam-splitter which mixes a field quadrature of the teleported optical state with the field of one of the entangled PDC modes in the teleportation scheme of .

A natural choice for the post-collision correlated (EPR-state) basis is the set of states associated with a sharp (well-defined) sum of the colliding particles’ momenta $`\widehat{p}_{x1}`$ $`+`$ $`\widehat{p}_{x2}`$ and sharp difference of their $`x`$ coordinates, $`\widehat{x}_1`$ $``$ $`\widehat{x}_2`$. At the time instant $`\tau `$, corresponding to the particles’ nearest approach in the absence of interaction (see below), the operators $`\widehat{x}_{}`$ and $`\widehat{p}_+`$, given by

$`\widehat{x}_{}`$ $``$ $`\widehat{x}_1(\tau )\widehat{x}_2(\tau )`$ (3)
$`\widehat{p}_+`$ $``$ $`\widehat{p}_{x1}(\tau )+\widehat{p}_{x2}(\tau )`$ (4)
are measured, with values $`x_{}`$, $`p_+`$ and precision $`\mathrm{\Delta }x_{\mathrm{meas}}`$ and $`\mathrm{\Delta }p_{\mathrm{meas}}`$, respectively (which are discussed below).

The last step required for the teleportation of the state of particle 2 is the shift of the $`x`$-coordinate and momentum of the “readout” particle 0, according to the inferred values of $`x_{}`$ and $`p_+`$ as $`\widehat{x}_0`$ $``$ $`\widehat{x}_0^{}`$ $`=`$ $`\widehat{x}_0x_{}`$, and $`\widehat{p}_{x0}`$ $``$ $`\widehat{p}_{x0}^{}`$ $`=`$ $`\widehat{p}_{x0}+p_+`$. Let the precision of the position and momentum shifters be $`\mathrm{\Delta }x_{\mathrm{shift}}`$ and $`\mathrm{\Delta }p_{\mathrm{shift}}`$, respectively. From the correlations of $`\widehat{x}_0`$ and $`\widehat{x}_1`$, or $`\widehat{p}_{x0}`$ and $`\widehat{p}_{x1}`$ discussed above, we find that (at the relevant times) the shifted variables satisfy

$`\widehat{x}_0^{}`$ $`=`$ $`\widehat{x}_2\pm \mathrm{\Delta }x_T,`$ (5)
$`\widehat{p}_{x0}^{}`$ $`=`$ $`\widehat{p}_{x2}\pm \mathrm{\Delta }p_T,`$ (6)
namely, they approximately reproduce the teleported variables of particle 2, to the accuracy given by $`\mathrm{\Delta }x_T`$ and $`\mathrm{\Delta }p_T`$ where

$`\mathrm{\Delta }x_T^2=\mathrm{\Delta }d^2+\mathrm{\Delta }x_{\mathrm{meas}}^2+\mathrm{\Delta }x_{\mathrm{shift}}^2,`$ (7)
$`\mathrm{\Delta }p_T^2=\mathrm{\Delta }P_x^2+\mathrm{\Delta }p_{\mathrm{meas}}^2+\mathrm{\Delta }p_{\mathrm{shift}}^2.`$ (8)
Equations (5), (6) imply that the Wigner $`xp`$ distribution $`W(x_0^{},p_{x0}^{})`$ of the readout particle 0 is given by that of the input particle 2, $`W_{in}(x_2,p_{x2})`$, coarse-grained by a smoothing function whose width is determined by Eqs. (7), (8).

The teleportation has non-classical properties provided: (i) $`\mathrm{\Delta }x_T`$ and $`\mathrm{\Delta }p_T`$ satisfy

$`\mathrm{\Delta }x_T\mathrm{\Delta }p_T<\mathrm{};`$ (9)
(ii) they are finer than the scales of the input wavepacket variation in position and momentum, respectively \[Fig. 2(c)-(e)\]. If the fluctuations produced by the measurements and shifts are negligible ($`\mathrm{\Delta }x_{\mathrm{meas}}`$, $`\mathrm{\Delta }x_{\mathrm{shift}}`$ $``$ $`\mathrm{\Delta }d`$, and $`\mathrm{\Delta }p_{\mathrm{meas}}`$, $`\mathrm{\Delta }p_{\mathrm{shift}}`$ $``$ $`\mathrm{\Delta }P_x`$), then the left-hand side of Eq. (9) is approximately $`\mathrm{\Delta }x_T\mathrm{\Delta }p_T`$ $`\mathrm{}/s`$. For teleported Gaussian wavepackets, the maximum attainable fidelity, i.e., the overlap of the input state with the teleported output, can be shown to be $`F_{\mathrm{max}}`$ $`=`$ $`(1+\mathrm{\Delta }x_T\mathrm{\Delta }p_T/\mathrm{})^1`$. As shown in , the maximal classical teleportation fidelity is 1/2, so that larger fidelity, as in the quantum-optical teleportation experiment  where $`F=0.58`$, attests to quantum correlations. \[“Classical teleportation” is based on simultaneous (though noisy) measurements of the non-commuting quantities $`\widehat{x}_2`$ and $`\widehat{p}_{x2}`$ and preparing a corresponding MUS wavepacket of particle 0.\]

We proceed to discuss the experimentally relevant aspects of the envisaged teleportation. The dissociated pairs of fragments (particles) 0 and 1 can be emitted in all directions, but the dissociating field polarization can impose directionality . A significant fraction of the fragment pairs enter the aperture of the lens focussing them onto the region where the collision of particles 1 and 2 takes place. The colliding particles are assumed to be prepared in uncorrelated wavepackets propagating with the same classical velocity along $`z`$, $`v_{z1}`$ $`=`$ $`v_{z2}`$, and opposite classical momenta along $`y`$ (Fig. 3), $`m_1v_{y1}`$ $`=`$ $`m_2v_{y2}`$, such that $`v_{z1,2}|v_{y1,2}|\mathrm{\Delta }v_{x1,2}`$. Their focussing (and laser-pulse deflection) is such that in the absence of interaction the wavepackets would cross each other in the $`xy`$ plane, where their size would be smallest (“contractive” wavepackets ). Let us take both particles 1 and 2 to be equally charged. Charged-particle focussing (by electrostatic and magnetic lenses used in high-resolution microscopy) would render the spread (“spot size”) of the colliding particles’ positions as sharp as it was at the relevant times (thus achieving “time-reversal”): for particle 1 this time corresponds to the completion of dissociation, whereas for the teleported particle 2 it is the time of the state preparation (e.g., by a double slit). The arrival of the two particles in the $`xy`$ collision plane at $`\tau `$ must be synchronized by laser pulses, causing the fast switching-on of their deflection towards the collision region in the $`xy`$ plane, provided the values of $`z_1`$, $`z_2`$, $`v_{z_1}=v_{z_2}`$ are appropriate. Synchronization accuracy of $``$10 ps and $`v_z`$ $``$ $`10^3`$ m/s would bring the particles well within the Coulomb collision range (see below). The $`yz`$ extent of the two wavepackets must be $`\mathrm{\Delta }d_c`$ $``$ $`D`$, so that they overlap in the absence of Coulomb repulsion (Fig. 3–inset). Among the dissociated atom pairs entering the input aperture $`L`$, only those participating in the synchronized collision events are counted as pairs that contribute to successful teleportation.

The success of teleportation hinges upon our ability to discriminate, when detecting the post-collision states, between correlated EPR states with different values of the relevant parameters, namely, the parameter $`x_{}`$ and the momentum sum $`p_+`$ in Eqs. (3) and (4). We can infer the parameter $`x_{}`$ by measuring the relative deflection angle of the colliding particles. A quantitative estimation can be made for repulsive Coulomb interaction of two particles with mass $`m`$ and charge $`q`$. Their relative deflection angle is given by (see, e.g., ) $`\theta `$ $`=`$ $`\pi `$ $``$ $`2`$ $`\mathrm{tan}^1`$ $`\left[(x_1x_2)/R_{\mathrm{col}}\right]`$, where $`R_{\mathrm{col}}`$ $`=`$ $`mq^2/(4\pi ϵ_0p_y^2)`$ is the characteristic collision range. To ensure maximal sensitivity of the deflection angle to $`x_{}`$ it is useful to choose the experimental parameters such that $`\mathrm{\Delta }d_c`$ $``$ $`R_{\mathrm{col}}`$ $``$ $`D`$ $`\mathrm{max}|x_{}|`$. From the collisional analysis it follows that the precision of inferring $`x_{}`$ is mostly affected by the fluctuation of its conjugate momentum component $`\mathrm{\Delta }(p_{x2}p_{x1})`$ $``$ $`\mathrm{}/\mathrm{\Delta }d_v`$. The corresponding resolution of $`x_{}`$, obtained from that of $`\theta `$, is $`\mathrm{\Delta }x_{\mathrm{meas}}`$ $``$ $`\pi ϵ_0\mathrm{}D^2p_y/(2q^2m\mathrm{\Delta }d_v)`$, provided that $`p_y`$ $``$ $`\sqrt{mq^2/(4\pi ϵ_0D)}`$. In particular, for lithium ions with velocities $`v_y300`$ m/s, $`R_{\mathrm{col}}220`$nm, the resolution of position difference measurements is $`\mathrm{\Delta }x_{\mathrm{meas}}`$ $`15`$nm. The momentum-sum $`p_+`$ is measurable by the Doppler shift of Raman transitions induced by two counterpropagating laser fields, with precision better than 1 mm/s . The dominant contribution to its resolution, $`\mathrm{\Delta }p_{\mathrm{meas}}`$, then stems from the fluctuation $`\mathrm{\Delta }P_x`$, which is, for the given parameters, $`\mathrm{\Delta }P_x/m`$ $``$ 30 mm/s.

In the last stage of teleportation, particle 0 must be “kicked” (momentum-shifted) and position-shifted by external fields. The spread incurred by these shifts is incorporated into Eqs. (7) and (8), and constitutes one of the factors in the anticipated resolution used in the teleportation simulations of Fig. 2(c). Yet the precision of ion optics is currently so high that the dominant contributions in $`\mathrm{\Delta }x_T`$ and $`\mathrm{\Delta }p_T`$ stem from $`\mathrm{\Delta }x_{\mathrm{meas}}`$ and $`\mathrm{\Delta }p_{\mathrm{meas}}`$, respectively. We then get for the teleportation of lithium ion wavepackets, under the conditions specified above, the estimation $`\mathrm{\Delta }x_T\mathrm{\Delta }p_T/\mathrm{}`$ $`0.08`$ $`1`$. Such spreads would allow for much higher teleportation fidelities of continuous variables than those presently achievable in quantum optical experiments.

The foregoing analysis shows that, despite the challenging nature of the proposed ideas, their experimental implementation does not pose insurmountable (rather than technical) difficulties. The teleportation of matter wavepackets by molecular dissociation and collisions is a novel concept, combining elements of molecular dynamics and ion (atom) optics, and having quantum optical analogues. Its realization is a challenging but viable goal to pursue, en route to quantum information exchange between complex material objects.

Acknowledgments: The support of US-Israel BSF and DFG is acknowledged. G. K. holds the G. W. Dunne Professorial Chair.


text:
# 1 𝑇_𝐾 as a function of the excited states kept in the scaling procedure. The calculations were done for assuming a heavy particle of mass 𝑀∼50×𝑚ₚᵣₒₜₒₙ, a barrier of height ∼300𝐾 and of width 0.5Å, and the width of the potential wells was taken to be 0.1Å.

Few comments on preprint cond-mat/0007430 by Aleiner et al.: It has been suggested in Ref.  that conduction electron generated virtual hopping to the excited states of a two-level system (TLS) could increase the amplitude of assisted tunneling and thus result in a measurable Kondo temperature. In a recent preprint, Aleiner et al. revisited the role of excited states and pointed out that in the calculations of Ref.  only the first few excited states have been taken into account and inclusion of all excited states induces a large decrease of $`T_K`$ instead of an increase.

While we think that the work of Aleiner et al. addresses an important issue, there are several points that we have to comment on:

(a) The reduction of $`T_K`$ heavily relies on the assumption of electron-hole symmetry. In fact, using the model of Ref.  and assuming a contact TLS-electron interaction of strength $`U`$, the assisted tunneling term generated by the reduction of the cutoff is proportional to:

$$\delta V_{ll^{}}^{\alpha \beta }=\underset{ϵ_j<D}{}\underset{l^{\prime \prime }}{}\varrho (D)V_{ll^{\prime \prime }}^{\alpha j}V_{l^{\prime \prime }l^{}}^{j\beta }\varrho (D)V_{ll^{\prime \prime }}^{j\beta }V_{l^{\prime \prime }l^{}}^{\alpha j},$$
(1)
where $`D`$ denotes the bandwidth cutoff, $`ϵ_j`$ is the energy of the j’th excited state of the TLS, and the matrix element between the TLS wave functions $`\phi _i`$, $`\phi _j`$ in the ”tower” model  and conduction electron states $`j_l(z)`$ $`j_l^{}(z)`$ is given by:

$$V_{ll^{}}^{ij}U𝑑z\phi _i(z)\phi _j(z)j_l(z)j_l^{}(z).$$
(2)
For large enough $`D`$’s the identity $`\phi _j(z)\phi _j(z^{})\delta (zz^{})`$ can be used and assuming $`\varrho (D)=\varrho (D)`$ the two terms in Eq. (1) cancel. However, for electronic bands usually $`\varrho (D)\varrho (D)`$ and the cancellation does not occur. In the Figure we show $`T_K`$ obtained from the complete solution of the leading logarithmic scaling equations. We assumed a density of states of the form

$$\varrho (ϵ)=\varrho (0)\left(1+\alpha ϵ/D_0\right),$$

where $`\alpha `$ is usually of the order of 1. As one can see, a small electron-hole asymmetry can increase $`T_K`$ by many orders of magnitude, bringing it into the physically measurable range. Note that while for $`\alpha =0`$ the inclusion of all excited states reduces $`T_K`$ by three orders of magnitude in agreement with Aleiner et al., for $`\alpha =0.30.4`$ there is very little reduction compared to the truncated model. \[The $`\alpha =0.4`$ case is special: Here five heavy particle states are still active at $`T_K`$.\]

(b) It is obvious from the figure that it is extremely difficult to give a ’first principles’ estimation for $`T_K`$, which is very sensitive to the exact model parameters such as the distance between the two wells, the mass of a heavy particle, or the value and structure of the heavy particle-conduction electron interaction. Within the Thomas-Fermi interaction approximation, e.g., $`T_K`$ changes many orders of magnitudes if one changes the effective charge of the heavy particle by a factor of 2 (Hall effect measurements indicate, e.g., that the valence of crystalline $`Cu`$ is about $`1.5`$. This value may considerably differ from the valence of the tunneling atoms that are most probably sitting in a more disordered region.)

In reality, one does not even know what tunneling centers look like. The simplistic model of Ref.  assumes a very specific TLS structure, which may be very far from that of real tunneling systems, possibly due to dislocation jogs, e.g..

(c) Our calculations clearly demonstrate that $`T_K`$ can be in the measurable range. It is a more delicate question whether the splitting of the TLS can be small enough to allow for the development of a two-channel Kondo behavior. In our model the bare splitting is $`\mathrm{\Delta }_00.4K`$, which is clearly smaller than $`T_K`$ for $`\alpha >0.2`$. However, the presence of electron-hole symmetry breaking increases $`\mathrm{\Delta }_0`$, which is, on the other hand considerably reduced under scaling. To decide which of these processes wins and whether the renormalized splitting can be less than $`T_K`$, a much more detailed study of the next leading logarithmic equations is needed.

Another possibility could be to have tunneling centers, where the degeneracy of the lowest-lying states is guaranteed by some local symmetry (rotational or possibly time reversal).

A. Zawadowski1 and G. Zaránd,1,2
1Institute of Physics, Technical University of Budapest, 2Lyman Physics Laboratory, Harvard University