text:
# Canonical and non-canonical equilibrium distribution

## I Introduction

A growing interest in non-extensive thermodynamics can be traced back to the 1988 pioneering work of Tsallis . In particular, non-extensive thermodynamics is the proper theoretical context for developing the foundations of non-canonical distribution functions and non-canonical equilibria in physical systems. On the basis of the subsequent development of the original work by Tsallis, it has been argued that canonical equilibria, the distributions that form the basis of equilibrium statistical mechanics, are not generic. Rather, canonical equilibria are singular in a more general form of equilibrium, called a generalized canonical equilibrium . A number of recent papers elaborate on this idea , but we find  to be of special interest addressing as it does the foundations of both canonical and non-canonical equilibria. Rajagopal and Abe  prove that the equilibrium described by the canonical distribution is not uniquely determined by the micro-canonical distribution, as one finds in text books. In fact if the phase space has a fractal, rather than smooth, structure, a non-canonical distribution will result. However, although canonical equilibria are not unique, and non-canonical equilibria are possible, one might conclude on the basis of the arguments in  that the form of the non-canonical distribution is uniquely of the form established by Tsallis and his co-workers . As attractive as the probabilistic and entropic arguments of Ref.  are, they are not indisputable, and in fact we find herein, using dynamical arguments, that this is not the case.

We show that the adoption of a dynamical approach to thermodynamical equilibrium yields a different conclusion. First of all, we argue that the non-extensive condition based on memory, and probably that resting on long-range correlations as well, has the striking effect of making the role of dynamics much more important than in the case of ordinary statistics. The dynamic approach, in these non-extensive conditions, generates a form on non-canonical equilibrium that, however, departs from the generalized canonical form prescribed by non-extensive thermodynamics. To make it easier for the reader to follow our arguments and to understand the purpose of the paper, illustrated in Section I D, we shall first discuss the different natures of the entropic, dynamic and stochastic approaches to equilibrium.

### A Non-extensive entropic indicator

The argument for the non-extensive thermodynamics of Tsallis goes as follows. First of all, the conventional entropy of Gibbs

$$S(\mathrm{\Pi })𝑑x\mathrm{\Pi }(x)ln\mathrm{\Pi }(x),$$
(1)
is replaced by the non-extensive entropic indicator

$$S_q(\mathrm{\Pi })𝑑x\frac{1\mathrm{\Pi }(x)^q}{q1}.$$
(2)
Secondly, we have to apply a method of entropy maximization under given physical constraints to determine the most plausible shape of the unknown probability density function $`\mathrm{\Pi }(x)`$ . The first constraint is on the normalization of the distribution $`\mathrm{\Pi }(x)`$:

$$𝑑x\mathrm{\Pi }(x)=1.$$
(3)
The second constraint is on the first moment of the variable $`x`$ itself. According to the most recent prescription of Ref. the constraints on $`x`$ must be applied on the mean value $`U_q`$ defined by

$$U_q\frac{𝑑xx\mathrm{\Pi }(x)^q}{𝑑x\mathrm{\Pi }(x)^q}.$$
(4)
It has to be pointed out that in Ref. the physical meaning of $`x`$ is that of energy. Here, we shall interpret $`x`$ as the “coordinate” of an overdamped particle driven by a fluctuation-dissipation process resulting from the interaction with a non-conventional “thermal bath”. When the friction term can be neglected we are expected to recover the results of the earlier work of Ref.. Thus, for the same reasons  as those illustrated in Ref. we set a constraint on the first moment of $`|x|`$. The result of entropy maximization subject to the imposed constraints of Eq. (4) yields

$$\mathrm{\Pi }(x)=\left[1\frac{(1q)\stackrel{~}{b}(xU_q)}{𝑑x\mathrm{\Pi }(x)^q}\right]^{1/(1q)}/Z_q,$$
(5)
where $`Z_q`$ is a normalization factor and $`\beta `$ is a constant value stemming from the Lagrange multiplier associated with the constraint on the variable $`x`$. Note that the expression of $`\mathrm{\Pi }(x)`$ provided by the authors of Ref. is not explicit. In fact, as shown by Eq.(5), it is a functional of $`\mathrm{\Pi }(x)`$. Thus, in practice the explicit form of $`\mathrm{\Pi }(x)`$ has to be established by means of an iteration procedure. It is remarkable, however, that the micro-canonical derivation from Ref. results in the same prescription as that of Eq.(5). From the point of view of the issues under discussion here what matters is the fact that at the end of the iterative procedure $`U_q`$ becomes a well-defined number. Consequently, the resulting expression for $`\mathrm{\Pi }(x)`$ is a simple analytical formula that in the asymptotic limit, $`|x|\mathrm{}`$, with the constraint on the first moment of $`|x|`$, has the same structure as that derived in Ref.,

$$\mathrm{\Pi }(x)=\frac{\stackrel{~}{b}(2q)}{[1+\stackrel{~}{b}(q1)|x|]^{1/(q1)}}.$$
(6)
Note that the adoption of a constraint on the second moment would lead to

$$\mathrm{\Pi }(x)=\frac{\left[\frac{\stackrel{~}{b}(q1)}{\pi }\right]\frac{\mathrm{\Gamma }(1/(q1))}{\mathrm{\Gamma }((3q)/[2(q1)])}}{[1+\stackrel{~}{b}(q1)x^2]^{1/(q1)}}.$$
(7)
Note that both Eq.(6) and Eq.(7) are the result of an earlier prescription, missing the normalization factor present in Eq.(4). Nevertheless, it has to be pointed out that both Eq.(6) and Eq.(5) share the characteristic of having long tails with Lévy statistics. In fact, it is well known  that the anomalous diffusion processes of the Lévy kind, in the one-dimensional case, are characterized by probability distributions $`p(x,t)`$ whose Fourier transform in the symmetric case reads

$$\widehat{p}(k,t)=exp(b|k|^\alpha t),$$
(8)
where $`\alpha `$ is the Lévy index ranging, in principle, in the interval $`0<\alpha 2`$, and $`b`$ denotes the diffusion intensity. The inverse Fourier transform of $`\widehat{p}(k,t)`$ of Eq.(8) is characterized by the tail

$$\underset{|x|\mathrm{}}{lim}p(x,t)\frac{t}{|x|^{1+\alpha }},$$
(9)
which would lead immediately to the anomalous entropy index

$$q=1+1/(1+\alpha ).$$
(10)
We have to stress, however, that except for the case $`\alpha =1`$, corresponding to the ballistic limit, the central part of the distribution produced by the Lévy process of Eq.(8) cannot be expressed in an analytical form and can significantly depart from the analytical form of both Eqs.(6) and (5). The adoption of the prescription of Eq.(7) is expected to yield better agreement, but, as we shall see subsequently, it does not seem to be compatible with the nature of the dynamical approach to equilibrium. In Fig.1 we compare the Lévy distribution to both the prediction of Eq.(6) and Eq.(7) and find that there exists a significant disagreement between Tsallis and Lévy statistics, even though some authors refer to them as equivalent. In this paper we focus our attention on the origin of the difference between the two kinds of statistics. Our dynamic approach leads to a form of equilibrium that is stable, in the sense of the Lévy-Gnedenko theorem , while the generalized canonical equilibrium of Eq.(5) is not. This is evident in the free case, due to the difference between the Tsallis and the Lévy structure. In Section III we prove that this is so also in the presence of a feedback with the system of interest of the generator of fluctuations.

### B From dynamics to thermodynamics

We are convinced that there are no incontrovertible reasons why the canonical distribution should be the unique form of thermodynamic equilibrium, and on this issue we essentially agree with the point of view of Rajagopal and Abe. However, we are equally convinced that the generalized canonical form of Eq.(5) does not satisfy the stable conditions necessary for Lévy statistics, and the results of the present paper can be thought of as providing plausible evidence of that. To substantiate this view with dynamical arguments, it is convenient to concisely review the results of an earlier work.

The ambitious purpose of this earlier work was that of reversing the path from thermodynamics to mechanics established by Boltzmann. The main idea behind Ref. is as follows. The Fokker-Planck equation is a well-known description of the evolution of the probability density in phase space. Many attempts have been made in the literature to derive this equation. However, all these attempts rest on the assumption that the bath, responsible for the Brownian behavior of the particle of interest, is given by a set of harmonic oscillators. This means that a statistical assumption must be made on the initial condition of the bath, which is arbitrarily given a canonical equilibrium distribution, corresponding to a given temperature $`T`$ of the thermal bath.

Consequently, this kind of approach to Brownian motion is only partially dynamical, since significant use is already made of statistical mechanics, and thus of thermodynamics . The authors of Ref. adopted a totally different approach. They assumed that a given oscillator, playing the role of a stochastic Brownian particle, interacts with another Hamiltonian system, which should play the role of a bath. They called this second system a *booster*, to stress that the ensuing approach has to rest only on the dynamical properties of this kind of bath, with no use whatsoever of thermodynamical arguments. After establishing the dynamical conditions ensuring the validity of the Fokker-Planck equation, it is found that the oscillator of interest reaches a canonical equilibrium distribution. Due to the nature of the procedure adopted the width of this canonical distribution depends only on the parameters of the Hamiltonian system under study. Consequently, it is possible to derive a mechanical expression for temperature. This is the key result of Ref., which reads

$$k_BT=\left[\frac{}{E}lnA(E)+\frac{}{E}ln\{\xi ^2_{eq}Re[\widehat{\mathrm{\Phi }}_\xi (\omega )]\}\right]^1.$$
(11)
Note that $`\xi `$ denotes the *doorway* variable, namely the variable of the booster through which the interaction between the booster and the oscillator of interest is established. The symbol $`\widehat{\mathrm{\Phi }}_\xi (\omega )`$ denotes the Laplace transform of the correlation function of $`\xi `$ evaluated at the oscillation frequency of the oscillator. The structure of this expression reflects the application of linear response theory. The correlation function whose Laplace transform is in Eq. (11) is evaluated assuming the booster to be in a micro-canonical equilibrium with energy $`E`$ and this condition is not affected by the interaction with the oscillator. The symbol $`A(E)`$ denotes the number of states of the booster in the same physical condition, and, consequently, obeys the ordinary prescription

$$A(E)E^{N/2},$$
(12)
where $`N`$ is the number of degrees of freedom in the booster.

The authors of Ref. note that for $`N\mathrm{}`$,

$$\frac{}{E}ln\{\xi ^2_{eq}Re[\mathrm{\Phi }_\xi (\omega )]\}<<\frac{}{E}lnA(E).$$
(13)
This means that in the limiting case of infinitely many degrees of freedom the predictions of Boltzmann are recovered: $`1/(k_BT)=\frac{}{E}lnA(E)`$. Note that the condition of Eq.(13) holds true if the doorway variable $`\xi `$ depends on a number of particles which is kept fixed with $`N\mathrm{}`$. In other words, the condition of Eq.(13), or equivalently, Boltzmann’s principle, rests on an interaction condition, which is crucial for the extensive statistical mechanics perspective to apply.

