# Broad Histogram Method for Continuous Systems : The XY-Model

## 1 Introduction

For thermodynamic systems in equilibrium one typically wants to calculate average values of quantities $`Q`$ like the energy, the magnetization, etc.. These mean values depend on one hand of the properties of the system himself and on the other hand on its interactions with the environment. For example, in the case of the canonical ensemble, the system is in equilibrium with a heat bath at temperature $`T`$, and the mean value $`<Q>_T`$ is given by

$$<Q>_T=\frac{\underset{E}{}<Q>_Eg(E)\mathrm{exp}(E/k_BT)}{_Eg(E)\mathrm{exp}(E/k_BT)},$$
(1)   
where the degeneracy $`g(E)`$ is the number of distinct states with energy $`E`$, $`<Q>_E`$ denotes the micro-canonical average of $`Q`$ at energy $`E`$ and $`k_B`$ is the Boltzmann constant (in the rest of this paper $`k_B=1`$ is taken). Here, $`g(E)`$ and $`<Q>_E`$ are intrinsic characteristics of the system, and $`\mathrm{exp}(E/k_BT)`$ is the Boltzmann factor.

Most Monte Carlo methods were designed to calculate particular averages like Eq. (1). For example, Metropolis algorithms  and similar canonical simulations accumulate sampling distributions that approximate Eq. (1) for a large number of samples and, therefore, the mean values $`<Q>_T`$ are estimated simply by averaging $`Q`$ on the sample set. However, if one wants to study the dependence of $`<Q>_T`$ on $`T`$ using this kind of algorithms it is necessary to run one entire simulation at each temperature value, with large computer efforts.

A well-proved strategy to avoid these multiple calculations is the histogram method (HMC), introduced by Salsburg  and popularized by Ferrenberg and Swendsen . This method reweights data of one canonical simulation at one temperature $`T_\mathrm{o}`$ to find mean values at other temperatures. However, a canonical distribution has a very narrow peak on the energy axis around its mean energy and, therefore, there are not enough samples in the tails of the distribution to have good statistics there. This restricts the reliability of the method to a relatively narrow temperature range around $`T_\mathrm{o}`$ .

The Broad Histogram Monte Carlo Method (BHMC) - was established by de Oliveira et. al. in order to overcome these restrictions. It is designed to determine the system density of states $`g(E)`$ and $`<Q>_E`$ directly. These quantities are calculated from micro-canonical averages and, therefore, samples with different energies can be considered independently. The samples can be obtained in many different ways, for instance, performing micro-canonical simulations  or performing a non-biased random walk along the energy axis, which gives rise to much broader histograms than the usual histogram method . Given $`g(E)`$ and $`<Q>_E`$, the canonical distribution and the desired mean values for any temperature can be obtained. Thus, it is possible to calculate mean values over the whole temperature range of interest in only one run. The BHMC method was first developed for systems with discrete degrees of freedom, and has shown to give results more accurately and using less computer time than the conventional histogram method for the 2d an 3d Ising model -, the Edwards Anderson spin glass  and the 3d Potts glass .

In the present paper we propose a way of extending the method to systems with continuous degrees of freedom, using the XY-model as testing ground. Our way of treatment is general, it can also be applied to other systems of this kind. The remaining part of the paper is organized as follows. In Sec. 2 we briefly summarize the basis of the BHMC method. In Sec. 3 we describe the XY-model and we use it as example to show how to apply the method to continuous systems. In Sec. 4 we discuss several ways of taking the samples. In Sec. 5 we give technical details about the implementation we did of the method. In Sec. 6 we present the results obtained with the BHMC method from a 3D XY-model on a $`10\times 10\times 10`$ cubic lattice with periodic boundary conditions, and we compare these results with those of Metropolis simulations and of the histogram method. On one side, we have obtained an excellent agreement between the BHMC method and Metropolis simulations in all the calculated quantities over a very broad temperature range, i.e. for temperatures between $`0.7`$ and $`4.7`$, using less computer time for the BHMC method than the time required for $`16`$ points by Metropolis simulations. On the other side, it reproduces with the same accuracy the results of the HMC method on the temperature range $`2.1<T<2.2`$ but, instead of diverging for temperatures out of this range as the HMC method does, it remains precise over the whole range $`1.2<T<4.7`$ using only $`2.2`$ times more computer effort than the HMC method. It justifies the name “Broad Histogram” for the method. Finally, in Sec. 7 we summarize the main steps of the proposed strategy to apply the BHMC method to continuous systems and we discuss other strategies that can be used to implement the method.

## 2 The BHMC Method

The strategy of the BHMC method is to calculate the degeneracy of energy states $`g(E)`$ directly. The method for systems with discrete degrees of freedom can be summarized as follows.

First, let us imagine a protocol of allowed movements in the space of states of the system such that changing from an $`X_{old}`$ to an $`X_{new}`$ configuration is allowed if and only if the reverse change is also allowed, i.e.

$`X_{old}X_{new}X_{new}X_{old}`$ (2)   
is allowedis allowed   
In other words, the protocol is micro-reversible. These movements are only virtual, in the sense that they are not performed. They are introduced only to estimate the density of states g(E).

