# A Monte Carlo Method for the Numerical Simulation of Tsallis Statistics

## 1 Introduction

In 1988 C. Tsallis proposed a thermostatistics formalism based on a non–extensive entropy definition. In the most recent formulation of the Tsallis Statistics (TS) in the canonical ensemble, at fixed temperature $`T`$, observables $`𝒪_q`$ are obtained as averages of microscopic functions $`O_i`$:

$$𝒪_q=\underset{i=1}{\overset{W}{}}O_iP_i$$
(1)   
while the entropy is given by

$$S_q=\frac{(_{i=1}^WP_i^{1/q})^q1}{1q}$$
(2)   
In these expressions, $`P_i`$ are the escort probabilities for configuration $`i=1,\mathrm{},W`$ with energy $`\epsilon _i`$:

$$P_i=\frac{A_i}{_{j=1}^WA_j}\frac{[1(1q)\epsilon _i/T^{}]^{\frac{q}{1q}}}{_{j=1}^W[1(1q)\epsilon _j/T^{}]^{\frac{q}{1q}}}$$
(3)   
and the additional rule that $`A_i=0`$ whenever $`1(1q)\epsilon _i/T^{}<0`$. Here, to simplify notation, one introduces the auxiliary parameter

$$T^{}=(1q)\underset{i=1}{\overset{W}{}}\epsilon _iP_i+T(\underset{i=1}{\overset{W}{}}P_i^{1/q})^q$$
(4)   
The probabilities $`P_i`$, the entropy $`S_q`$ and the mean value defining the observable $`𝒪_q`$ depend, besides the temperature $`T`$, on a parameter $`q`$, which measures the degree of non–extensivity of the TS. It is not possible, in general, to solve the previous equations to give explicit expressions for the probabilities $`P_i`$ as a function of $`q`$ and $`T`$. An exception being the limit $`q1`$ in which one recovers the Boltzmann–Gibbs Statistics (BGS): $`P_i\mathrm{exp}(\epsilon _i/T)`$ and $`S_1=_{i=1}^WP_i\mathrm{ln}(P_i)`$. For systems with a small number $`W`$ of configurations it is possible to solve equations (3,4) iteratively starting from an initial ansatz for $`P_i`$, e.g. the Boltzmann-Gibbs expression. However, this method is not useful for a system with a moderately large number of configurations.

The fact that one can not give explicit expressions for the probabilities $`P_i`$ has hampered the development of numerical methods to perform the usual Metropolis Monte Carlo or Molecular Dynamics simulations of the TS for interacting systems. Very recently, however, methods based on the numerical calculation of the number of configurations with a given energy have allowed the direct calculation of the necessary averages. In this paper, we introduce a new and more direct method, based on the standard Metropolis algorithm combined with a numerical integration, which can be used in many cases to perform the thermodynamic averages involved in the TS. In the next section we describe the method in some detail, and in section 3 we show the results of the application of the method to the Ising model as well as some of the difficulties encountered.

## 2 The Monte Carlo method

As mentioned before, the main problem to perform a numerical simulation of a system described by the TS at fixed temperature $`T`$ is that we do not have at hand the solution for the probabilities $`P_i`$, since the nonlinear equations (3,4) have no explicit solution for $`q1`$. For $`q=1`$ (BGS), it is $`P_i=𝒵^1\mathrm{exp}(\epsilon _i/T)`$, and one can use a variety of Monte Carlo techniques for the numerical calculation of the averages, Eq.(1). For example: in the Metropolis Monte Carlo algorithm, one generates a change in the configuration $`ij`$ and the new configuration $`j`$ is accepted with a probability $`\mathrm{min}(1,P_j/P_i)=\mathrm{min}[1,\mathrm{exp}((\epsilon _j\epsilon _i)/T)]`$. Notice that the partition function $`𝒵`$ cancels out in the calculation of the acceptance probabilities. Unfortunately, since for $`q1`$ the probabilities $`P_i`$ are not known as a function of $`T`$, there is no trivial generalization of the Monte Carlo method to perform the averages (1) at fixed temperature $`T`$. The method we propose in this paper works in two steps: (i) we perform Monte Carlo simulations at a fixed value of the auxiliary parameter $`T^{}`$; (ii) we then use Eq.(4) in order to determine the physical temperature $`T`$. We describe now in detail both steps.

(i) To perform a Metropolis Monte Carlo simulation of the TS at a fixed “fictitious temperature” $`T^{}`$, one proposes a change in the configuration $`ij`$ and accepts this change with probability $`\mathrm{min}(1,P_j/P_i)`$. Using Eq.(3), one notices that the normalizing factor cancels out:

$$P_{accep}=\mathrm{min}[1,\left[\frac{T^{}(1q)\epsilon _j}{T^{}(1q)\epsilon _i}\right]^{\frac{q}{1q}}]$$
(5)   
(it is also understood that the acceptance probability is zero if the configuration $`j`$ is such that $`T^{}(1q)\epsilon _j<0`$). By using this Monte Carlo algorithm, one generates a sequence of $`M`$ representative configurations which are distributed according to the probability Eq.(3). The statistical averages, Eq.(1), are then approximated by sample averages $`𝒪_q=_{k=1}^MO_k/M`$ and the errors computed in the standard way.

(ii) To perform the $`T^{}T`$ transformation, we invert Eq.(4) using Eqs. (1,2):

$$T=\frac{T^{}(1q)U_q(T^{})}{1+(1q)S_q(T^{})}$$
(6)   
In this expression, the energy $`U_q(T^{})`$ is obtained in the Monte Carlo simulation which is performed at fixed $`T^{}`$. In order to compute the entropy $`S_q(T^{})`$ we make use of the thermodynamic relation $`1/T=S_q/U_q`$ which, using Eq.(6), can be integrated between two equilibrium states characterized by values $`T_0^{}`$ and $`T^{}`$ of the parameter, to yield:

$$\frac{1}{1q}\mathrm{ln}\left[\frac{1+(1q)S_q(T^{})}{1+(1q)S_q(T_0^{})}\right]=\underset{T_0^{}}{\overset{T^{}}{}}\frac{dU_q}{T^{}(1q)U_q}$$
(7)   
The temperature $`T`$ is finally given by:

$$T=\frac{T^{}(1q)U_q(T^{})}{1+(1q)S_q(T_0^{})}\mathrm{exp}\left[(q1)\underset{T_0^{}}{\overset{T^{}}{}}\frac{dU_q}{T^{}(1q)U_q}\right]$$
(8)   
In summary, one performs a Monte Carlo simulation using the Metropolis acceptance probability corresponding to states with a fixed value of $`T^{}`$, starting from some initial value $`T_0^{}`$ up to the desired valued of $`T^{}`$. In these simulations one computes, using the standard Monte Carlo procedure, the sample value of the energy $`U_q(T^{})`$. Finally, in order to obtain the physical temperature $`T`$ one uses the previous expression Eq.(8) where the integral is computed numerically. The initial value $`T_0^{}`$ must be some limiting value in which the entropy $`S_q(T_0^{})`$ is known. This depends on the problem, although obvious candidates are the very high or very low temperature configurations. One could think that a disadvantage of the method is that the integration in Eq.(8) requires a large number of simulations at different fictitious temperatures to be able to perform the $`T^{}T`$ transformation accurately. However, we will show that the use of the reweighting techniques reduces drastically the number of simulations.

## 3 The $`2`$-dimensional Ising Model

To illustrate the method described in the previous section we present results for the short range Ising Model, defined by the following Hamiltonian:

$$=\underset{<i,j>}{}(1S_iS_j)$$
(9)   
where each of the $`N=L\times L`$ spin variables $`S_i`$ can take the values $`\pm 1`$, and the sum $`_{<i,j>}`$ runs over all nearest neighbor sites on a $`2`$-dimensional lattice with periodic boundary conditions.

The Metropolis Monte Carlo simulation is performed, as usual, by randomly choosing one of the spins, $`S_i`$, and proposing a change in configuration in which the spin $`S_i`$ flips its value, $`S_iS_i`$. This change is accepted with a probability given by Eq.(5) whereas, if rejected, the spin keeps its old value.

In the Fig.(1) we show the raw Monte Carlo data for $`L=4,10,20,30`$, in the case $`q=0.8`$. In order to produce this plot, we have used the Metropolis algorithm at the points marked by symbols in the plot. The lines joining the points in the Fig. (1) have been obtained by a reweighting technique, which allows the calculation of mean values of observables for a $`T^{}`$ different from the actual $`T_a^{}`$ where the simulation was performed. In the case of the TS and the probabilities given by Eq.(2), the reweighting is based upon the exact relation:

$$𝒪_q(T^{})=\frac{_\epsilon O(\epsilon )H(\epsilon ;T_a^{})[\frac{1(1q)\epsilon /T^{}}{1(1q)\epsilon /T_a^{}}]^{\frac{q}{1q}}}{_\epsilon H(\epsilon ;T_a^{})[\frac{1(1q)\epsilon /T^{}}{1(1q)\epsilon /T_a^{}}]^{\frac{q}{1q}}}$$
(10)   
where $`H(\epsilon ;T_a^{})`$ is the histogram of all the energy values generated in the Monte Carlo run at $`T_a^{}`$. Here $`O(\epsilon )`$ is the microcanonical mean value at energy $`\epsilon `$ of the observable $`O`$. Notice that in Eq.(10) the sums run over the energy levels $`\epsilon `$ at variance with Eq.(1) where the sums run over all configurations $`i`$.

Once we have obtained the averages as a function of $`T^{}`$ as shown in the figure (1) for the energy, the value of the temperature $`T`$ is obtained from Eq.(8) by performing the indicated integration from the chosen temperature $`T_0^{}`$. In the Ising model, we have used the limit $`T_0^{}=0`$ for which the representative configurations are the two ground states, $`U_q=0`$, and the entropy is: $`S_q(T_0^{}=0)=\frac{2^{(1q)}}{1q}`$.

In Fig.(2) we plot the function that has to be integrated in order to perform the $`T^{}T`$ transformation. The main goal of this figure is to show that the function to integrate is, at least for these values of the parameters and system sizes, a smooth function. For the actual integration, we have used Simpson’s $`3/8`$ rule. Finally, in Fig.(3) we plot the resulting internal energy $`U_q`$ as a function of the actual temperature $`T`$, for several values of the parameter $`q`$ and different system sizes. We note that for $`q<1`$ one obtains a hysteresis–like loop that induces an ambiguity in the actual value of the energy. This is a generic feature of the TS, which can be resolved by applying a minimum free energy criterion to choose the most stable solution. We have observed that for $`q>1`$ and large system sizes (this is not the case in Fig. 3) the Monte Carlo simulations in some temperature ranges near the ferromagnetic transition take a very long time to equilibrate, thus preventing us from performing a very accurate measurement. We believe this is a generic feature of any dynamical updating scheme one could use to simulate the TS.

Finally, we want to remark that Eq.(8) can be combined with any other simulation technique which performs a sampling of the configuration space according to the probability Eq.(3). For instance, one could use the Molecular Dynamic methods using thermostats at $`T^{}`$  to study the dynamic behavior of systems with a large number of degrees of freedom.

Acknowledgments

We acknowledge financial support from DGES, grants PB94-1167 and PB97-0141-C02-01.