Let us see this aspect in detail. In Ref. the booster is the the well known Fermi-Pasta-Ulam (FPU) system. This is a chain of particles interacting with one another via nonlinear interactions. Let us consider two opposite conditions. In the former, which is the one adopted in Ref., the oscillator of interest interacts only with the first particle of the FPU chain. In the latter the oscillator of interest interacts with all the particles of the FPU chain, with interaction strength of comparable intensities. The former condition refers to a short-range interaction, confined to the positions of the oscillator of interest and to the first particle of the FPU chain. In the latter condition, the interaction extends over the whole FPU chain. It is evident that the former condition fits the inequality of Eq.(13), whereas the latter does not. This means that the latter condition results in a non-extensive form of dynamics, with a consequent breakdown of the prescriptions of ordinary statistical mechanics. The former condition, on the contrary, for $`N\mathrm{}`$, recovers ordinary statistical mechanics. This suggests that the dynamical corrections to the Boltzmann principle, recorded in Ref. for relatively small values of $`N`$, are a manifestation of incipient non-extensive statistical mechanics, and so are very close to the breakdown of the Fokker-Planck treatment on which the analysis Ref. rests. Note that when we assign to the oscillator of interest a very low frequency, compared to the booster frequencies, the quantity $`Re[\widehat{\mathrm{\Phi }}_\xi (\omega )]`$ turns out to virtually coincide with the time scale $`T_{micro}`$ of the variable $`\xi `$, defined by

$$T_{micro}_0^{\mathrm{}}\mathrm{\Phi }_\xi (t)𝑑t<\mathrm{}.$$
(14)
Herein we focus our attention on the case where

$$lim_t\mathrm{}\mathrm{\Phi }_\xi (t)=const/t^\beta ,$$
(15)
with

$$0<\beta <1.$$
(16)
The inverse power law form of the correlation function means that we select time memory as the source of violation of extensivity, rather than long-range spatial interactions. The parameter $`T_{micro}`$ denotes the correlation time of the fluctuating variable $`\xi `$. Here we consider the dichotomous case, where the variable $`\xi `$ has two distinct values, $`W`$, and $`W`$, with fluctuating time durations. The correlation function $`\mathrm{\Phi }_\xi (t)`$ is proven  to be proportional to the second time derivative of the distribution of waiting times in the two states of the variable $`\xi `$. This function, as we shall see in Section I C, depends on another time, $`\overline{T}`$, as well as on $`\beta `$.

In conclusion, the work of Ref. establishes the dynamical conditions necessary to derive canonical equilibrium. The authors prove that canonical equilibrium implies an interaction with a booster with a finite time scale and a short-range interaction. We can also observe that in this case the use of dynamical arguments yields the same conclusions as the very simple argument based on the law of large numbers. The advocates of the law of large numbers for the foundation of statistical mechanics might judge the dynamical perspective to be of limited use. We prefer to interpret the conclusion of Ref. as proof that in the extensive case the canonical distribution description of thermodynamic equilibrium is unique. The numerical work of Ref. shows that it is very difficult to detect the dynamical corrections to the Boltzmann principle illustrated by Eq.(11) from within an extensive perspective. This is so because the dynamical corrections become significant when the booster is small, and, as clearly pointed out by the recent work of Gross, a small system is essentially non-extensive. Thus, the numerical findings of Ref. refer to a condition very close to the breakdown of the extensive condition and consequently of the canonical equilibrium on which Eq.(11) rests. This means that the dynamical approach can be of great utility. In the case of non-extensive statistical mechanics, the dynamic approach is probably the only non-ambiguous way to address the subtle issues posed by the entropic and probabilistic methods.

### C Stochastic dynamics

In this section we explain the nature of our dynamical approach to the non-canonical equilibrium distribution. This approach rests on a stochastic method adapted to the need of realizing a dichotomous variable $`\xi `$ with the two possible values, $`W`$ and $`W`$, and with the waiting time distribution

$$\psi (t)=\frac{(\beta \overline{T})^{\beta +1}(\beta +1)}{(\beta \overline{T}+t)^{2+\beta }},$$
(17)
where $`\beta `$ ranges in the interval of Eq.(16). Note that the parameter $`\overline{T}`$ is of crucial importance to define the time scale of our process and corresponds to the mean residence time in either of the two states of the velocity variable $`\xi `$. In accordance with the prescriptions of Ref. this waiting-time distribution yields the kind of correlation function $`\mathrm{\Phi }_\xi (t)`$ that we plan to study herein (see Eq.(15)). In fact, as shown in Ref., the form of Eq.(17) yields the correlation function

$$\mathrm{\Phi }_\xi (t)=\frac{(\beta \overline{T})^\beta }{(\beta \overline{T}+t)^\beta },$$
(18)
fitting the asymptotic time limit of Eq.(15). This condition can be considered as the natural one-dimensional counterpart of the two-dimensional billiards of Zaslavsky . From this point of view, our commitment to the adoption of a merely dynamical approach is not broken, since the stochastic approach that will be adopted in Sections II and III is statistically equivalent to the adoption of the intermittent map of Ref., on which the theoretical work of Ref. is based.

We stress that the equivalence between a dynamical map and a stochastic process is the reason why contact can be established between the dynamical and the entropic approaches. In fact, as shown in Ref. the non-extensive Tsallis entropic indicator serves the purpose of guessing the most convenient form for the transition probability $`\mathrm{\Pi }(x)`$, which is then related to the waiting-time distribution $`\psi (t)`$ through the basic property

$$\mathrm{\Pi }(x)=\frac{1}{2}\psi (x/W)/W.$$
(19)
The factor of $`1/2`$ takes into account that the probability of making the jump $`x`$ is equal to that of making the jump $`x`$, namely, the jump in the opposite direction. The authors of Ref. show that the left-hand term of this equality, Eq.(19), can be predicted using entropic arguments, while the right-hand term of the same equation is dictated by dynamical arguments based on the intermittent map of Ref.. These dynamical arguments are supplemented by the assumption of random injection of the trajectory from the chaotic into the laminar region of the intermittent map, an argument leading to an analytical prediction for $`\psi (t)`$ in complete agreement with the numerical observation of dynamics.

The earlier work of Ref. established that the adoption of the method of entropy maximization applied to the non-extensive entropy of Eq.(2) results in a form of $`\mathrm{\Pi }(x)`$ which is compatible with the birth of Lévy statistics. However, the $`\mathrm{\Pi }(x)`$ thus derived is not the equilibrium distribution of the variable $`x`$. Rather it is the probability for the random walker to make a jump of length $`|x|`$. This is a stationary property determined by the special kind of booster here under study. The ensuing diffusion process yields a Lévy form as a result of the Lévy-Gnedenko theorem . As shown in Fig. 1, the shape of this distribution departs from the form of the generalized canonical distribution of Eq.(5), even if we adopt the constraint on the second moment, in spite of the fact that it does not fit the nature of the dynamical approach here illustrated.

To account for this discrepancy we might make the conjecture that the comparison between the Lévy statistics and the Tsallis generalized canonical distribution is not appropriate. The former refers to a diffusion process, and the latter to an allegedly equilibrium condition. Actually, after exploring this possibility we shall conclude that the latter, at least in the case of the dynamic model of the present paper, cannot reflect an equilibrium property. However, at the present stage, we are forced to develop a picture comparable to that of the earlier work of Ref.. We have to study a case where the variable $`x`$ not only undergoes the influence of the diffusion producing fluctuations, but it produces a feedback on its own “bath”, balancing the diffusion process, so as to create an equilibrium condition. This is the condition to compare to the generalized canonical distribution of Eq.(5). In other words, to address, from a dynamical perspective, the issues recently dealt with by Rajagopal and Abe we cannot disregard the feedback from the system to the “bath.”

### D Purpose of the present paper

At this stage it is much easier for us to illustrate the main purpose of the present paper. We want to explore a condition where the booster does not fulfill the key condition of Eq.(14), that is, the microscopic relaxation time diverges, and we want to assess whether or not this divergence leads to the generalized canonical distribution of Eq.(5). We aim at answering the question: Which is the form of the equilibrium distribution reached by the system of interest when the time scale of the *booster* is infinite?

To answer this question we do not adopt, as done in Ref. , a Hamiltonian approach. The latter approach is difficult for obvious reasons. A numerical simulation check of the theoretical prediction would imply technical difficulties caused by the slowness of the booster itself. This means that the anomalous booster is replaced, as done in earlier papers, by a generator of a dichotomous fluctuation with an inverse power law distribution of waiting times. The essential ingredient of the approach of Ref., not yet present in the dynamical derivation of a free Lévy diffusion, is a feedback of the diffusing variable on the generator of the fluctuations. This aspect has already been considered in the earlier work of Ref.. However, in that paper the dynamical approach to the feedback was replaced by a phenomenological friction, which did not allow the authors of that paper to keep the promise of resting solely on dynamical arguments at any level. In conclusion, we adopt the program of Ref., based on the observation of the fluctuation-dissipation process caused by the interaction of a particle with a booster having no finite time scale. The presence of feedback serves the purpose of balancing the diffusion process with dissipation so as to result eventually in an equilibrium condition.

With the program of Ref. in mind, we have to refer ourselves to the correlation function, Eq.(18), with the index $`\beta `$ fulfilling the condition of Eq.(16), and so implying the breakdown of the time scale separation between the macroscopic and the microscopic levels. The realization of the program of statistical mechanics requires an accurate definition of the process of memory erasure associated with the transition from one to the other velocity state. This is more conveniently defined by the waiting-time distribution $`\psi (t)`$ of Eq.(17), than by the correlation function of Eq.(18). We see that even if the condition of Eq.(16) applies, the time $`\overline{T}=_0^{\mathrm{}}t\psi (t)𝑑t`$ remains finite. Memory of microscopic dynamics is lost in times $`t>>\overline{T}`$.

At this stage it is convenient to support our claim concerning memory erasure with arguments borrowed from the earlier investigation of Gaspard and Wang. These authors studied the Kolmogorov complexity of the Manneville map and found that in the regime corresponding to the dynamical foundation of Lévy processes the Kolmogorov complexity is a linear function of time. This means that at a given time $`t`$ the number of transitions $`M`$ from the one to the other laminar region is given by $`Mt/\overline{T}`$. This means that for $`t\mathrm{}`$ the conditions for the realization of the generalized version of the central limit theorem are fulfilled, since the function defined by

$$p(x,M)\mathrm{\Pi }(x)\mathrm{\Pi }(x)\mathrm{}\mathrm{\Pi }(x)p_{in}(x),$$
(20)
where the asterix denotes a convolution, and for $`M\mathrm{}`$, Eq. (20) tends to the Lévy distribution. One might be tempted to make the conjecture that in the presence of feedback the number $`M`$ cannot increase beyond some limit, and that equilibrium is reached with a relatively low value of $`M`$ so as to allow the distribution to maintain the structure of a generalized canonical distribution *á la* Tsallis. In this paper we limit our analysis to the case where the linear response theory of Section III holds true. This forces us to adopt a feedback so weak that this possibility does not emerge from our simulations.

## II Dynamical model

The dynamical model studied here depends on a process of free diffusion with feedback, established through control of a dynamical parameter of the booster. This feedback process has the role of balancing diffusion so as to realize dynamically the equilibrium condition. We illustrate first the dynamical model used as the generator of free diffusion, then we show how the feedback is realized, and finally we illustrate the numerical technique adopted.

### A Free diffusion

Formally free diffusion is realized using the equation of motion

$$\frac{dx(t)}{dt}=\xi (t).$$
(21)
Here the dynamical variable $`x(t)`$ denotes either a spatial coordinate or a velocity. In the former case the variable $`\xi (t)`$ has to be considered a fluctuating velocity, while in the latter case it is regarded as a fluctuating acceleration. The results are equivalent, and the reader can adopt either of them, even if we consider $`x(t)`$ to be a spatial coordinate so as to make the connection with earlier work  more natural.