Next, let us choose a fixed amount of energy change $`\mathrm{\Delta }E_{fix}`$ and compute $`N_{up}(X)`$ ($`N_{dn}(X)`$) for the configuration X as the number of allowed changes that increases (decreases) the energy of the configuration by $`\mathrm{\Delta }E_{fix}`$. Let $`<N_{up}(E)>`$ ($`<N_{dn}(E)>`$) be the micro-canonical average of $`N_{up}(X)`$ ($`N_{dn}(X)`$) at energy $`E`$.

Due to Eq. (2), the total number of ways to go up, summed over all states with energy $`E`$, is equal to the total number of ways to go down, summed over all states with energy $`E+\mathrm{\Delta }E_{fix}`$, i.e.

$$g(E)<N_{up}(E)>=g(E+\mathrm{\Delta }E_{fix})<N_{dn}(E+\mathrm{\Delta }E_{fix})>.$$
(3)   
Eq. (3) can be rewritten as:

$$\mathrm{ln}g(E+\mathrm{\Delta }E_{fix})\mathrm{ln}g(E)=\mathrm{ln}\frac{<N_{up}(E)>}{<N_{dn}(E+\mathrm{\Delta }E_{fix})>}.$$
(4)   
Hence, knowing $`<N_{up}(E)>`$ and $`<N_{dn}(E)>`$ allows to calculate the right hand side of Eq. (4) and obtains $`\mathrm{ln}g(E)`$ for all values of $`E`$ in steps of $`\mathrm{\Delta }E_{fix}`$ by adding up differences from a given initial point $`(E_\mathrm{o},\mathrm{ln}g(E_\mathrm{o}))`$.

Similarly, it can be observed that Eq. (4) divided by $`\mathrm{\Delta }E_{fix}`$ estimates $`\beta (E)`$, i.e.

$`\beta (E){\displaystyle \frac{d\mathrm{ln}g(E)}{dE}}`$ $``$ $`{\displaystyle \frac{\mathrm{ln}g(E+\mathrm{\Delta }E_{fix})\mathrm{ln}g(E)}{\mathrm{\Delta }E_{fix}}}`$ (5)   
$``$ $`{\displaystyle \frac{1}{\mathrm{\Delta }E_{fix}}}\mathrm{ln}{\displaystyle \frac{<N_{up}(E)>}{<N_{dn}(E+\mathrm{\Delta }E_{fix})>}}.`$   
Then, $`\beta (E)`$ can be obtained for a set of energy values $`\{E_1,E_2,\mathrm{},E_n\}`$ by estimating $`<N_{up}(E_i)>`$ and $`<N_{dn}(E+\mathrm{\Delta }E_{i}^{}{}_{fix}{}^{})>`$ for each value $`E_i`$. Eq. (4) can be used to estimate $`\beta (E)`$ in more accurate ways, as it will be discussed in Sec. 4. $`\mathrm{ln}g(E)`$ can be obtained, thereafter, by integrating $`\beta (E)`$ numerically from a given initial point $`(E_\mathrm{o},\mathrm{ln}g(E_\mathrm{o}))`$.

The value of $`\mathrm{ln}g(E_\mathrm{o})`$ for the initial point does not affect the calculation of mean values $`<Q>_T`$, as can be shown in Eq. (1). If $`\mathrm{ln}g(E_\mathrm{o})=0`$ is taken, the method gives $`g(E)/g(E_\mathrm{o})`$ for each energy.

Finally, to carry out calculations with this method four histograms are required, i.e. $`N_{up}(E)`$, $`N_{dn}(E)`$, $`Q(E)`$ and the number of visits $`V(E)`$. Estimations for $`<N_{up}(E)>`$ and $`<N_{dn}(E)>`$ and, in consequence, an estimation of $`g(E)`$ are obtained by dividing $`N_{up}(E)`$ and $`N_{dn}(E)`$ by $`V(E)`$. $`<Q>_E`$ is calculated by dividing $`Q(E)`$ by $`V(E)`$. Combining these two results Eq. (1) can be used to calculate $`<Q>_T`$ at any desired temperature.

It can be shown that Eq. (3) is exact for any statistical model . Monte Carlo sampling processes are used only to estimate the micro-canonical averages $`<N_{up}(E)>`$ and $`<N_{dn}(E)>`$. Many different sampling strategies can be employed, as will be discussed in Secs. 4 and 7. Nevertheless, all states with the same energy have to be sampled with the same probability and the correlations between samples have to be small enough to obtain good estimations of such micro-canonical averages. Otherwise, problems can appear .

## 3 The BHMC method and the XY-Model

### 3.1 The XY-model

In order to describe our technique for extending the BHMC method to systems with continuous degrees of freedom we have chosen the XY-model as testing ground. The XY-model -, consists of a set of spins $`\stackrel{}{\sigma }`$ of length unity arranged on a lattice. Each spin is allowed to rotate in a plane, characterized by an angle $`\theta `$. These variables define the state of the system and can change continuously in the range $`\theta [\pi ,\pi ]`$. The Hamiltonian with zero external field is given by

$$=J\underset{<\mathrm{𝑖𝑗}>}{}\stackrel{}{\sigma }_i\stackrel{}{\sigma }_j=J\underset{<\mathrm{𝑖𝑗}>}{}\mathrm{cos}(\theta _i\theta _j)$$
(6)   
where the summations $`<ij>`$ are taken over all pairs of nearest-neighbor sites, $`\stackrel{}{\sigma }_i\stackrel{}{\sigma }_j`$ is the scalar product between $`\stackrel{}{\sigma }_i`$ and $`\stackrel{}{\sigma }_j`$, and $`J`$ denotes the maximal energy per bond (in the rest of this paper we will take $`J=1`$). The order parameter chosen for this model  is the average magnetization $`M`$, defined by:

$$M\underset{i}{}A_i,$$
(7)   
with

$$A_i\sqrt{\left(\underset{j}{}\mathrm{cos}\theta _j\right)^2+\left(\underset{j}{}\mathrm{sin}\theta _j\right)^2}$$
(8)   
where the summations $`j`$ are taken over all nearest neighbors of the site $`i`$. The XY-model so defined is one of the simplest thermodynamic systems with continuous degrees of freedom.

### 3.2 The BHMC method for continuous systems

In order to extend the BHMC method for continuous systems, it is necessary to redefine the quantities involved in Eq. (3) and to find a condition for the protocol of movements equivalent to Eq. (2) such that Eq. (3) holds.

First, let us imagine a protocol of random movements in the space of states of the system, such that for each allowed movement the probability to perform it equals the probability to revert it. These movements are again only virtual, in the sense that they are not executed. They are introduced only to estimate the density of states g(E). For the XY-model, for example, we use the following protocol: we choose one spin at random with angle $`\theta _{old}`$, and then we choose for it a new angle $`\theta _{new}`$ according to a uniform probability distribution over the range $`\theta _{new}[\pi ,+\pi ]`$. It has to be noted that $`\theta _{new}`$ is a random variable, with a probability density function (p.d.f) $`f_{\theta _{new}}(\theta )`$ defined as the probability of obtaining a value $`\theta _{new}`$ between $`\theta `$ and $`\theta +d\theta `$ given by