We assume the variable $`\xi (t)`$ to be dichotomous, namely, we shall assign to this variable only two distinct values, either $`W`$ or $`W`$. The motivation for this choice is not only simplicity, but also has to do with the main purpose of this paper, that of assessing the specific form the non-canonical equilibrium. Also, of course, and prior to this purpose, that of proving that deviations from canonical equilibrium can be generated as a consequence of the breakdown of the ordinary condition of time-scale separation between macroscopic and microscopic processes. In the limiting case where the variable $`\xi (t)`$ is Gaussian, and consequently produces a Gaussian equilibrium distribution , the resulting diffusion process is Gaussian and the final equilibrium distribution can turn out to be Gaussian as well , as a result of having established a seed of ordinary statistics at the microscopic level. For this reason the choice of microscopic statistics, departing from the ordinary statistical condition is crucial, and the adoption of the dichotomous assumption serves the purpose of eventually establishing a non-canonical equilibrium.

The variable $`\xi (t)`$ keeps one of the two possible values for times with a random duration. Thus, a statistical treatment must be adopted and a waiting-time distribution $`\psi (t)`$ is used. The time intervals of sojourn in a given state are labelled by an integer index $`k`$ running from $`1`$ to $`+\mathrm{}`$. At the end of any sojourn the variable can either change its value, and thus make a transition from $`W`$ to $`W`$ or from $`W`$ to $`W`$, or it can keep its original value. The probability of changing values and that of keeping the same value are equal, and thus equal to $`1/2`$. We can interpret the variable $`\xi (t)`$ as a function of the continuous time $`t`$ by adopting the following prescription:

$$\xi (t)=\underset{k}{}\xi _k\chi _{t_k,t_{k+1}}(t),$$
(22)
where the times $`t_k`$ and $`t_{k+1}`$ denote the beginning and the ending time of the $`kth`$ sojourn, respectively, and the values $`\xi _k`$ are randomly assigned either the value $`W`$ or the value $`W`$ with equal probability and $`\chi _{t_k,t_{k+1}}(t)=1`$ if $`t_kt<t_{k+1}`$. The sequel of the times $`t_k^{}`$ is fixed by the waiting time density distribution $`\psi `$, by selecting a given value $`\tau _k`$, of the interval $`[\tau _k,\tau _k+ϵ]`$, with $`ϵ<<1`$ with the probability $`\psi (\tau _k)ϵ`$. The waiting time distribution $`\psi `$ is assigned an inverse power law form determined by the constraint of yielding the correlation function

$$\mathrm{\Phi }_\xi (t)=\frac{(\beta \overline{T})^\beta }{(\beta \overline{T}+t)^\beta },$$
(23)
where $`\beta `$ is a positive number determining the integrability of the correlation function. If $`\beta >1`$ the correlation function, Eq.(23), is integrable and the microscopic time scale

$$T_{micro}_0^{\mathrm{}}\mathrm{\Phi }_\xi (t)𝑑t=\frac{\beta \overline{T}}{\beta 1}$$
(24)
can be defined. According to Refs. the form of $`\psi (t)`$ is determined by that of $`\mathrm{\Phi }_\xi (t)`$ through the relation

$$\psi (t)=\overline{T}\frac{d^2}{dt^2}\mathrm{\Phi }_\xi (t),$$
(25)
which yields

$$\psi (t)=\frac{(\beta +1)(\beta \overline{T})^{\beta +1}}{(\beta \overline{T}+t)^{\beta +2}}.$$
(26)
We note that the parameter $`\overline{T}`$ appearing in both Eq.(23) and Eq.(26) is the mean sojourn time. Thus, in some sense, when the microscopic time of Eq.(24) becomes infinite, the role of the microscopic time scale is played by $`\overline{T}`$.

By integration of Eq.(21) we get

$$x(t)=\xi _n(tt_{n1})+\underset{k=0}{\overset{n1}{}}\xi _k\tau _k,,$$
(27)
where $`t`$ is a time located in the interval $`[t_{n1},t_{n1}+\tau _n]`$ with $`t_{n1}=_{k=0}^{n1}\tau _k`$, in accordance with the earlier prescriptions. Furthermore, the transition probability of Section I A is obviously related to the waiting time distribution by Eq.(19). Note that, as discussed in detail in Refs., there are two distinct basins of attraction for the diffusion process resulting from the repeated occurrence of the transition process of Eq.(19), the Gauss basin for $`\beta >1`$ and the Lévy basin for $`\beta <1`$.

### B The booster linear response and the diffusion control

To control diffusion in such a way so as to generate a known equilibrium distribution we adopt a model based on replacing Eq.(21) with

$$\frac{dx(t)}{dt}=\xi _x(t),$$
(28)
which means that the motion of the fluctuating variable $`\xi _x(t)`$ also depends on the value of the dynamic variable $`x`$. We make the simplifying assumption that the *relaxation time*, or response time, of the booster, $`T_B`$, is shorter than that of the particle of interest, $`T_R`$, so as to ensure for the booster a *time dependent equilibrium condition* determined by the state of the system of interest. In other words, the booster is assumed to be in a condition of equilibrium determined by the variable $`x`$. As the variable $`x`$ moves from an initial condition $`x(0)`$ out of equilibrium state, i.e. from a position larger than that fluctuations of $`x`$ itself, to the final equilibrium, the booster correspondingly moves from an equilibrium to another equilibrium condition. All this is illustrated by the following mathematical arguments. These arguments do not aim at providing a rigorous mathematical treatment, but rather a heuristic treatment reflecting the physical condition that we are assigning to the booster.

To make the analysis in terms of a strictly stochastic process, we build up the following regularized process $`\stackrel{~}{x}_{T_\rho }(t)`$, defined by:

$`{\displaystyle \frac{d\stackrel{~}{x}_{T_\rho }(t)}{dt}}`$ $``$ $`{\displaystyle \frac{x(t)x(tT_\rho )}{T_\rho }}={\displaystyle \frac{1}{T_\rho }}{\displaystyle _{tT_\rho }^t}\xi _x(t^{})𝑑t^{}`$ (29)
$``$ $`{\displaystyle \frac{1}{T_\rho }}{\displaystyle _{tT_\rho }^t}\xi _{\stackrel{~}{x}}(t^{})𝑑t^{},`$ (30)
where $`T_\rho `$ is a regularization time and in the last equivalence it has been explicited the spatial dependence of the noise on $`\stackrel{~}{x}`$. This new process $`\stackrel{~}{x}`$ has the characteristics of $`x(t)`$ in the mean, but its fluctuations are smoothed by using the temporal mean.

If the response time $`T_B`$ of the bath is fast and if $`T_\rho `$ is long enough, but always smaller than relaxation time $`T_R`$ of the variable $`x`$, one can suppose that the phase space mean is almost equivalent to the temporal mean:

$$\frac{d\stackrel{~}{x}_{T_\rho }(t)}{dt}=\frac{1}{T_\rho }_{tT_\rho }^t\xi _{\stackrel{~}{x}}(t^{})𝑑t^{}\xi _{\stackrel{~}{x}}.$$
(31)
Now, assuming linear response is valid, the mean value of the thermal bath variable depends on $`\stackrel{~}{x}`$ as:

$$\xi _{\stackrel{~}{x}}=\gamma \stackrel{~}{x}.$$
(32)
With the linear response hypothesis inserted into Eq. (31) the equation for the regularized process

$$\frac{d\stackrel{~}{x}_{T_\rho }(t)}{dt}\gamma \stackrel{~}{x}_{T_\rho }(t)$$
(33)
is obtained. The solution $`\stackrel{~}{x}_{T_\rho }(t)=x(0)e^{\gamma t}`$ shows a relaxation toward the equilibrium $`x0`$, on a time scale $`T_R1/\gamma `$.

The analogy with a stochastic equation, such as the Ornstein-Uhlenbeck or its extension to Lévy noise , is completed by adding a fictitious stochastic noise $`d\stackrel{~}{\xi }_{T_\rho }(t)`$ in (33), that summarizes the fluctuations on the time scale $`T_\rho `$ and that compensates for the approximation inserted into Eq. (31):

$$d\stackrel{~}{x}_{T_\rho }(t)=\gamma \stackrel{~}{x}_{T_\rho }(t)dt+d\stackrel{~}{\xi }_{T_\rho }(t)$$
(34)
Hence, with the assumptions of linear response theory and relaxation times of the thermal bath faster than of that of the variable of interest (partial equilibrium), a diffusive process becomes an equilibrium process. We emphasize that the analogy (34) is true only if:

$$T_B<<T_\rho <<T_R=\gamma ^1.$$
(35)
By using this simplified analysis it is also possible to make qualitative predictions about the equilibrium distribution.

When $`x0`$, the feedback $`\gamma \stackrel{~}{x}_{T_\rho }(t)`$ can be neglected with respect to the fluctuations $`\stackrel{~}{\xi }_{T_\rho }(t)`$. This limit corresponds to free diffusion for which the central limit theorem (generalized or not) can be used. The distribution is Gaussian or Lévy’s according to whether the correlation function is integrable or not.

When $`|x|>>0`$ the feedback term (32) prevails over fluctuations in Eq.(34). Then the dynamical variable is pushed back to values near the equilibrium $`x0`$ as an attractive field and, so, quenching the diffusion. This description corresponds, at least for the region $`|x|<<W/\gamma `$, to what one could obtain if in Eq. (34) the fictitious noise $`\stackrel{~}{\xi }_{T_\rho }(t)`$ is simply stochastic .

The condition (32) is satisfied when the time evolution of $`\xi _x`$ is replaced by the stochastic variable $`\xi _k`$ appearing in Eq. (22) with the the stochastic variable $`\xi _k(x_k)`$, which has two distinct values, with the same probability:

$$\xi _k(x_k)\{\begin{array}{c}W_+(x_k)=W\gamma x_k\\ W_{}(x_k)=W\gamma x_k\end{array}$$
(36)
In the time scale $`T_\rho `$, we can replace the variable $`\xi _x(t)`$ of Eq.(28) with its average over the faster bath fluctuations, and thus with

$$\xi _x=\frac{1}{2}[W\gamma x(W+\gamma x)]=\gamma x.$$
(37)
In Fig. 2 we show a sample trajectory corresponding to the dynamic prescriptions of Eq. (36). In the lateral region, after the end of the state that has pushed the trajectory outside the central region, the slope assigned to the trajectory by the ensuing states have always the same sign, either negative, for the top region, or positive, for the bottom region. This is a consequence of the definition itself of side regions and central stripe, the former being those characterized by $`|x|>W/\gamma `$ and the latter referring $`|x|<W/\gamma `$. The final effect is that the trajectories that have exited the central region, are steadily brought back to it. We see that, on the contrary, in the central region the trajectory fluctuates back and forth more or less as it would do in the case of free diffusion. This is perhaps the intuitive reason why the linear response theory of Section III, referring to a central region of extremely large size, makes the system’s statistics fall in the Lévy basin of attraction as it does in the free case .

### C Computer simulation

The computer simulations are done generating the trajectories of Eq.(27) with the assumption that the values $`W_+(x_k)=W\gamma x_k`$ and $`W_{}(x_k)=W\gamma x_k`$ have the same probability. The distribution of sojourn times in one of these two states is given by the function $`\psi (t)`$ of Eq.(26), which, in turn, is realized by a suitable nonlinear transformation of a variable with a uniform distribution in the interval $`[0,1]`$ obtained by a standard routine for the generation of random numbers. It has to be pointed out that the value of either $`W_+(x_k)`$ or $`W_{}(x_k)`$ is determined by the value that the variable $`x`$ has at the moment when a new state is established and it is kept fixed until the end of the sojourn in this state, namely its value is fixed in the interval $`[t_k,t_{k+1}]`$. Each trajectory is generated according to the integral of the equation (28) with the initial condition $`x(0)=0`$ and it is observed until the time $`T_f20/\gamma >>T_R`$. At this time the value is recorded in the form of a histogram. The distribution of the variable $`x`$ at this final time, assumed to correspond to equilibrium, is obtained by generating a sufficiently large number of trajectories and constructing a histogram.

## III Theoretical prediction

We are now in a position to demonstrate that in the case of extremely weak feedback the resulting equilibrium distribution is Lévy. There are no compelling theoretical predictions in the case of strong feedback, except the heuristic arguments used in section I.

According to the strategy adopted in Ref. the theoretical approach to friction implies the treatment of the response of the booster to an abruptly applied external perturbation. Using the notation adopted in Ref.  the solution to the equation of motion

$$\frac{}{t}\rho (\xi ,t)=\mathrm{\Gamma }(K)\rho (\xi ,t),$$
(38)
where $`\mathrm{\Gamma }(K)`$ denotes the operator driving the motion of the probability distribution $`\rho (\xi ,t)`$ in the presence of a perturbation of intensity $`K`$. We assume that this operator can be expressed as the sum of an unperturbed and a perturbed part,

$$\mathrm{\Gamma }(K)=\mathrm{\Gamma }_0+K\mathrm{\Gamma }_1,$$
(39)
where $`\mathrm{\Gamma }_0`$ denotes the operator driving the unperturbed motion of the variable $`\xi `$ and $`K\mathrm{\Gamma }_1`$ denotes the operator corresponding to the perturbation. The intensity of this perturbation can be freely changed by modifying the size of the parameter $`K`$. We apply the perturbation abruptly at time $`t=0`$, to the booster, assumed to be in the equilibrium state, $`\rho _{eq}`$. The perturbation strength is assumed to be so weak as to make it possible to solve Eq.(38) at the first order perturbation,

$$\frac{}{t}\rho _1(\xi ,t)=\mathrm{\Gamma }_0\rho _1(\xi ,t)+K\mathrm{\Gamma }_1\rho _0,$$
(40)
where $`\rho _0`$ and $`\rho _1`$ denote the zero-th and the first-order distribution density, respectively. We assume that the system before the abrupt perturbation is in thermodynamic equilibrium. This means that we identify $`\rho _0`$ with $`\rho _{eq}`$, which is assumed to fulfill the following equation:

$$\mathrm{\Gamma }_0\rho _{eq}=0.$$
(41)
The formal solution to Eq.(40) is

$`\rho _1(\xi ,t)`$ $`=`$ $`K{\displaystyle _0^t}𝑑t^{}exp[\mathrm{\Gamma }_0(tt^{})]\mathrm{\Gamma }_1\rho _0(\xi ,t^{})=`$ (42)
$`=`$ $`K{\displaystyle _0^t}𝑑t^{}exp(\mathrm{\Gamma }_0t^{})\mathrm{\Gamma }_1\rho _0(\xi ,tt^{})=`$ (43)
$`=`$ $`K{\displaystyle _0^t}𝑑t^{}exp(\mathrm{\Gamma }_0t^{})\mathrm{\Gamma }_1\rho _{eq}.`$ (44)
We assume that at equilibrium the mean value of $`\xi `$ vanishes. As a consequence, Eq.(44) yields

$$\xi (t)=K_0^tC(t^{})𝑑t^{},$$
(45)
where

$$C(t)\xi exp(\mathrm{\Gamma }_0t)\mathrm{\Gamma }_1_{eq}.$$
(46)
This concise expression means a trace over $`\xi `$ of the distribution obtained by applying to $`\rho _{eq}`$ first the perturbation operator $`\mathrm{\Gamma }_1`$, then the unperturbed time evolution operator $`exp(\mathrm{\Gamma }_0t)`$ and finally the variable $`\xi `$ itself.

In the case of ordinary statistical mechanics one usually makes the assumption that the variable $`\xi `$ has Gaussian statistics. This means a continuous variable driven by both the unperturbed stochastic environment and by the perturbation of intensity $`K`$. To define the action of the perturbation on the time evolution of the $`\xi `$-trajectory, let us switch off the influence of the stochastic environment. In this case the perturbation action would be expressed by $`d\xi /dt=K`$. Moving from the Heisenberg-like to the Schrödinger-like representation, we would obtain $`\frac{}{t}\rho (\xi ,t)=K\frac{}{\xi }\rho (\xi ,t)`$ which implies that the perturbation operator is

$$\mathrm{\Gamma }_1=\frac{}{\xi }.$$
(47)
Note that in the Gaussian case $`\rho _{eq}(\xi )`$ is a Gaussian distribution proportional to $`exp(\xi ^2/(2\xi ^2_{eq}))`$. Consequently $`C(t)`$ of Eq.(46) becomes

$$C(t)=\mathrm{\Phi }_\xi (t).$$
(48)
Note that

$$dx(t)/dt=\xi (t)$$
(49)
and that, in the absence of perturbations, with all the trajectories starting at $`x=0`$,

$$x^2(t)_0=2\xi ^2_{eq}_0^t𝑑t^{}_0^t^{}𝑑t^{\prime \prime }\mathrm{\Phi }_\xi (t^{\prime \prime }).$$
(50)
We adopt the subscript $`0`$ to indicate that the time evolution of the second moment takes place in the absence of perturbations. Using Eq.(49), Eq.(50), and Eq.(45) supplemented by Eq.(48), and setting equal to zero an arbitrary integration constant so as to ensure a response proportional to the perturbation strength, we arrive at

$$x(t)=Kx^2(t)_0/(2\xi ^2_{eq}).$$
(51)
This equation is a generalized Einstein relation , more recently discussed by Barkai and Fleurov .

It is evident that the ordinary condition, (51), cannot apply in the case here under study, corresponding to the conditions Eqs. (15) and (16). In fact, in this latter case the adoption of the ordinary linear response relation would yield an infinite current. Earlier investigation  has established that the ordinary linear response theory, breaks down, and a new form of linear response to perturbation shows up, which, however, cannot be expressed in terms of unperturbed dynamics. This means that the response depends on the dynamical model adopted and must be established through numerical integration.

In conclusion, using the strategy of Ref. we find that in the presence of feedback the free diffusion process

$$dx/dt=\xi $$
(52)
must be replaced with

$$dx/dt=\xi (t)_x+\xi ,$$
(53)
where we have divided the current, in accordance with linear response theory into the sum of two terms. The zero-th order term is given by the fluctuation $`\xi `$ assumed to be in the same unperturbed condition as in Eq.(52). The first-order term is $`\xi (t)_x`$. Although the variable $`\xi `$ sojourns for long times in one of the two-velocity states, we assume the feedback to be so weak as to be compatible with $`\xi `$ making many jumps from one velocity state to the other while the value of the variable $`x`$ remains essentially unchanged. On the other hand, setting

$$K=\mathrm{\Delta }^2x$$
(54)
where $`\mathrm{\Delta }^2`$ is a constant, we obtain

$$dx/dt=\gamma x(t)+\xi (t),$$
(55)
and

$$\gamma =\mathrm{\Delta }^2\chi ,$$
(56)
the susceptibility $`\chi `$ being defined by

$$\chi _0^{\mathrm{}}C(t)𝑑t.$$
(57)
We are now in a position to address the problem under discussion here via three approximation steps:

(i) We replace Eq.(55) with

$$dx/dt=\gamma x(t)+\eta (t),$$
(58)
where $`\eta `$ is a Lévy noise. This means that the variable $`\eta `$ in one single time step produces jumps proportional to those produced by the variable $`\xi `$ sojourning for a time $`t=\eta /W`$ in one of the two velocity states.

(ii) We make a numerical simulation of Eq. (55)

(iii) We make a completely dynamical treatment of the whole fluctuation-dissipation process.

We shall refer to approximation (i) as the *stochastic approximation*. The advantage of this approximation is that Eq.(55) is made equivalent to an equation studied years ago by West and Seshadri.The stochastic force results in a phase space operator equivalent to a fractional derivative of order $`\alpha `$, with

$$\alpha =\beta +1.$$
(59)
Thus, it is possible to show (see also Refs. for a more recent discussion of the same problem) that the Fourier transform of the distribution density $`\sigma (x,t)`$, $`\widehat{\sigma }(k,t)`$,reads:

$$\frac{}{t}\widehat{\sigma }(k,t)=\left(b|k|^\alpha +\gamma k\frac{}{k}\right)\widehat{\sigma }(k,t),$$
(60)
which yields the following characteristic function

$$\text{ln}\widehat{\sigma }(k,t)=\frac{(1e^{\alpha \gamma t})}{\alpha \gamma }b|k|^\alpha $$
(61)
which for no dissipation has the familiar form

$$\underset{\gamma 0}{lim}\text{ln}\widehat{\sigma }(k,t)=tb|k|^\alpha .$$
(62)
However, asymptotically, in time we obtain the equilibrium condition

$$\widehat{\sigma }(k,\mathrm{})=exp\left(\frac{b|k|^\alpha }{\alpha \gamma }\right)=e^{b_\gamma |k|^\alpha }.$$
(63)
The second equality defines the equilibrium parameter $`b_\gamma `$ of the Lévy characteristic function. Approximation (ii) is referred to as the *dissipative Lévy walk*, and has already been discussed in an earlier work, where it was shown that in the limiting case of $`\gamma 0`$ leads to an equilibrium distribution equivalent to that of Eq.(63).

## IV The numerical simulation

An interesting result of the numerical simulation of the dynamical model is that the equilibrium distribution is not uniquely determined. The condition of weak feedback yields a distribution which has properties different from those of the equilibrium distribution stemming from the condition of strong feedback.

The numerical simulation of the dynamical model can be done both in the case of weak and strong feedback. Interesting new effects are revealed by the simulation of the case of strong feedback. However, these are left as a subject for future theoretical discussions. Here we illustrate only the simulation results concerning the case of weak feedback. According to the program of Ref., we have developed a theory resting on the linear response theory, although in the non-conventional form of Section III, and consequently on the assumption of a very weak feedback. Therefore, the simulation results here discussed refer to a case of feedback so weak as to ensure that the requirements of Section III are fulfilled.

Fig. 3 refers to the Lévy basin of attraction. We notice that the dissipative Lévy walk results in pronounced peaks. These peaks are produced by the fact that the trajectories cannot overcross the values $`x=\pm W/\lambda `$. These peaks signal the region within which a good agreement among the fully dynamic treatment, the dissipative Lévy walk and the stochastic approximation is expected. Note that, as remarked in Section III, the stochastic approximation is equivalent to the equilibrium of the theory of West and Seshadri (WS) , where the parameter $`b_\gamma `$ of Eq. (63) is calculated with the following formula :

$$b_\gamma =\frac{W(\beta \overline{T}W)^\beta }{\gamma (\beta +1)}\mathrm{sin}\left(\frac{\pi }{2}\beta \right)\mathrm{\Gamma }(1\beta ),$$
(64)
where $`\overline{T}`$ and $`\beta `$ are defined in Eq.(23) and $`\mathrm{\Gamma }()`$ is the Euler gamma function. We see that in the region enclosed by the peaks a very good accordance among these three distinct approaches is found indeed.

Fig. 4 illustrates the situation in the Gauss basin of attraction. In this case the side peaks of the dissipative walk disappear. The counterpart of the WS statistics here is the Ornstein-Uhlenbeck statistics , already used in an earlier publication  to predict the form of equilibrium distribution in the Gauss basin of attraction. We use the same prediction where the variance $`\sigma ^2`$ is:

$$\sigma ^2=\frac{\beta W^2\overline{T}}{\gamma (\beta 1)}$$
(65)
and we find that the agreement with the other two treatments is again very good.

We are convinced that the simulation results of Fig.3 are a reliable numerical evidence of the fact that the dynamic approach to non-canonical equilibrium yields the analytical form proposed 18 years ago by West and Seshadri , here referred to as West-Seshadri (WS) statistics. However, Tsallis statistics and WS statistics have in common the analytical shape of the tails, which is an inverse power law. This might suggest that the accordance between our simulation results and Tsallis statistics is as satisfactory as, or more satisfactory than, the accordance with the WS statistics. We now show that it is not so, and that the accordance of our results with Tsallis statistics is much less satisfactory than with the WS statistics. According to the spirit of our dynamic approach, the comparison should be done with the Tsallis form of equilibrium corresponding to a constraint on the first moment of $`|x|`$. This would produce a more significant departure from the WS statistics, and from our simulation results, the WS statistics being, in fact, a Lévy form of equilibrium statistics. The reader can easily convince him/herself about this property observing Fig 1. Thus, we decided to discuss the comparison between our simulation results and Tsallis statistics adopting a condition more favorable to the Tsallis statistics, namely, the analytical proposal of Eq.(7).

In Fig.5 we establish an exact accordance between Tsallis and WS statistics in the region of large distances. We see that the WS statistical prescription yields a very satisfactory accordance with the results of simulation also in the central region, while the Tsallis non-canonical equilibrium does not. In Fig.6 we organize the comparison in such a way as produce the best fitting between Tsallis statistics and our simulation in the central part of the equilibrium distribution. We see that this has the effect of making the discrepancy between Tsallis statistics and WS statistics much worse in the region of large distances. We see that our simulation results are much closer to the WS equilibrium than to the Tsallis equilibrium. We have to point out, furthermore, that the region of large distances of Fig. 6, where the Tsallis form of non-canonical equilibrium apparently yields a better agreement with our numerical results, is probably a non-stationary region, which depends on the observation time. We do not have yet any theory concerning this region but the result of the numerical observation. The numerical simulation shows that this region triggers non-stationary effects, due to the fact that the trajectories reaching large distances from the central diffusion region, illustrated in Fig. 2, tend to escape forever. This escape process gives rise, in a time scale that becomes infinitely large upon decrease of the feedback strength, to a diffusion process, namely to a form of non-equilibrium distribution.

We are convinced that the WS is closer to the results of numerical simulation than the Tsallis non-canonical equilibrium. We see, however, that in the case of Fig. 6 the departure of the Tsallis prediction from the results of the numerical simulation is not so marked as to rule out this theoretical proposal. In conclusion, if we use only fitting arguments, we can select the constraint on the second moment rather than that on the first, as the adoption of theoretical arguments would suggest us to do. Eventually, the situation would appear as favorable to Tsallis statistics as illustrated in Fig. 6.

## V Concluding remarks

The present paper is not just about obtaining the best fit to a numerical simulation. We aim at deriving equilibrium from dynamics with no use of thermodynamic arguments. The result of our investigation is that equilibrium is dictated by the WS statistics rather than by the the Tsallis statistics. This yields a satisfactory agreement with the results of numerical simulation with no fitting parameters. The advocates of Tsallis statistics might reach a result as satisfactory as that provided by the dynamic theory (even if we think that actually it is less satisfactory) only using fitting parameters. This paper results in two important facts. The first is a numerical simulation yielding a non-canonical equilibrium. The second is a theory to account for this non-canonical equilibrium. The theoretical arguments lead us to conclude that the non-canonical equilibrium must correspond to the WS statistics.

If only fitting arguments were to be used, it would be more difficult to make a choice between the Tsallis and the WS statistics. However, if we invoke also theoretical arguments, we find that no room is left for the Tsallis statistics in the case of the dynamic model studied in the present paper. At least within the range of this dynamic model, we show that the non-canonical equilibrium is possible in nature, but it must correspond to the theoretical proposal of West and Seshadri  rather than to the predictions of non-extensive thermodynamics . It is interesting to remark that in the limiting case of very weak feedback our model becomes identical to the Lévy flights subject to the Hookean force of the recent work of Jespersen, Metzler and Fogedby , and that both models are equivalent to that originally studied by West and Seshadri . Also these authors  found the Tsallis statistics to be incompatible with non-Gibbsian nature of the corresponding stationary solution.

In a lucid discussion Lebowitz  has recently restated the point of view of Boltzmann, establishing the microscopic origin of irreversible macroscopic behavior. In his view the adoption of the laws of large numbers is essential, and the role of deterministic chaos becomes important only if it applies to a macroscopic number of non-interacting particles. According to Lebowitz, mixing and ergodicity are notions that are “unnecessary, misguided and misleading.” In other words, this opinion reflects the conviction, mirrored by the handbooks of statistical mechanics, that the unification of mechanics and thermodynamics rests on the condition $`N\mathrm{}`$, where $`N`$ denotes the number of degrees of freedom of the system under study. These statements do not, however, conflict with the results of Ref. reported in Section I B. In fact, the signature of the dynamical origin of thermodynamics, as expressed by Eq.(11), is lost in the limiting case of the thermodynamic limit $`N\mathrm{}`$, thereby making the controversy between the advocates of mixing and the advocates of $`N=\mathrm{}`$ difficult, if not impossible, to substantiate with experimental arguments. This means, in other words, that the canonical equilibrium distribution can be derived using simple arguments, based only on probabilistic concepts and the law of large numbers, or, if we wish, also on the dynamical arguments of Ref.. If the criterion of simplicity is adopted, one might be tempted to choose the former approach which leads to the wanted result with little or no effort, while the second approach yields the same conclusion after many complicated calculations based on assumed dynamical properties, that only in a few cases have been rigorously proved.

The present paper shows how to extend to the case of boosters with no finite time scale the program of Ref. . In this new case the dynamical approach to equilibrium yields a non-canonical equilibrium which is that advocated many years ago by West and Seshadri . As pointed out by the theoretical discussion of Section III, it must be remarked that Eq. (58), yielding the form of equilibrium of West and Seshadri, is the consequence of a form of linear response, departing from the traditional wisdom behind the generalized Einstein relation . It is remarkable that the only plausible form of linear response, resting on dynamics, yields the same prescription as that suggested by a phenomenological approach. Thus, while we do agree with Rajagopal and Abe about the fact that the canonical equilibrium is not the only acceptable form of equilibrium, we depart from them on the specific form that this equilibrium will take, since, as we have seen in Section I A, the Lévy statistics must not be confused with the generalized canonical distribution of Tsallis.


format: text

text:
# Effective Field Theory in Nuclear Many-Body Physics

## I Overview

Reference \[REFERENCES\] is a presentation entitled Relativistic Nuclear Many-Body Theory given at the Seventh International Conference on Recent Progress in Many-Body Theories held in Minneapolis, Minnesota, in August, 1991. This was a report on a long-term effort to understand the nuclear many-body system in terms of relativistic quantum field theories based on hadronic degrees of freedom , a topic we refer to as quantum hadrodynamics (or QHD). An extensive, more recent review of work in this area is contained in Ref. \[REFERENCES\], and a text now exists  that provides background material.***Extensive references to other work in this field are contained in Refs. \[REFERENCES\] through \[REFERENCES\]. There has been significant recent progress in this area , and the goal of this contribution is to summarize briefly what has transpired since the presentation in Ref. \[REFERENCES\].

The only consistent framework we have for discussing the relativistic many-body system is relativistic quantum field theory based on a local lagrangian density. In any lagrangian approach, one must first decide on the generalized coordinates, and hadronic degrees of freedom—baryons and mesons—are the most appropriate for ordinary nuclear systems (QHD). Early attempts involved simple renormalizable models, which reproduced some basic features of the nuclear interaction . The advantage of such models is that in principle, one can consistently investigate and relate all aspects of nuclear structure to a small number of renormalized coupling constants and masses. The disadvantage, in addition to the strong coupling constants that make reliable approximation schemes difficult to come by, is that limiting the discussion to renormalizable lagrangians is too restrictive. Despite these drawbacks, the simple models led to interesting insights. In relativistic mean-field theory (MFT), nuclear densities, the level structure of the nuclear shell model, and the spin dependence of nucleon–nucleus scattering are reproduced . The simplest model (QHD–I) consists of baryons and isoscalar scalar and vector mesons. A basic feature of all these models is that there are strong scalar and vector mean fields present in the nucleus, which cancel in the binding energy but which add to give the large spin-orbit interaction .

There is now overwhelming evidence that the underlying theory of the strong interaction is quantum chromodynamics (QCD), a Yang–Mills non-abelian gauge theory built on an internal color symmetry of a system of quarks and gluons. If mass terms for the $`u`$ and $`d`$ quarks are absent in the lagrangian, QCD possesses chiral symmetry in the nuclear domain; although spontaneously broken in manifestation, this symmetry should play an essential role in nuclear dynamics. The challenge  was to understand the theoretical basis of QHD, the successes that it had, and its limitations, in terms of QCD. Indeed, as we say in our summary in Ref. \[REFERENCES\]:

> More generally, it is probable that at low energies and large distances, QCD can be represented by an effective field theory formulated in terms of a few hadronic degrees of freedom. All possible couplings must be included in the low-energy effective lagrangian, which is then to be used at tree level. The underlying assumption of QHD is that of a local relativistic theory formulated in terms of baryons and the lightest mesons. The theory is assumed to be renormalizable, and one then attempts to extract predictions for long-range phenomena by computing both tree-level diagrams and renormalized quantum loop corrections. In the end, it may turn out that this assumption is untenable, and that the only meaningful interpretation of QHD is as an effective theory, to be used at the tree or one-loop level. The limitation to renormalizable couplings may then be too restrictive. Nevertheless, the phenomenological success of the MFT of QHD–I in the nuclear domain implies that whatever the effective field theory for low-energy, large-distance QCD, it must be dominated by linear, isoscalar, scalar and vector interactions.

The major progress since Ref. \[REFERENCES\], in addition to the multitude of applications discussed in Ref. \[REFERENCES\], has been the following :

* The understanding of QHD as a low-energy, effective lagrangian for QCD, which can be used to improve MFT calculations systematically;
* The understanding of the way spontaneously broken chiral symmetry is realized in QHD;
* The development of a consistent, controlled expansion and approximation scheme that allows one to compute reliable results for bulk nuclear properties;
* The relation of relativistic MFT to density functional theory and Kohn–Sham potentials, placing it on a sounder theoretical basis;
* The understanding of the robustness of many of the QHD–I results.

In section 2 we discuss the relation to density functional theory  and Kohn–Sham potentials . Section 3 contains a brief presentation of the effective lagrangian, and section 4 summarizes some recent results.

## II Density Functional Theory

We begin with a discussion of nonrelativistic density functional theory (DFT) and generalize later to include relativity. The basic idea behind DFT is to compute the energy $`E`$ of the many-fermion system (or, at finite temperature, the grand potential $`\mathrm{\Omega }`$) as a functional of the particle density. DFT is therefore a successor to Thomas–Fermi theory, which uses a crude energy functional, but eliminates the need to calculate the many-fermion wave function.

The strategy behind DFT can be seen most easily by working in analogy to thermodynamics . For a uniform system in a box of volume $`V`$ at temperature $`T`$, one first computes the grand potential $`\mathrm{\Omega }(\mu ,T,V)`$, where $`\mu `$ is the chemical potential. It then follows that the number of particles $`N`$ is determined by

$$N=\widehat{N}=\mathrm{\Omega }/\mu .$$
(1)
The convexity of $`\mathrm{\Omega }`$ implies that $`N`$ is a monotonically increasing function of $`\mu `$, so this relation can be inverted for $`\mu (N)`$. Finally, one makes a Legendre transformation to the Helmholtz free energy $`F(N,T,V)=\mathrm{\Omega }(\mu (N),T,V)+\mu (N)N`$ to discuss systems with a fixed density $`n=N/V`$.

For a finite system, we replace the chemical potential with an external, single-particle potentialIn fact, one can absorb $`\mu `$ into the definition of $`v`$. We suppress all spin dependence at this point. $`_iv(𝐫_i)`$. The grand potential is now a functional: $`\mathrm{\Omega }([v(𝐫)],T)`$, and a functional derivative with respect to $`v`$ gives the particle density:Higher variational derivatives yield various correlation functions.

$$n(𝐫)=\widehat{n}(𝐫)=\delta \mathrm{\Omega }/\delta v(𝐫).$$
(2)
The convexity of $`\mathrm{\Omega }`$ allows us (in principle) to invert this relation and find $`v(𝐫)`$ as a (complicated) functional of $`n(𝐫)`$. Finally, we make a functional Legendre transformation to define the Hohenberg–Kohn free energy, which is a functional of $`n(𝐫)`$:

$$F_{\mathrm{HK}}[n(𝐫)]=\mathrm{\Omega }[v(𝐫)]d𝐫n(𝐫)v(𝐫).$$
(3)
($`T`$ is suppressed.) The variational derivative of this free energy functional with respect to $`n`$ now gives

$$\delta F_{\mathrm{HK}}/\delta n(𝐫)=v(𝐫).$$
(4)
If we now restrict consideration to $`T=0`$ and $`v(𝐫)=0`$, then the Hohenberg–Kohn theorem follows : If the functional form of $`F_{\mathrm{HK}}[n(𝐫)]`$ is known exactly, the ground-state expectation value of any observable is a unique functional of the exact ground-state density. Moreover, it follows immediately from Eq. (4) that the exact ground-state density can be found by minimizing the energy functional. Although we have assumed here that the ground state is non-degenerate, this assumption can be easily relaxed .

The generalization of DFT to relativistic systems is straightforward . The energy functional $`F_{\mathrm{HK}}`$ now becomes a functional of both scalar and vector densities (or more precisely, vector four-currents). Extremization of the functional gives rise to variational equations that determine the ground-state densities.

Significant progress in solving these equations was made by Kohn and Sham , who introduced a complete set of single-particle wave functions. In our case, these wave functions allow us to recast the variational equations as Dirac equations for occupied orbitals. The single-particle hamiltonian contains local, density-dependent, scalar and vector potentials, even when the exact energy functional is used. Moreover, one can introduce auxiliary (scalar and vector) fields corresponding to the local potentials, so that the resulting equations resemble those in a relativistic MFT calculation .

The strength of the approach rests on the following theorem:

> The exact ground-state scalar and vector densities, energy, and chemical potential for the fully interacting many-fermion system can be reproduced by a collection of (quasi)fermions moving in appropriately defined, self-consistent, local, classical fields.

The proof is straightforward . Start with a collection of noninteracting fermions moving in an externally specified, local, one-body potential. The exact ground state for this system is known: just calculate the lowest-energy orbitals and fill them up.§§§For simplicity, we assume that the least-bound orbital is completely filled, so the ground state is non-degenerate. Therefore, if one can find a suitable local, one-body potential based on an exact energy functional, the exact ground state of that system can be determined. But this potential is precisely what one obtains by differentiating the interaction parts of $`F_{\mathrm{HK}}`$ with respect to $`n(𝐫)`$ . The resulting one-body potential will generally be density dependent and thus must be determined self-consistently.

Several points are noteworthy. As noted by Kohn , the single-particle basis constructed as described above can be considered “density optimal”, in contrast to the Hartree (or Hartree–Fock) basis, which is “total-energy optimal”. Thus the exact scalar and vector densities are given by sums over the squares of the Dirac wave functions, with unit occupation probability. Moreover, since these densities are guaranteed to make the energy functional stationary \[the external $`v(𝐫)=0`$\], the exact ground-state energy is also obtained. The proof that the eigenvalue of the least-bound state is exactly the Fermi energy is given in Ref. \[REFERENCES\]. Note, however, that aside from this association, the exact Kohn–Sham wave functions (and remaining eigenvalues) have no known, directly observable meaning.

If one knows the exact functional form of the energy on the density, one can describe the observables noted in the theorem exactly (and easily) in terms of the Kohn–Sham basis. Observables of this type are typically the ones calculated in relativistic MFT. Moreover, it has been known for many years  that the mean-field contributions dominate the single-particle potentials at ordinary densities. Thus, by parametrizing the energy functional in a mean-field (or “factorized”) form, and by fitting the parameters to empirical bulk and single-particle nuclear data, one should obtain an excellent approximation to the exact energy functional in the relevant density regime. This is the key to the success of relativistic MFT calculations, as we will verify below, using the effective lagrangian constructed in the next section.

## III Effective Lagrangian

We cannot give a detailed derivation and discussion of the effective lagrangian of QHD in this short article, but we can illustrate the basic principles. To exhibit how spontaneously broken chiral symmetry is incorporated quite generally into the hadronic theory, consider the linear $`\sigma `$-model with an additional linear coupling of an isoscalar $`V^\mu `$ to the baryon current, the so-called “chiral $`(\sigma ,\omega )`$ model” . Define right- and left-handed nucleon fields by $`\psi _{R,L}(1\pm \gamma _5)\psi /2`$ and the $`SU(2)`$ matrix $`\text{\$̱mkern4.0muU\$}\mathrm{exp}(i\stackrel{}{\tau }\stackrel{}{\pi }/s_0)`$, where $`\stackrel{}{\pi }`$ is the isovector pion field. If $`M`$ is the nucleon mass, determined by the spontaneous breaking of chiral symmetry, and $`s_0=M/g_\pi `$, then the lagrangian for the chiral $`(\sigma ,\omega )`$ model can be written as the $`s_0\mathrm{}`$ limit of the following generalized lagrangian

$``$ $`=`$ $`i\left[\overline{\psi }_R\gamma _\mu (^\mu +ig_\mathrm{v}V^\mu )\psi _R+\overline{\psi }_L\gamma _\mu (^\mu +ig_\mathrm{v}V^\mu )\psi _L\right]`$ (7)
$`g_\pi s_0\left(1{\displaystyle \frac{\sigma }{s_0}}\right)\left[\overline{\psi }_R\text{\$̱mkern4.0muU\$}^{}\psi _L+\overline{\psi }_L\text{\$̱mkern4.0muU\$}\psi _R\right]+{\displaystyle \frac{1}{2}}(_\mu \sigma ^\mu \sigma )+{\displaystyle \frac{1}{4}}s_0^2\mathrm{tr}(_\mu \text{\$̱mkern4.0muU\$}^\mu \text{\$̱mkern4.0muU\$}^{})`$
$`𝒱(\text{\$̱mkern4.0muU\$},_\mu \text{\$̱mkern4.0muU\$};\sigma )+{\displaystyle \frac{1}{4}}m_\pi ^2s_0^2\mathrm{tr}(\text{\$̱mkern4.0muU\$}+\text{\$̱mkern4.0muU\$}^{}2){\displaystyle \frac{1}{4}}F_{\mu \nu }F^{\mu \nu }+{\displaystyle \frac{1}{2}}m_\mathrm{v}^2V_\mu V^\mu .`$
For $`m_\pi ^2=0`$, this lagrangian is evidently invariant under chiral $`SU(2)_L\times SU(2)_R`$ transformations of the form ($`\sigma `$ and $`V^\mu `$ are unchanged)

$`\psi _L\text{\$̱L\$}\psi _L,\psi _R\text{\$̱R\$}\psi _R,\text{\$̱mkern4.0muU\$}\text{\$̱L\$}\text{\$̱mkern4.0muU\$}\text{\$̱R\$}^{}.`$ (8)
Here $̱L$ and $̱R$ are independent, global $`SU(2)`$ matrices, and the generalized potential $`𝒱`$ is chosen to be invariant, with the limit $`𝒱m_\sigma ^2\sigma ^2/2+O(1/s_0)`$. Conventional notation is recovered with the identification

$`s_0=M/g_\pi f_\pi .`$ (9)
The change of variables $`\text{\$̱mkern4.0muU\$}\text{\$̱xi\$}\text{\$̱xi\$},N_L\text{\$̱xi\$}^{}\psi _L,N_R\text{\$̱xi\$}\psi _R`$ reduces the fermion terms in the preceding lagrangian to

$`_{\mathrm{fermion}}`$ $`=`$ $`\overline{N}\left[i\gamma ^\mu (_\mu +i\text{\$̱v\$}_\mu +ig_\mathrm{v}V_\mu )+\gamma ^\mu \gamma _5\text{\$̱a\$}_\mu M+g_\pi \sigma \right]N,`$ (10)
$`\text{\$̱v\$}_\mu `$ $``$ $`{\displaystyle \frac{i}{2}}(\text{\$̱xi\$}^{}_\mu \text{\$̱xi\$}+\text{\$̱xi\$}_\mu \text{\$̱xi\$}^{}),\text{\$̱a\$}_\mu {\displaystyle \frac{i}{2}}(\text{\$̱xi\$}^{}_\mu \text{\$̱xi\$}\text{\$̱xi\$}_\mu \text{\$̱xi\$}^{}).`$ (11)
This lagrangian is invariant under the following nonlinear chiral transformation:

$$\text{\$̱xi\$}(x)\text{\$̱L\$}\text{\$̱xi\$}(x)\text{\$̱h\$}^{}(x)\text{\$̱h\$}(x)\text{\$̱xi\$}(x)\text{\$̱R\$}^{},N(x)\text{\$̱h\$}(x)N(x),$$
(12)
where $`\text{\$̱h\$}(x)`$ is a local $`SU(2)`$ matrix. It follows that $̱\mkern4.0muU$ still transforms globally according to Eq. (8). Additional mesons and interactions can now be introduced requiring only invariance under the local isospin transformations of Eq. (12). While illustrated within the framework of a simple model, this nonlinear realization of $`SU(2)_L\times SU(2)_R`$ is, in fact, quite general, and can be used as a basis for constructing the most general QHD lagrangian .

The effective lagrangian, which reflects the underlying spontaneously broken chiral symmetry of QCD, and from which the energy functional of the previous section is obtained, is constructed from the following series of steps :

1. A baryon field and low-mass meson fields that concisely describe the important interaction channels, namely, $`\pi (0^{},1),\varphi (0^+,0),V_\mu (1^{},0)`$, and $`\rho _\mu (1^{},1)`$, are the generalized coordinates of choice. The pion, a Goldstone boson, is treated as in the example above. Higher mass meson fields are assumed to be “integrated out” and their contributions contained in the effective coupling constants.
2. Dimensional analysis is first used to characterize the various terms in the effective lagrangian. Briefly, this is done as follows. The initial couplings of the meson fields to the baryon fields are linear, with a strong coupling constant $`g`$. The dimensionless form of this combination is $`g\varphi /M=\varphi /f_\pi `$ \[see Eq. (9)\]; non-Goldstone boson fields are assumed to enter in this dimensionless form. From the mass term of the meson fields $`m^2\varphi ^2`$, with $`m^2M^2`$, one then deduces an overall scale factor in the lagrangian density of $`f_\pi ^2M^2`$. From the baryon mass term $`M\overline{\psi }\psi `$, one concludes that the appropriate dimensionless form of the baryon densities is $`\overline{\psi }\psi /Mf_\pi ^2`$. This “naive” dimensionless analysis (NDA) then implies that, after appropriate combinatorial factors are included, the various terms in the effective lagrangian enter with dimensionless coefficients of order unity.
3. The various interaction terms allowed by the $`SU(2)_L\times SU(2)_R`$ symmetry of QCD are then constructed using the nonlinear realization of chiral symmetry illustrated above. Simply writing down all possible terms does not get one very far unless there is an organizational principle, and the following provides the crucial insight:
4. Although the mean scalar and vector field energies are large compared to the nuclear binding energy, the dimensionless combinations $`g_\mathrm{s}\varphi _0/M\varphi _0/f_\pi `$ and $`g_\mathrm{v}V_0/MV_0/f_\pi `$ are roughly 1/3 and thus provide convenient expansion parameters. Furthermore, spatial variations of the meson fields and of the baryon densities in the nucleus are observed to occur over the scale of the nuclear surface region, and hence the dimensionless ratio $`/M`$ also provides a useful expansion quantity (as does the characterization of chiral symmetry violation at the lagrangian level, $`m_\pi /M`$).
5. A combination of these observations allows one to construct a hierarchy of decreasing contributions to the effective lagrangian for the nuclear many-body system characterized by an integer $`\nu `$ defined by

$`\nu =d+{\displaystyle \frac{n}{2}}+b,`$ (13)
where $`d`$ is the number of derivatives, $`n`$ is the number of nucleon fields, and $`b`$ is the number of non-Goldstone boson fields present in the interaction term. The effective lagrangian at various levels of $`\nu `$ is given in Refs. \[REFERENCES,REFERENCES\].The extension of the effective lagrangian to include electromagnetic interactions as an expansion in powers of derivatives is also discussed in these references.

The effective lagrangian with mean meson fields then determines the energy functional of the previous section, and a representation in terms of Dirac–Hartree orbitals leads to local, nonlinear Hartree equations, which can be solved numerically. The extent to which nuclei exhibit this hierarchy of interactions, and to which this effective lagrangian indeed describes the nucleus, is discussed in the next section.

## IV Results and Summary

The QHD mean-field energy functional discussed above is given in several different forms in Refs. \[REFERENCES,REFERENCES,REFERENCES,REFERENCES\], as are the field equations that result from extremization. The equations are solved self-consistently for the closed-shell nuclei 16O, 40Ca, 48Ca, 88Sr, and 208Pb, and also in the nuclear matter limit. The parameters are then best-fit to empirical properties of the charge densities, the binding energies, and various splittings between energy levels near the Fermi surface using a figure of merit ($`\chi ^2`$) defined by a weighted, squared deviation between the 29 calculated and empirical values. When working at the highest order of truncation (essentially $`\nu =4`$), the calculated results are very accurate, as we illustrate shortly, but they are too numerous to reproduce here .

The critical question is whether the hierarchal organization of interaction terms is actually observed. This is illustrated in Fig. 2, where the nuclear matter energy/particle is shown as a function of the power of the mean fields, which is called $`b`$ in Eq. (13). (There are no gradient contributions in nuclear matter and $`\widehat{\pi }=0`$.) The crosses and error bars are estimates based on NDA and naturalness, that is, overall coefficients are of order unity. It is clear that each successive term in the hierarchy is reduced by a roughly factor of five, and thus for any reasonable desired accuracy, the lagrangian can be truncated at a low value of $`\nu `$. Derivative terms and other coupling terms are analyzed in Ref. \[REFERENCES\], with similar conclusions.

The quality of the fits to finite nuclei and the appropriate level of truncation is illustrated in Fig. 2, where the figure of merit is plotted as a function of truncation order and of various combinations of terms retained in $``$. The full calculations ($``$) retain all allowed terms at a given level of $`\nu `$, while the other two choices keep only the indicated subset. There is clearly a great improvement in the fit (more than a factor of 35) in going from $`\nu =2`$ to $`\nu =4`$, but there is no further improvement in going to $`\nu =5`$. Speaking chronologically, the $`\nu =2`$ results show the level of accuracy obtained almost 20 years ago, while the $`\nu =4`$ results were obtained four years ago . Moreover, the $`\varphi ^n`$ only results at $`\nu =4`$ (   ) show the state of the situation in the late 1980s, as discussed in Ref. \[REFERENCES\]. Recent work  shows that the full complement of parameters at order $`\nu =4`$ is underdetermined, and that only six or seven are determined by this data set, which explains the great success of these earlier models with a restricted set of parameters.

In summary, the hadronic theory of QHD is truly a manifestation of low-energy, strong-coupling QCD. The modern viewpoint of QHD based on effective field theory and density functional theory explains the accurate description of bulk and single-particle nuclear properties. Corrections to the mean-field parametrization of the energy functional can be calculated systematically using the effective hadronic lagrangian . An important goal for the future is finding an efficient, tractable, nonperturbative way to match the low-energy, strong-coupling, effective field theory of QHD to the underlying QCD lagrangian.

## Acknowledgements

We thank Dick Furnstahl for useful comments. This work was supported in part by the US Department of Energy under contracts DE–FG02–87ER40365 and DE–FG02–97ER41023.


text:
# Modelling sublimation and atomic layer epitaxy in the presence of competing surface reconstructions

## Abstract

We present a solid-on-solid model of a binary AB compound, where atoms of type A in the topmost layer interact via anisotropic interactions different from those inside the bulk. Depending on temperature and particle flux, this model displays surface reconstructions similar to those of (001) surfaces of II-VI semiconductors. We show, that our model qualitatively reproduces many of the characteristic features of these materials which have been observed during sublimation and atomic layer epitaxy. We predict some previously unknown effects which might be observed experimentally.
Keywords: Monte Carlo simulations, Surface relaxation and reconstruction, Models of non-equilibrium phenomena, Evaporation and sublimation, Growth

Apart from a few exceptions (eg. ), theoretical models of crystal surfaces in nonequilibrium situations have not taken into account the effects of surface reconstructions, so far. Consequently, a theoretical understanding of the interplay of this phenomenon with the dynamics of growing or sublimating surfaces has not yet been achieved. In this paper we will address this important question in the context of a compound material that displays a competition of different vacancy structures in the terminating layer.

Within the last years, there have been intensive experimental studies concerning the surfaces of II-VI semiconductors, see eg.  for an overview. These studies have revealed the existence of a fairly small number of surface reconstructions, which makes them promising candidates for a theoretical modelling. In the following, we will refer mostly to experimental results for CdTe(001) \[2-9\]. Under vacuum conditions, the CdTe(001) surface is metal terminated. The Cd atoms arrange in a vacancy structure, where half of the potential Cd sites are empty. At low temperatures, a $`c(2\times 2)`$ reconstruction is observed, in which Cd atoms arrange in a checkerboard like configuration. Often a contribution of $`(2\times 1)`$ arrangement in rows along the (110) direction is observed. At a temperature of $`300^0C`$, a reordering occurs, above which a $`(2\times 1)`$ ordering dominates on a sublimating surface. An additional Cd flux stabilizes the $`c(2\times 2)`$ ordering even at high temperatures. Under a Te flux, the surface is Te-terminated, with a $`(2\times 1)`$ reconstruction. At small Te fluxes, the Te-coverage is 1, at low temperatures and high Te fluxes a higher coverage of $`1.5`$ is observed . In this Letter, we present a 2+1-dimensional solid-on-solid model of a binary compound material to get insight into the following problems: (1) Is the claim  of an effective thermal equilibrium of the surface layer justified under the non-equilibrium conditions of sublimation? (2) What is the effect of an external particle flux on the surface reconstruction? (3) How do reconstructions influence techniques of crystal growth like atomic layer epitaxy (ALE)?

In the terminating layer of a CdTe(001) surface, the potential Cd sites form a regular square lattice, where the simultaneous occupation of nearest neighbour (NN) sites in the $`[0\overline{1}0]`$-direction (denoted as y in the following) is forbidden by electron counting rules . In our model, there is an attractive interaction between nearest neighbours in $`[110]`$(x)-direction (coupling $`\epsilon _x`$) and diagonal neighbours (coupling $`\epsilon _d`$). Measuring energy in units of $`|\epsilon _d|`$, we set $`\epsilon _d=1`$. A model, where the underlying crystal was assumed to be fixed, has been studied elsewhere . If the energy difference between a perfect $`c(2\times 2)`$ arrangement and a perfect $`(2\times 1)`$ arrangement is small and positive ($`2<\epsilon _x1.9`$), this model explains the $`c(2\times 2)`$-$`(2\times 1)`$ reordering as an accompanying effect of an order-disorder transition in thermal equilibrium. At low temperature $`T`$, the system is in an ordered phase, where $`c(2\times 2)`$ dominates, above the transition it is in a globally disordered phase, where the $`(2\times 1)`$ rows dominate the local environment of a particle.

In this paper, we extend this planar lattice gas to a model of a three-dimensional crystal. Currently no values of diffusion energy barriers from ab initio or semi-empirical calculations are available, and the existing experimental data are insufficient for a parameter fit. Therefore, a realistic modelling of CdTe is beyond the scope of this work, which aims at an understanding of fundamental properties of nonequilibrium crystal surfaces. Instead, we consider a model of a compound “AB” with a cubic lattice and a comparatively simple potential energy surface which, although not being quantitative, reproduces essential features observed in experiments on CdTe(001). Here, the A (B) atoms correspond to Cd (Te). Within the framework of a solid-on-solid model, we describe the crystal by a two-dimensional array $`\{h_{xy}\}_{x,y=1}^N`$ of integers which denote the height of a column of atoms. We use periodic boundary conditions. The interactions between the particles are described by phenomenological parameters which are supposed to include the effects of surface strain and other interactions. To model the layered structure of the CdTe crystal, we assign the odd heights to A particles, and the even heights to B particles. Inside the bulk of the crystal, there is an attractive interaction between the particles and their NNs which is isotropic in directions parallel to the surface. While there is no difference in the interaction of B particles between the bulk and the surface, A particles on the surface interact with the anisotropic interactions $`\epsilon _x`$ and $`\epsilon _d`$. The other parameters of the model are defined in Figure 1. Thus, the Hamiltonian of the crystal is

$$H=\epsilon _{AB}n_{AB}+\epsilon _{BB}n_{BB}+\epsilon _bn_{AA}^b+\epsilon _hn_{AA}^h+\epsilon _dn_{AA}^d+\epsilon _xn_{AA}^x$$
(1)
Here, $`n_{AB}`$ is the number of AB bonds between the layers, $`n_{BB}`$ is the number of NN pairs of B atoms, $`n_{AA}^b`$ the number of A NN pairs inside the bulk, $`n_{AA}^h`$ the number of surface A atoms next to a higher column, $`n_{AA}^d`$ ($`n_{AA}^x`$) the number of bindings of surface A atoms to diagonal neighbours (NNs in x-direction) at the same height. The stay of particles in the wrong sublattice and the formation of NN pairs of surface A atoms in y-direction is forbidden. We allow diffusion both to nearest neighbour sites and to diagonal neighbours. Additionally, we permit the diffusion of an AB pair with the B atom on top of the A atom, if the diffusion of the B atom alone would end in the wrong sublattice. This process is required to preserve the ergodicity of the system. We simulate a Kawasaki dynamics, i.e. the rate of diffusion events is $`\nu _0\mathrm{exp}((B_0+\mathrm{\Delta }H)/T)`$ if $`\mathrm{\Delta }H>0`$ and $`\nu _0\mathrm{exp}(B_0/T)`$ else. The desorption of an A (B) particle requires an additional activation energy $`E_A`$ ($`E_B`$). This is done using rejection free Monte-Carlo techniques employing a binary search tree . Particles from an external source arrive at each lattice site with equal probability. If the adsorption of the particle at this site would lead to a forbidden state, it is destroyed.

This dynamics does not depend on the parameter $`\epsilon _{AB}`$, so it needs not to be specified, and $`\epsilon _{BB}`$ and $`\epsilon _b`$ enter only via the sum $`\epsilon ^{}:=\epsilon _{BB}+\epsilon _b`$. The constant prefactor $`1/t_0:=\nu _0\mathrm{exp}(B_0/T)`$ sets the timescale of the model and needs not to be specified, if we measure time in units of $`t_0`$.

In the following, we will use the parameter set $`\epsilon ^{}=\epsilon _d=1`$, $`\epsilon _h=0`$, $`\epsilon _x=1.97`$, $`E_A=E_B=1.5`$ (in units of $`|\epsilon _d|`$). The fact that the difference in the surface energies between perfectly $`c(2\times 2)`$ and $`(2\times 1)`$ reconstructed surfaces is small is consistent with DFT calculations . However, precise values of binding energies as a quantitative input to our model are not available. We have verified, that the behaviour of our model remains qualitatively the same for a wide range of parameters, as long as the basic property of a small energy difference between the reconstructions is preserved.

To characterize the surface reconstruction quantitatively, we measure the A-coverage $`\theta _A=n_A^{surf.}/N^2`$, and the normalized correlations $`C_{AA}^d:=n_{AA}^d/(4n_A^{surf.})`$ and $`C_{AA}^x:=n_{AA}^x/(2n_A^{surf.})`$, which quantify the fraction of A-atoms on the surface which are incorporated in locally $`c(2\times 2)`$ or $`(2\times 1)`$ reconstructed areas, respectively. Additionally, we calculate LEED intensities to model diffraction experiments. Neglecting multiple scattering, this is done by calculating the fourier transform of the array $`\{\mathrm{exp}(i\mathrm{\Delta }kh_{xy})\}_{x,y=1}^N`$, where equal diffraction intensities for A and B atoms have been assumed. The antiphase condition between different layers is $`\mathrm{\Delta }k=\pi /2`$.

Figure 2a shows the evolution of an initially B-terminated, flat surface under vacuum at a temperature $`T=0.5`$. Initially, $`\theta _A`$ increases from zero to an asymptotic value close to 0.5. We find, that this increase can be fitted with an exponential relaxation, $`\theta _A=\theta _A^{\mathrm{}}\left(1\mathrm{exp}(t/\tau )\right)`$. Here, $`\tau `$ is the time constant of the decay of the B-terminated surface. We find, that the temperature dependence of the corresponding desorption probability $`p1/\tau `$ follows an Arrhenius law, $`p\mathrm{exp}(E_{act.}/T)`$, with an activation energy $`E_{act.}=B_0+3.0\pm 0.1`$ After this onset, A and B evaporate stoichiometrically, and sublimation proceeds in layer-by-layer mode, which reflects in the oszillations of the specular LEED intensity $`I_{spec}`$ in antiphase (dotted curve in Figure 2a). The dependence of the sublimation rate on $`T`$ follows an Arrhenius law with a higher activation energy $`B_0+6.0\pm 0.2`$. Thus, our model reproduces the observation of Cd-terminated CdTe surfaces under vacuum. Studies of surface lifetimes  and QMS measurements of evaporation rates  have also revealed significantly smaller activation energies for the decay of the Te-terminated surface ($`0.7eV`$) than for congruent evaporation ($`1.94eV`$ in layer-by-layer mode). We note however, that in our model the microscopical energy barriers for desorption of A and B are equal. Therefore the difference in the macroscopic activation energies is solely due to the stabilizing effect of the surface reconstruction.

Surprisingly, $`C_{AA}^d`$ and $`C_{AA}^x`$ oszillate during sublimation, which proceeds layer by layer. Each time a complete layer has desorbed and the surface is atomically flat (maxima of $`I_{spec}`$) the $`c(2\times 2)`$ reconstructed fraction of the surface is maximal. On the contrary, a rough surface with a large number of islands (minimal $`I_{spec}`$) seems to prefer the $`(2\times 1)`$ reconstruction. This can be understood from the fact that the attractive lattice gas interactions are present only between particles in the same layer. Thus, the island edges impose open boundary conditions to the lattice gas of A atoms on the island. In contrast to a $`c(2\times 2)`$ reconstructed domain, a $`(2\times 1)`$ terminated island can reduce the energy of its boundary by elongating in x-direction. Since the ground state energies of both structures are nearly degenerate, the formation of $`(2\times 1)`$ may reduce the surface free energy. This picture is confirmed by the fact that islands on sublimating surfaces are indeed elongated (See Figure 4a,b) .

Experiments have revealed a high density of steps on CdTe surfaces due to the intrinsic surface morphology , which leads to a dominant contribution of step flow sublimation. To incorporate this effect in our model, we have performed simulations of vicinal surfaces . We find, that the oszillations in the correlations dissappear for terrace widths smaller than $`50`$ lattice constants, where step flow sublimation dominates. Then, the correlations become stationary apart from statistical fluctuations. Figure 3a shows the $`T`$-dependence of $`\theta _A`$, $`C_{AA}^d`$, $`C_{AA}^x`$ and the sublimation rate $`r_{sub}`$ at a step distance of $`32\sqrt{2}`$ lattice constants. The steps were oriented at an angle of 45o to the x-axis, which corresponds to the preferential (100)-orientation of steps observed on CdTe(001) . Snapshots of the surfaces are shown in Figure 4. Clearly, at $`T=0.55`$ there is a transition from a $`c(2\times 2)`$ configuration at low temperature to a high-temperature regime where the $`(2\times 1)`$ ordering dominates and the material is sublimating. This is the conterpart of the phase transition observed in our two-dimensional lattice gas model. In the investigated temperature range $`\theta _A`$ decreases only slightly from $`\theta _A=0.498`$ at $`T=0.4`$ to $`\theta _A=0.41`$ at $`T=0.8`$.

To get insight into the behaviour of a material in an MBE environment, it is important to understand how the structure of the surface changes when it is exposed to a particle beam. To this end, we have simulated surfaces which were exposed to a flux of pure A or B. Figure 3b,c shows the quantities $`\theta _A`$, $`C_{AA}^d`$ and $`C_{AA}^x`$ as functions of the adsorption rate at a temperature $`T=0.4`$, which is far below the $`c(2\times 2)`$-$`(2\times 1)`$ reordering under vacuum, and at $`T=0.57`$, which is slightly above the transition. Applying B fluxes, we can regulate $`\theta _A`$ to values smaller than those close to 0.5 observed under vacuum. We find, that a decrease of $`\theta _A`$ leads to a significant increase in $`C_{AA}^x`$ (Fig. 3b,c). This effect is analogous to the higher fraction of $`(2\times 1)`$ ordering on surfaces with a large number of islands: at low $`\theta _A`$, A-covered islands form on the surface, whose boundaries favour the arrangement of A atoms in rows. On the contrary, an A flux increases $`\theta _A`$ and $`C_{AA}^d`$ and leads to the formation of a dominant $`c(2\times 2)`$ reconstruction even at temperatures above the reordering under vacuum. The reappearance of a $`c(2\times 2)`$ reconstruction under a Cd-flux at high temperature and the formation of Te-terminated surfaces with coverage 1 even at low Te fluxes are well known from experiments on CdTe . In these experiments however, the formation of Te dimers on the surface was observed, an effect which is not present in our simple model assuming isotropic interactions between B atoms.

It has been argued , that the CdTe surface is close to thermal equilibrium during sublimation, since the bulk serves as a particle reservoir which defines chemical potentials for both elements. Comparing our solid-on-solid model with the two-dimensional anisotropic lattice gas, we calculate the equilibrium values of $`C_{AA}^d`$, $`C_{AA}^x`$ in the planar lattice gas model at the value of $`\theta _A`$ we measure for AB at a given temperature, which can be done using transfer matrix techniques . Although the non-equilibrium conditions of sublimation enhance the dominance of $`C_{AA}^x`$ over $`C_{AA}^d`$ in the high temperature regime, we find qualitative agreement in the behaviour of $`C_{AA}^d`$ and $`C_{AA}^x`$. However, a simple mapping of the surface layer of the sublimating system on a 2D lattice gas in thermal equilibrium is not possible.

Atomic layer epitaxy (ALE) provides a scenario for the investigation of our model of the II-VI semiconductor surface in a situation present in technical applications. The idea is to obtain self-regulated growth by alternate deposition of pure A and B. In the absence of reconstructions, one would expect the formation of complete monolayers of A and B during deposition of the elements, yielding growth at a speed of one monolayer (ML) per cycle. However, this does not necessarily apply in the presence of surfaces terminated by vacancy structures with submonolayer coverage. Experimentally, in CdTe one finds growth rates of $``$ 1 ML/cycle only at $`T<260^oC`$, and growth at a speed of $``$ 0.5 ML/cycle at higher temperatures . Figure 2b shows the evolution of $`\theta _A`$, $`C_{AA}^d`$ and $`C_{AA}^x`$ during ALE at $`T=0.3`$. Each cycle consists of two phases of length $`t_{cycl}=410^3t_0`$. In the first phase, an A flux of $`510^3\text{ML}/t_0`$ is applied, in the second phase B is deposited at the same rate. Growth was started from a B terminated flat surface. During the first A phase the surface becomes A-terminated with $`\theta _A1/2`$. Since $`T`$ is far below the transition, the reconstruction is $`c(2\times 2)`$. At the onset of the B-phase, $`\theta _A`$ decreases rapidly. This leads to an increase of $`C_{AA}^x`$, as reported above for surfaces under a small stationary B-flux. At the end of the cycle, $``$50% of the surface are covered with B-terminated islands, since only half of the A atoms needed for a closed monolayer is present. The following A-phase deposits half a monolayer of A again. However, now the reconstruction is preferentially $`(2\times 1)`$. This is analogous to the maximum of $`C_{AA}^x`$ we observe in layer-by-layer sublimation after half a monolayer has desorbed, which is due to the influence of the island edges on the reconstruction. During the following B phase, A atoms and AB pairs on top of the islands diffuse into the gaps between the islands, which leads to the formation of a closed B-terminated monolayer. In our simulations, we observe 2-3 repetitions of this cycle . There are two effects, which keep it from repeating infinitely. First, we obtain a growth rate of 0.44 ML/cycle, which is smaller than the ideal value of 0.5 ML/cycle. Second, the diffusion of particles from the islands into the gaps is not complete. Both effects hinder a perfect closure of the monolayer at the end of the even cycles, which leads to a damping of the characteristic behaviour. These observations closely resemble the ideas developed to explain CdTe growth at a speed of $``$0.5 ML/cycle in .

The experimental observation of growth rates of 1 ML/cycle in ALE at low temperatures is inconsistent with our model in its present form, since the hardcore repulsion between the A atoms hinders the deposition of more than 0.5 ML of A. However, this effect might be explained with the existence of weakly bound precursor states, where particles could stay close to the surface until the deposition of the opposite species allows for an incorporation into the crystal. At higher temperatures this reservoir would be inactive due to the weakness of the binding. Future research might also address the influence of an interlayer diffusion barrier on ALE and MBE in our model. Additionally, one further step towards a more realistic modelling of CdTe should be the simulation of a zinc-blende lattice which, however, requires the specification of a greater number of parameters.

In summary, we have presented a simple model of a binary compound material which reproduces many effects observed in experiments on CdTe(001) qualitatively. We conclude suggesting some experiments to test our model. (1) The preparation of CdTe(001) surfaces with an extremely low intrinsic step density might allow for the observation of oscillations in the correlations during layer-by-layer sublimation. (2) The deposition of small amounts of Te on a CdTe(001) surface should lead to a preferential arrangement of the remaining Cd atoms in a $`(2\times 1)`$ ordering. (3) During Cd deposition in ALE of CdTe a $`c(2\times 2)`$ reconstruction is expected in cycles starting from flat Te-terminated surfaces, while in cycles starting from rough surfaces $`(2\times 1)`$ should be present.

We thank W. Kinzel and M. Sokolowski for stimulating discussions and a critical reading of the manuscript. M.A. was supported by the Deutsche Forschungsgemeinschaft.