$$f_{\theta _{new}}(\theta )=\{\begin{array}{ccc}1/2\pi & :& \pi <\theta _{new}<\pi \\ 0& :& otherwise\end{array}.$$
(9)   
Thus, the protocol of random movements can be defined by giving the p.d.f.’s for the new values of the system variables. This seems a natural implementation for continuous systems.

The next step consist in estimating $`N_{up}`$ and $`N_{dn}`$ for a given configuration with energy $`E_{old}`$. The problem arises because there is a continuum non-numerable set of configurations that can be reached. However, it has to be noted that each one of these new configurations has a well-defined energy value $`E_{new}`$ and that exactly one of them is chosen at random. Therefore, we can define $`E_{new}`$ and the energy change $`\mathrm{\Delta }EE_{new}E_{old}`$ as random variables. The p.d.f. $`f_{\mathrm{\Delta }E}(\mathrm{\Delta }E)`$ is defined by the probability of obtaining an energy change between $`\mathrm{\Delta }E`$ and $`\mathrm{\Delta }E+d\mathrm{\Delta }E`$. We propose to redefine $`N_{up}`$ and $`N_{dn}`$ as

$`N_{up}f_{\mathrm{\Delta }E}(\mathrm{\Delta }E_{fix})`$ ; $`N_{dn}f_{\mathrm{\Delta }E}(\mathrm{\Delta }E_{fix})`$ (10)   
From this point, the method proceeds the same as in the discrete case. Eqs. (10) extend the BHMC method to systems with continuous degrees of freedom. Therefore, the problem reduces to finding the function $`f_{\mathrm{\Delta }E}`$ for a given configuration $`X_{old}`$.

### 3.3 Finding $`f_{\mathrm{\Delta }E}`$

The procedure of finding $`f_{\mathrm{\Delta }E}`$ is straightforward. First, the protocol of random movements is defined by giving the p.d.f.’s for the new values of the system variables, as in Sec. 3.2. Then, the energy $`E_{new}`$, is expressed in terms of such new values. Next, its p.d.f. $`f_{E_{new}}`$ is determined by using the usual rules of finding the p.d.f. for a function of random variables . Finally $`f_{\mathrm{\Delta }E}(\mathrm{\Delta }E)`$ is found, according to these rules, replacing $`E`$ in $`f_{E_{new}}(E)`$ by $`\mathrm{\Delta }E+E_{old}`$.

To continue with our example, we analyze first what happens with the turn of one specific spin. Suppose that we choose to turn spin $`i`$. Its bonds are the only ones that change the energy value. We define $`\epsilon _i`$ as sum of the energies of these bonds, i.e.

$$\epsilon _i\underset{j}{}\mathrm{cos}(\theta _i\theta _j).$$
(11)   
where the summations are taken over all nearest neighbors $`j`$ of the site $`i`$. When $`\theta _i`$ changes from its previous value $`\theta _{i}^{}{}_{old}{}^{}`$ to the new value $`\theta _{i}^{}{}_{new}{}^{}`$, $`\epsilon _i`$ changes from $`\epsilon _{i}^{}{}_{old}{}^{}`$ to $`\epsilon _{i}^{}{}_{new}{}^{}`$. It is clear that $`\mathrm{\Delta }E=\epsilon _{i}^{}{}_{new}{}^{}\epsilon _{i}^{}{}_{old}{}^{}`$ and that $`\epsilon _{i}^{}{}_{new}{}^{}`$ is a function of the random variable $`\theta _{i}^{}{}_{new}{}^{}`$. Therefore, $`\epsilon _{i}^{}{}_{old}{}^{}`$ and $`\epsilon _{i}^{}{}_{new}{}^{}`$ can be used instead of $`E_{new}`$ and $`E_{old}`$ for the general treatment described above.

In order to express $`\epsilon _{i}^{}{}_{new}{}^{}`$ in terms of $`\theta _{i}^{}{}_{new}{}^{}`$, Eq. (11) can be rewritten as

$$\epsilon _{i}^{}{}_{new}{}^{}=A_i\mathrm{cos}(\theta _{i}^{}{}_{new}{}^{}\delta \theta _i),$$
(12)   
where

$$\delta \theta _i=\mathrm{arctan}\frac{\underset{j}{}\mathrm{sin}\theta _j}{_j\mathrm{cos}\theta _j}+\pi $$
(13)   
and $`A_i`$ is defined by Eq. (8).

If $`\theta _{i}^{}{}_{new}{}^{}`$ is a random variable with p.d.f. given by Eq. (9), the usual rules of finding the p.d.f. for a function of random variables can be used to show that $`\epsilon _{i}^{}{}_{new}{}^{}`$ is a random variable with p.d.f. given by 

$$f_{\epsilon _{i}^{}{}_{new}{}^{}}(\epsilon )=\{\begin{array}{ccc}\frac{1}{\pi }\frac{1}{\sqrt{A_i^2\epsilon ^2}}& :& |\epsilon _i|<A_i\\ 0& :& otherwise\end{array}.$$
(14)   
Finally, we replace $`\epsilon `$ in Eq. (14) by $`\mathrm{\Delta }E+\epsilon _{i}^{}{}_{old}{}^{}`$, to obtain the p.d.f. $`f_{\mathrm{\Delta }E}^i(\mathrm{\Delta }E)`$ of $`\mathrm{\Delta }E`$ corresponding to a turn of the spin $`i`$, i. e.

$$f_{\mathrm{\Delta }E}^i(\mathrm{\Delta }E)=\{\begin{array}{ccc}\frac{1}{\pi }\frac{1}{\sqrt{A_i^2(\epsilon _{i}^{}{}_{old}{}^{}+\mathrm{\Delta }E)^2}}& :& |\mathrm{\Delta }E+\epsilon _{i}^{}{}_{old}{}^{}|<A_i\\ 0& :& otherwise\end{array}.$$
(15)   
Since all spins are equally probable, the probability of obtaining an energy change $`\mathrm{\Delta }E`$ between $`\mathrm{\Delta }E`$ and $`\mathrm{\Delta }E+d\mathrm{\Delta }E`$ is equal to the sum of $`f_{\mathrm{\Delta }E}^i(\mathrm{\Delta }E)`$ over all spins divided by $`N`$, the number of spins, i.e.

$$f_{\mathrm{\Delta }E}(\mathrm{\Delta }E)=\frac{1}{N}\underset{i}{}f_{\mathrm{\Delta }E}^i(\mathrm{\Delta }E).$$
(16)   
$`N_{up}`$ and $`N_{dn}`$ for a given configuration are obtained finally using Eqs. (10) as

$`N_{up}={\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}N_{up}^i`$ ; $`N_{dn}={\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}N_{dn}^i`$ (17)   
where we have defined

$`N_{up}^if_{\mathrm{\Delta }E}^i(\mathrm{\Delta }E_{fix})`$ ; $`N_{dn}^if_{\mathrm{\Delta }E}^i(\mathrm{\Delta }E_{fix})`$ (18)   
When Eqs. (17) are replaced in Eq.(4), the factor $`1/N`$ in Eq. (16) as well as the factor $`1/\pi `$ in Eq. (15) will eventually cancel out.

The p.d.f. $`f_{\mathrm{\Delta }E}^i(\mathrm{\Delta }E)`$ is shown in Fig. 1. The function diverges at points $`\mathrm{\Delta }E=\pm A_i\epsilon _{i}^{}{}_{old}{}^{}`$ and takes its minimum value between these two points at $`\mathrm{\Delta }E=\epsilon _{old}`$ with $`f_{\mathrm{\Delta }E}^i(\epsilon _{i}^{}{}_{old}{}^{})=(\pi A_i)^1`$. Outside this range the function equals zero. The values $`N_{up}^i`$ and $`N_{dn}^i`$ are also shown.

Summarizing, Eqs. (17) give us the way to calculate $`N_{up}`$ and $`N_{dn}`$ for the XY-model and, therefore, a way to apply the BHMC method to it. The same strategy can be used to apply the method to other continuous systems.

## 4 How to take the samples

A second idea behind the BHMC method is to take the samples in such a way that they distribute almost homogeneous over the whole energy range of interest. There are many different strategies that can be used for this purpose. We have used two of them, namely, a micro-canonical sampling process that maintains the system inside a narrow energy window (BHMC-$`\mu `$C) , and a random walk on the energy axis using Metropolis  steps (BHMC-M) .

The BHMC-$`\mu `$C sampling strategy consist in taking a small window centered at each energy of interest, starting from a configuration inside the window and performing random movements never surpassing the limits of the window. It can be implemented, for instance, by taking at random a new configuration and accepting it only if it falls inside the window, or by using the Creutz energy bag method . According to Eq. (5), two micro-canonical simulations are required to calculate $`\beta (E)`$ at one energy value. However, the data of one simulation can be used for two energy values if their difference equals $`\mathrm{\Delta }E_{fix}`$. For the XY-model, we take one spin and generate at random a new angle $`\theta [\pi ,\pi ]`$. If the total energy of the system falls inside the window, the change is accepted, otherwise it is rejected. It is clear that this sampling process maintains a detailed balance condition with equal probabilities for the old an new configurations and, therefore, it samples all states inside the window with the same probability. A whole lattice sweep is performed by repeating this procedure for all spins of the system, and a new configuration is taken after every fixed number of lattice sweeps.

The BHMC-M sampling strategy consist first in performing a Markovian process with symmetric probability distributions for increasing and decreasing the energy, i.e. a non biased random walk on the energy axis. Second, since the canonical distribution is almost symmetric on the energy axis around its mean energy value $`<E>_T`$, one way to do it (but not the only one) consists of taking the last configuration of the sample and performing on it a fixed number of Metropolis  steps at such a temperature that its canonical distribution appears centered on the energy value $`E`$. The resulting configuration is taken as the next sample. In order to find this temperature $`T(E)`$ Eq. (5) can be used by remembering that $`\beta (E)1/T(E)`$. Furthermore, Eq. (4) can be used to estimate $`\beta (E)`$ and $`T(E)`$ in more accurate ways. In the present paper we use the centered difference estimator

$`\beta (E)`$ $``$ $`{\displaystyle \frac{\mathrm{ln}g(E+\mathrm{\Delta }E_{fix})\mathrm{ln}g(E\mathrm{\Delta }E_{fix})}{2\mathrm{\Delta }E_{fix}}}`$ (19)   
$``$ $`{\displaystyle \frac{1}{2\mathrm{\Delta }E_{fix}}}\mathrm{ln}{\displaystyle \frac{<N_{up}(E)><N_{up}(E\mathrm{\Delta }E_{fix})>}{<N_{dn}(E)><N_{dn}(E+\mathrm{\Delta }E_{fix})>}}.`$   
These values are estimated from the data accumulated in the histograms $`N_{up}(E)`$, $`N_{dn}(E)`$ and $`V(E)`$. It can happens that the system reaches an energy value $`E`$ such that some of the neighboring energies have yet not been sampled and, therefore, $`T(E)`$ cannot be estimated. In such a case, the Metropolis steps are performed at the same temperature used to produce the last configuration.

## 5 Implementation

In order to test our strategy numerically we performed computations for the three-dimensional XY-model on a $`10\times 10\times 10`$ cubic lattice with periodic boundary conditions using the BHMC-M, the BHMC-$`\mu `$C, the HMC and the Metropolis method.

For the Metropolis method the spin system was initialized at each temperature value from a random configuration performing $`N_{ini}=500`$ entire Metropolis lattice sweeps before sampling, and then $`N_s=1500`$ samples were taken separated by $`N_i=10`$ lattice sweeps among them in order to decrease their correlations. Here and in the following the correlation between successive samples was determined by looking at the correlation of their magnetization values. For the Metropolis method we obtained a correlation between successive samples of $`56.5\%`$ at $`T_\mathrm{o}=2.159`$, i.e. at the pseudo critical temperature we have observed for this finite system. Some additional Metropolis simulations were performed as reference data for the comparison between the BHMC-$`\mu `$C and the HMC methods. These simulations were performed in the same way as before, but taking $`N_s=10000`$ samples at each temperature value.

For the BHMC-$`\mu `$C method we divided the whole energy range of positive temperatures, i.e. energies per bond between $`1.0`$ and $`0.0`$, in adjacent windows of equal size and chose $`\mathrm{\Delta }E_{fix}`$ to be equal to this bin size. So $`\mathrm{\Delta }E_{fix}`$ is equal to the difference of the energy values at which consecutive micro-canonical simulations are performed and, therefore, the data of each simulation can be used twice (see Sec. 4). To acquire the samples the spin system was initialized for the lowest-energy window from a random configuration by turning some spins at random until the window is reached and performing $`N_{ini}`$ micro-canonical lattice sweeps before sampling. Then, we took $`N_s`$ samples separated by $`N_i`$ lattice sweep among them. After that, we move to the next window by turning some spins at random starting from the previous sample. Then, $`N_{ini}`$ micro-canonical lattice sweeps were performed before sampling, and so on.

To compare with the Metropolis simulations, we divided the whole energy range of positive temperatures in $`1225`$ adjacent windows, i.e. $`\mathrm{\Delta }E_{fix}2.45`$. To handle the critical slowing down in the relaxation time for the micro-canonical simulations, larger values of $`N_{ini}`$ have been taken in the vicinity of the critical point, i.e. $`N_{ini}=250`$ for energies per bond between $`0.400`$ and $`0.320`$ and $`N_{ini}=70`$ for energies outside this range. To obtain the same correlation between samples as from Metropolis simulations we used $`N_i=3`$ in the energy range $`(0.429,0.282)`$ and $`N_i=2`$ outside this range. The maximal correlation in the vicinity of the critical point so obtained was $`56.4\%`$. We took $`N_s=50`$ for the whole energy range.

To compare with the HMC method, the whole energy range of positive temperatures was divided in $`500`$ adjacent windows, i.e. $`\mathrm{\Delta }E_{fix}=6.0`$. To handle the critical slowing down, we took $`N_{ini}=200`$ for energies per bond between $`0.400`$ and $`0.320`$, and $`N_{ini}=80`$ for energies outside this range. $`N_i=2`$ was enough to obtain lower correlations in the vicinity of the critical point than from Metropolis simulations. The maximal correlation so obtained was $`46.3\%`$. In order to obtain the same accuracy than the HMC method on the temperature range $`2.1<T<2.2`$ more samples have been taken in the vicinity of the critical point, i.e. $`N_s=380`$ for energies per bond between $`0.504`$ and $`0.248`$ and $`N_s=80`$ for energies outside this range.

For the BHMC-M method the histograms were constructed by dividing the whole energy range of positive temperatures in $`1225`$ boxes of equal size. Ten independent systems that accumulated data in common histograms were simulated simultaneously, in order to estimate $`T(E)`$ more accurately. These systems were initialized from random configurations thermalized by $`N_{ini}=500`$ Metropolis lattice sweeps at the critical temperature $`T_\mathrm{c}=2.20196`$ . We took a total number of $`N_s=10\times 10000`$ samples performing $`N_i=1`$ lattice sweeps between them. In order to avoid wasting samples at very low temperatures, we restricted our systems to energies between $`0.90`$ and $`0.0`$ per bond; spin systems going out of this range were reinitialized by copying the values of an additional spin system, which acts as a replacement. This additional system was initialized in the same way as the other ones, and then twenty entire Metropolis lattice sweeps at $`T_\mathrm{c}`$ were applied on it after each time it was taken as replacement.

The HMC method was implemented by dividing the whole energy range of positive temperatures in $`1225`$ and $`500`$ boxes of equal size, and accumulating the data of a Metropolis simulation at $`T_\mathrm{o}=2.159`$ performed with $`N_{ini}=500`$,$`N_i=10`$ and $`N_s=10000`$.

For all methods, we stored $`\mathrm{cos}(\theta _i)`$ and $`\mathrm{sin}(\theta _i)`$ in tables instead of the angles $`\theta _i`$ in order to accelerate the calculations. For the BHMC-$`\mu `$C and BHMC-M methods, we stored in addition $`A_i^2`$ and $`\epsilon _i`$ for each spin for the same reason. The functions $`\mathrm{cos}(x)`$, $`\mathrm{sin}(x)`$, $`\sqrt{x}`$, $`1/\sqrt{x}`$, $`\mathrm{exp}(x)`$, were stored in tables for all methods.

The mean values and error bars for each method were made by repeating eight times the whole simulation with different initial seeds of the random number generator. We used the “minimal” Park and Miller random number generator combined with a Marsaglia shift sequence from Ref. .

## 6 Results

The results are shown in Fig. 2-10. Figs. 2-8 display data obtained by using $`1225`$ divisions of the energy range of positive temperatures, i.e. $`\mathrm{\Delta }E_{fix}2.45`$ for the BHMC-$`\mu `$C, the BHMC-M and the HMC methods.

Fig. 2 shows how samples distribute on the energy axis using the BHMC-$`\mu `$C, the BHMC-M and the HMC methods. For the BHMC methods a large energy range is sampled, in contrast to the narrow window sampled by the classical canonical histogram method. In addition, with the BHMC-$`\mu `$C it is possible to take exactly the same number of samples at each energy value.

Fig. 4 shows $`\beta (E)`$ obtained from the BHMC-$`\mu `$C and the BHMC-M methods. It has to be noted that both curves are almost indistinguishable, in spite of the completely different sampling strategies of both methods.

Fig. 4 shows $`\mathrm{ln}g(E)`$ obtained from the BHMC-$`\mu `$C and the BHMC-M methods. Since the curve from the BHMC-M method was obtained adding the differences of Eq. (4) starting from the point $`(E_\mathrm{o}=0.90,\mathrm{ln}g(E_\mathrm{o})=0)`$, its values $`g(E)`$ are the ratio of the probability to obtain at random an energy per site between $`E`$ and $`E+dE`$ to the probability to obtain an energy per site between this value of $`E_\mathrm{o}`$ and $`E_\mathrm{o}+dE`$, when all states are equally probable. The same applies to the curve obtained from the BHMC-$`\mu `$C method, with $`E_\mathrm{o}=0.94`$. It can be noted, that both curves are only shifted in a constant value, namely $`\mathrm{ln}(g(0.90)/g(0.94))`$.

Figs. 6-8 compare results from the BHMC methods with Metropolis simulations. Figs. 6 and 6 show the average energy and the average magnetization obtained with the BHMC-$`\mu `$C and the BHMC-M methods compared with the same quantities obtained using Metropolis simulations. The agreement between the three Monte Carlo methods is excellent. The data obtained by the BHMC methods fit so well with those obtained by Metropolis simulations, that their curves and their error bars are almost indistinguishable.

In Fig. 8 and 8 we show the heat capacity and the magnetic susceptibility obtained by using the BHMC-$`\mu `$C, the BHMC-M and the Metropolis simulations. The curves calculated by the BHMC-M method fit very well within their error bars with them of the Metropolis simulations for the temperature range between $`T=4.7`$ and $`T=0.7`$, and deviate for lower temperatures. The curves calculated by the BHMC-$`\mu `$C method fit as well as those from the BHMC-M method in the same temperature range, begins to oscillate at $`T=0.7`$ and deviates for a lower temperature, namely $`T=0.45`$. This is due to a lack of good statistics at low energies. In the case of BHMC-$`\mu `$C method most of the possible changes at low energy values increase the system energy and, therefore, it is difficult to obtain estimations for $`<N_{dn}(E)>`$ different from zero at this energy values. This difficulty increases when $`\mathrm{\Delta }E_{fix}`$ increases. In the case of the BHMC-M method we have not sampled the energy range below $`E=0.9`$ at all. It can be observed that this fact only affects the results at low temperatures. For all methods, the errors bars become larger around the critical temperature $`T_\mathrm{c}`$, as it would be expected.

For all these curves, we present error bars of $`1.0`$ standard deviations. For eight runs, the confidence level of $`99\%`$ corresponds to $`3.5`$ standard deviations.

Since all sets of data obtained from the BHMC-$`\mu `$C method and the Metropolis simulations present error bars of similar size using the same correlation between samples, we are allowed to compare the speeds of the two algorithms. The Metropolis method took $`28.5`$ seconds on a Digital Alpha workstation to make one simulation at one temperature point. The BHMC-$`\mu `$C method took $`440`$ seconds to perform one calculation on the whole temperature range using the same machine. It means that our method calculates the whole temperature range $`0.7<T<4.7`$ using approximately $`15.5`$ times the computer effort required by a Metropolis simulation for one point. In addition, the BHMC-M method took $`409`$ seconds per run.

Figs. 10 and 10 compare the heat capacity obtained from the BHMC-$`\mu `$C and the HMC methods using Metropolis simulations as reference data. Both figures display data obtained from the BHMC-$`\mu `$C and the HMC methods using $`500`$ divisions of the energy axis, i.e. $`\mathrm{\Delta }E=6.0`$.

Fig. 10 presents the usual error bars of $`1.0`$ standard deviations and Fig. 10 presents error bars of $`3.5`$ standard deviations corresponding to a confidence level of 99% for eight runs, i. e. the probability of finding the true value outside these error bars due only to statistical variations is equal to $`1\%`$. It can be observed that all methods coincide in the temperature range between $`T=2.1`$ and $`T=2.2`$. This is the temperature range of validity of the HMC method according to the criteria that the mean energy value calculated at temperature $`T`$ should be less than $`\sigma _E`$ away of the mean energy value at $`T_\mathrm{o}`$, where $`T_\mathrm{o}`$ is the temperature used to generate the samples and $`\sigma _E`$ is the standard deviation on the energy of these samples . For temperatures out of this range, the HMC method loses its precision but the BHMC-$`\mu `$C method is still precise on the whole temperature range $`1.2<T<4.7`$. Since the error bars in the range $`2.1<T<2.2`$ are similar in size and the correlations between samples are also similar, we can compare the speeds of both methods. The HMC method took $`188`$ seconds per run on a Digital Alpha workstation. The BHMC-$`\mu `$C method took $`411`$ seconds per run using the same machine. It means that our method gives precise results for the temperature range $`1.2<T<4.7`$, i.e. $`35`$ times larger than the range $`2.1<T<2.2`$ of precise results given by the HMC method, using only $`2.2`$ times more CPU time.

## 7 Conclusions

We have proposed a way of applying the BHMC method to continuous systems. We have found results for the three-dimensional SC XY-model in excellent agreement with the HMC method and Metropolis simulations. Our method calculates the temperature range $`0.7<T<4.7`$ using less computer time than the Metropolis simulations for $`16`$ points, and calculates with precision the temperature range $`1.2<T<4.7`$ using only $`2.2`$ times more computer effort than the HMC method for the range $`2.1<T<2.2`$.

The strategy proposed could also be applied to other thermodynamic systems with continuous degrees of freedom, for example, the Heisenberg model and fluid models. Our strategy can be summarized as follows:

* Choose a protocol of random movements in the space of states of the system, such that for each allowed movement the probability to perform it equals the probability to revert it. These movements are only virtual, in the sense that they are not executed. They are introduced only to calculate $`Nup`$, $`Ndn`$ and, therefore, to estimate the density of states g(E). The protocol would change the system from an $`old`$ to a $`new`$ state. Express the protocol giving the probability density function (p.d.f.) for the values of the system variables at the $`new`$ state. These new values are, therefore, random variables.
* Find the p.d.f. $`f_{\mathrm{\Delta }E}(\mathrm{\Delta }E)`$, i.e. the probability of obtaining an energy change between $`\mathrm{\Delta }E`$ and $`\mathrm{\Delta }E+d\mathrm{\Delta }E`$ using the chosen protocol. This can be done as follows.

1. Express the energy of the $`new`$ state $`E_{new}`$ as a function of the new values of the system variables. $`E_{new}`$ is, therefore, a random variable.
    2. Use the usual rules of finding the p.d.f. for a function of random variables to obtain $`f_{E_{new}}(E)`$, i.e. the p.d.f. of $`E_{new}`$.
    3. Replace $`E`$ in $`f_{E_{new}}(E)`$ by $`\mathrm{\Delta }E+E_{old}`$ to obtain, $`f_{\mathrm{\Delta }E}(\mathrm{\Delta }E)`$. Here $`E_{old}`$ denotes the energy of the old state.
* Take $`N_{up}`$ ($`N_{dn}`$) for the $`old`$ state as proportional to $`f_{\mathrm{\Delta }E}(\mathrm{\Delta }E_{fix})`$   
($`f_{\mathrm{\Delta }E}(\mathrm{\Delta }E_{fix})`$), and continue with the BHMC method as usual. $`\mathrm{\Delta }E_{fix}`$ denotes the fixed energy interval of the BHMC method.

The BHMC method calculates the degeneracy $`g(E)`$ either by adding terms of the form $`\mathrm{ln}g(E+\mathrm{\Delta }E_{fix})\mathrm{ln}g(E)`$ or by integrating $`\beta (E)`$, starting in both cases from a given initial point $`(E_\mathrm{o},\mathrm{ln}g(E_\mathrm{o}))`$. If $`\mathrm{ln}g(E_\mathrm{o})=0`$ is taken, $`g(E)`$ can be reinterpreted as the ratio of the probability to obtain at random an energy between $`E`$ and $`E+dE`$ to the probability to obtain an energy between $`E_\mathrm{o}`$ and $`E_\mathrm{o}+dE`$. This interpretation of the values $`g(E)`$ obtained from the BHMC method is valid both for discrete and continuous systems.

Many different strategies can be used to take the samples in the BHMC method. In the present paper we used two of them, namely, a micro-canonical sampling process that maintains the system inside a narrow energy window (BHMC-$`\mu `$C) , and a random walk on the energy axis using Metropolis  steps (BHMC-M) . Between these two methods, we prefer the first one, because of the following reasons. First, it resembles better the micro-canonical character of the averages $`<N_{up}(E)>`$ and $`<N_{dn}(E)>`$. Second, it offers a fine control on the correlations between samples and the relaxation time at each energy value. Third, two micro-canonical simulations are enough to estimate $`\beta (E)`$ at any desired energy and this fact is independent of the system size. Therefore, it seems that the computer time will grow linearly with the system size but, of course, the possibility of using data of one simulation for two energy points discussed in Sec. 4 will be lost. This will be studied in a future paper. Fourth, many other kinds of micro-canonical simulations can be employed to take the samples, different from the accepting-rejecting process used in this paper but maintaining also a detailed balance condition with equal probabilities for all configurations inside a window. For instance, the new configurations can be produced to be always inside the window, and this will decrease relaxation and correlation times . Another option would be to use Creutz’s demons to take the samples.

Other possible sampling strategy is the Entropic Sampling of Ref. . To use it in conjunction with the BHMC method is a very interesting idea, because the BHMC method gives direct estimations of $`\mathrm{ln}g(E+\mathrm{\Delta }E_{fix})\mathrm{ln}g(E)`$ and, therefore, no previous estimation of the entropy is needed to perform the sampling process. Something similar has been employed in the original work of Ref. , but estimating these differences only from the last sample.

In conclusion, the strategy proposed here in order to apply the BHMC method to systems with continuous degrees of freedom gives excellent results for the three dimensional XY-model and seems feasible of being applied to other thermodynamic systems. It has to be noted, that the BHMC method seems a very precise method to find $`\beta (E)`$ and $`\mathrm{ln}g(E)`$ for every thermodynamic system and it seems also extensible to the calculation of $`g(E)`$ for non-thermodynamic systems.

## 8 Acknowledgments

We thank P.M.C. de Oliveira for helpful discussions. One of the authors (Muñoz) thanks Fernec Kun. The work of one of the authors (Muñoz) is supported by the Deutscher Akademischer Austauschdienst (DAAD) through scholarship A/96/0390.