text:
# Measurement of the 𝑝̄𝑝→𝐾_𝑆𝐾_𝑆𝜂 cross section at beam momenta in the regions of 1.45 and 1.7 GeV/c

## 1 INTRODUCTION

The PS185 experiment at the Low Energy Antiproton Ring, LEAR, at CERN was designed to study strangeness production in $`\overline{p}p`$ collisions. Most attention has been drawn to the two-body $`\overline{p}p\overline{\mathrm{\Lambda }}\mathrm{\Lambda }`$ reaction . Several observables, such as total and differential cross sections, polarisations and spin correlations, have been measured over the momentum region from threshold to 2 GeV/c. PS185 has also investigated other two–body final states in the $`\overline{p}p\overline{\mathrm{\Sigma }}^0\mathrm{\Lambda }+c.c.`$ , $`\overline{p}p\overline{\mathrm{\Sigma }^+}\mathrm{\Sigma }^+`$ and $`\overline{p}p\overline{\mathrm{\Sigma }^{}}\mathrm{\Sigma }^{}`$  and $`\overline{p}pK_SK_S`$  reactions.

These studies are now extended to many–body final–states with two neutral strange particles. We present in this paper a measurement of the $`\overline{p}pK_SK_S\eta `$ cross section from $`\overline{p}p`$ annihilation in flight. This constitutes the first cross section measurement of the reaction since the only other available data on this final state is the branching ratio from $`\overline{p}p`$ annihilation at rest, obtained in a bubble chamber experiment .

A study of the $`\overline{p}pK_SK_S\eta `$ reaction is interesting because all three final–state particles contain strangeness and this may provide additional information on strangeness production. In a constituent quark model, the process can be viewed as the annihilation of quark–antiquark pairs with the subsequent creation of strange–antistrange quark pairs. Such models are depicted in Fig. 1a. The process could also proceed via the dissociation of a strangeness component in the (anti-)proton, as shown in Fig. 1b, an idea which has already been applied to the $`\overline{p}p\overline{\mathrm{\Lambda }}\mathrm{\Lambda }`$ reaction . More conventionally, the reaction can be described in a one-baryon-exchange model via charged $`\mathrm{\Sigma }`$ exchange. The final state can then either be directly produced (Fig. 1c) or proceed via the excitation of a resonance (Fig. 1d). In the latter case one would expect the $`K_3^{}(1780)`$ resonance to be important in flight since it is the only known $`K^{}`$ resonance in this mass region with a substantial decay branching ratio into the $`K_S\eta `$ channel $`(30\pm 13)\%`$ . However, to the best of our knowledge, no theoretical calculation has been made for this reaction.

## 2 EXPERIMENT AND ANALYSIS

In this study, we used data from two PS185 runs. The first involved six beam momentum settings between 1.440 and 1.477 GeV/c, whereas the second was performed at beam momenta of 1.657, 1.662 and 1.774 GeV/c. We will use $`1.455\pm 0.020`$ GeV/c and $`1.715\pm 0.060`$ GeV/c as momenta for the data sets.

The PS185 experimental set-up is shown in Fig. 2. It employs a segmented target consisting of five cylindrical cells each with a length and diameter of 2.5 mm. The first cell is made from pure carbon whereas the others consist of polyethylene (CH2). The 5 cells are embedded and separated by scintillators which are used as vetoes to allow for triggering on events with neutral final particles. A two-layer scintillator hodoscope detects the charged decay products. The sequence of charged $``$ neutral $``$ charged states forms the signature of the desired events and is the basis for the triggering.

A chamber stack, consisting of multiwire proportional chambers (MWPC) and drift chambers (DC), is located between the target and the hodoscope. A set of three drift chamber planes, contained inside a magnetic field and placed downstream of the hodoscope, allows a charge determination. The backward region is covered by an arrangement of limited streamer tubes.

Charged particle tracks are reconstructed in the chamber stack and combined pairwise to form neutral decay $`V^0`$’s. A kinematical fit is then performed on the $`V^0`$’s assuming that they correspond to $`K_S\pi ^+\pi ^{}`$ decays. It is based on a procedure described in Ref. , which has been extended to handle the case where only angles and no momenta are measured (cf. Ref. ). Monte Carlo simulations show that the resulting uncertainty in the reconstructed momenta is about 10%.

Both the $`K_SK_SX`$ and the $`\mathrm{\Lambda }\overline{\mathrm{\Lambda }}X`$ reaction channels lead to events with two $`V^0`$’s in the final state and it is not possible to separate these channels using the kinematical fit alone since no momenta are measured. However, the information from the energy loss in the hodoscope can be used to differentiate between them. Fig. 3a shows a scatter plot of the energy losses of particles in the hodoscope versus their momenta obtained from the kinematical fit assuming a $`K_S`$ decay. Bands belonging to pions and (anti-)protons are clearly visible. The protons originate from $`\mathrm{\Lambda }`$($`\overline{\mathrm{\Lambda }}`$) decays and constitute background events. As a consistency check, these background events are refitted assuming $`\mathrm{\Lambda }`$($`\overline{\mathrm{\Lambda }}`$) decays. The results given in Fig. 3b confirm the (anti-)proton assignment for these events.

We obtained a total sample of 31,732 $`K_SK_SX`$ events originating from either $`\overline{p}p`$ or $`\overline{p}`$C interactions from the two data-sets. Due to the relatively sharp lower edge of the Landau distribution, the contamination in the event sample arising from the hyperon background is estimated to be less than 2%. As a check of the final event sample, we determined the lifetime for the assigned kaons. Based on the measured decay length and the fitted momentum, we derived a value of $`(93\pm 4)`$ ps, which is in good agreement with the PDG value of 89.35 ps .

Fig. 4 illustrates the missing mass distributions for the extracted $`\overline{p}AK_SK_SX`$ events. They show a broad distribution stretching from 0 to 1.2 and 1.4 GeV/c2, respectively, which is due to $`K_SK_SX`$ production on 12C, involving the production of $`K^{}`$’s and $`\pi ^0`$’s, as well as combinatorial background. In addition, three peaks are visible in the spectra from the CH2 cells at the masses of the $`\pi ^0`$, $`\eta `$ and $`\omega `$ mesons (lower row) which are not found in the spectra from the carbon cell (upper row). These peaks are attributed to the $`\overline{p}pK_SK_S\pi ^0`$, $`\overline{p}pK_SK_S\eta `$ and $`\overline{p}pK_SK_S\omega `$ reactions.

## 3 RESULTS AND CONCLUSIONS

We have made fits to the missing mass distributions from the CH2 cells to obtain the number of events with $`\pi ^0`$’s, $`\eta `$’s and $`\omega `$’s in the final states arising from $`\overline{p}p`$ interactions. These involve a 5th order polynomial for the background and three gaussians; the results for the two data-sets are shown in Fig. 4.

The central values of the $`\eta `$ and $`\omega `$ peak are found to be $`(560\pm 14)`$ MeV/c2 and $`(786\pm 27)`$ MeV/c2 for the 1.455 GeV/c data and $`(553\pm 17)`$ MeV/c2 and $`(786\pm 19)`$ MeV/c2 for the 1.715 GeV/c data. The ratios of the number of $`\eta `$ to $`\omega `$ events are $`0.691\pm 0.139\pm 0.069`$ (1.455 GeV/c) and $`0.560\pm 0.160\pm 0.056`$ (1.715 GeV/c), where the first error is statistical and the second systematical. The latter, which is due to the uncertainty in the background subtraction, is estimated to be about 10%. The $`\pi ^0`$ case requires a slightly different treatment which will be described later.

The detector has the highest efficiency for forward-going (charged) particles. Due to the influence of the $`\eta `$/$`\omega `$ mass difference on the kinematics, we expect a larger acceptance for $`\overline{p}pK_SK_S\eta `$ than for $`\overline{p}pK_SK_S\omega `$ events. This difference is less important at higher energies.

The acceptances have been estimated by generating events with phase space distributions and applying cuts representing the geometrical acceptance and reconstruction efficiency. From these considerations we obtain relative efficiencies $`ϵ(\omega )/ϵ(\eta )=0.791\pm 0.040`$ at 1.455 GeV/c and $`0.882\pm 0.044`$ at 1.715 GeV/c. These results are consistent with the expectations mentioned above.

Different sets of cuts were tried and it was found that the relative efficiencies are not very sensitive to the specific set. This is not unexpected since the topologies for the two types of events are very similar provided that the differential distributions of the two reactions are rather alike. We investigated this by looking at the distribution of the c.m. scattering angles. We defined six regions in the missing mass spectrum: one covering each of the $`\eta `$ and $`\omega `$ peaks and one on either side of them. The c.m. scattering angle distributions for these regions are not significantly different from each other. Furthermore, the experimental momentum and scattering angle distributions and those obtained from Monte Carlo simulations based on phase space distributions are quite similar.

We rely on the known $`\overline{p}pK_SK_S\omega `$ cross section to obtain the absolute normalisation of the $`\overline{p}pK_SK_S\eta `$ cross section. Data on the $`\overline{p}pK_SK_S\omega `$ cross section are found in Refs.  and summarised in Ref. . It is worth noting that all the data, except those of Ganguli et al. , have to be corrected for unseen decay modes before the total cross section can be obtained. For this analysis, we confined ourselves to the data in the $`\overline{p}`$ momentum range from 1.0 to 2.2 GeV/c. The cross sections given in Ref.  have not been included since they are averages of values already given in Ref. . The remaining 18 data points shown in Fig. 5 have been corrected for the unseen $`\omega `$ decay modes using the branching ratio BR$`(\omega \pi ^+\pi ^{}\pi ^0)=0.89`$ . A linear fit is sufficient to describe the data in this region ($`\chi _{red}^2=0.65`$). From the fit we obtain $`\sigma \left(\overline{p}pK_SK_S\omega \right)=(50.0\pm 2.5)\mu `$b at $`p_{\overline{p}}=1.455`$ GeV/c and $`(39.7\pm 2.0)\mu `$b at $`p_{\overline{p}}=1.715`$ GeV/c.

The cross section for the $`\overline{p}pK_SK_S\eta `$ reaction is finally deduced using the relation

$$\sigma \left(\overline{p}pK_SK_S\eta \right)=\frac{\text{number of }\eta \text{events}}{\text{number of }\omega \text{events}}\times \frac{ϵ(\omega )}{ϵ(\eta )}\times f_{neut}\times \sigma \left(\overline{p}pK_SK_S\omega \right),$$
(1)
where $`ϵ(\omega )`$ and $`ϵ(\eta )`$ denote the acceptances of the corresponding reactions and $`f_{neut}`$ is the relative suppression of $`\omega `$ events compared to $`\eta `$ events due to the different neutral branching ratios. The branching ratio for neutral decays is $`(9.0\pm 0.7)`$% for the $`\omega `$ and $`(71.6\pm 0.4)`$% for the $`\eta `$  leading to $`f_{neut}=0.126\pm 0.010`$.

Inserting the numbers quoted above into Eq. 1 leads to cross sections for $`\overline{p}pK_SK_S\eta `$ of $`(3.44\pm 0.69\pm 0.50)\mu `$b at $`1.455\pm 0.020`$ GeV/c and $`(2.47\pm 0.71\pm 0.36)\mu `$b at $`1.715\pm 0.060`$ GeV/c, where the first uncertainty is statistical and the second systematical.

Many systematic uncertainties cancel when forming ratios as in Eq. 1. As a cross check of the method, we derived also the $`\overline{p}pK_SK_S\pi ^0`$ cross section in a similar way. We obtained the number of $`K_SK_S\pi ^0`$ events from the missing mass plot. The sample also contains some $`\overline{p}pK_SK_S`$ events lying under the $`\pi ^0`$ peak. The cross section for this reaction is small, about 2 $`\mu `$, and this, together with the lower acceptance, allows us to neglect its contribution within the precision of the present experiment.

The ratios of the number of $`\pi ^0`$ to $`\omega `$ events are $`6.23\pm 0.86\pm 0.62`$ (1.455 GeV/c) and $`3.85\pm 0.61\pm 0.38`$ (1.715 GeV/c). The relative efficiency for measuring the $`\overline{p}pK_SK_S\omega `$ and $`\overline{p}pK_SK_S\pi ^0`$ reactions becomes $`0.80\pm 0.04`$ and $`0.93\pm 0.05`$ at the two momenta and $`f_{neut}=0.091\pm 0.007`$. With these values, we find $`(22.7\pm 3.1\pm 3.3)\mu `$b at $`1.455\pm 0.020`$ GeV/c and $`(12.9\pm 2.0\pm 1.9)\mu `$b at $`1.715\pm 0.060`$ GeV/c for the $`\overline{p}pK_SK_S\pi ^0`$ cross section. The good agreement between these results and the published data  shown in Fig. 6 gives us further confidence in the analysis.

The only published data on $`\overline{p}pK_SK_S\eta `$ is the branching ratio for annihilation at rest, BR$`(\overline{p}pK_SK_S\eta )=(0.25\pm 0.04)\times 10^3`$ . The corresponding value for the $`K_SK_S\omega `$ annihilation channel is BR$`(\overline{p}pK_SK_S\omega )=(1.16\pm 0.06)\times 10^3`$ . In order to compare these data with the results of our measurement, we extrapolated the ratio of the $`K_SK_S\eta `$ and $`K_SK_S\omega `$ cross sections from the threshold value to the momentum range between 1.2 and 2.0 GeV/c. To do this, we retained only phase space factors, assuming that the ratio of the matrix elements for both transitions is constant. Such considerations predict cross section ratios of $`0.043\pm 0.007`$ at 1.455 GeV/c and $`0.039\pm 0.006`$ at 1.715 GeV/c. Our data points lie within one standard deviation above this estimate.

In conclusion, we have made the first measurement of the cross section for the $`\overline{p}pK_SK_S\eta `$ reaction in flight in two regions of $`\overline{p}`$ momenta, around 1.455 and 1.715 GeV/c. The magnitude of the cross section is in reasonable agreement with an extrapolation from the known $`\overline{p}pK_SK_S\eta `$ and $`\overline{p}pK_SK_S\omega `$ branching ratios at rest using phase space arguments.

We hope that these results will stimulate further theoretical efforts to shed light on the process of strangeness production.

The PS185 team wants to thank the LEAR accelerator team for providing excellent experimental conditions. We would also like to thank Professor Colin Wilkin for a very helpful critique of our manuscript. We also gratefully acknowledge the financial and material support from the German Bundesministerium für Bildung und Forschung, the Swedish Natural Science Research Council, the United States Department of Energy and the United States Science Foundation. This work is based in part on the dissertation of S. Pomp submitted to the Uppsala University in partial fulfillment of the requirements for the PhD degree.


format: text

text:
# Determining the density of states for classical statistical models: A random walk algorithm to produce a flat histogram

## I Introduction

Computer simulation now plays a major role in statistical physics , particularly for the study of phase transitions and critical phenomena. The standard Monte Carlo (MC) method has been the Metropolis importance sampling algorithm , but more recently new, efficient algorithms have begun to play a role in allowing simulation to achieve the resolution which is needed to accurately locate and characterize phase transitions . For example, cluster flip algorithms, beginning with the seminal work of Swendsen and Wang  and extended by Wolff , have been used to reduce critical slowing down near 2nd order transitions. Similarly, the multicanonical ensemble method  was introduced to overcome the tunneling barrier between coexisting phases at 1st order transitions and has general utility for systems with a rough energy landscape . In both situations, histogram re-weighting techniques  can be applied in the analysis to increase the amount of information that can be gleaned from simulational data, but the applicability of re-weighting is severely limited in large systems by the statistical quality of the “wings” of the histogram. This latter effect is quite important in systems with competing interactions for which short range order effects might occur over very broad temperature ranges or even give rise to frustration that produces a very complicated energy landscape and limits the efficiency of “standard” methods.

One of the most important quantities in statistical physics is the density of states $`g(E)`$, i.e. the number of all possible states (or configurations) for an energy level $`E`$ of the system, but direct estimation of this quantity has not been the goal of simulations. Instead, most conventional Monte Carlo algorithms  such as Metropolis importance sampling, Swendsen-Wang cluster flipping, etc. generate a canonical distribution $`g(E)e^{E/k_BT}`$ at a given temperature. Such distributions are so narrow that, with conventional Monte Carlo simulations, multiple runs are required if we want to know thermodynamic quantities over a significant range of temperatures. In the canonical distribution, the density of states does not depend on the temperature at all. If we can estimate the density of states $`g(E)`$ with high accuracy for all energies, we can then construct canonical distributions at essentially any temperature. For a given model in statistical physics, once the density of states is known we can calculate the partition function as $`Z=\underset{E}{}g(E)e^{\beta E}`$, and the model is essentially “solved” since most thermodynamic quantities can be calculated from it. Though computer simulation is already a very powerful method in statistical physics , it seems that there is no efficient algorithm to calculate the density of states very accurately for large systems. Even for exactly solvable models such as the 2-dim Ising model, $`g(E)`$ is impossible to calculate exactly for a large system .

The multicanonical ensemble method  proposed by Berg et al. has been proved to be very efficient in studying first-order phase transitions where simple canonical simulations have difficulty overcoming the tunneling barrier between coexisting phases at the transition temperature . The method also has been successfully applied to some complex systems, such as spin glass models  and protein folding problems , for which the energy landscape is very rough and the conventional canonical Monte Carlo simulation gets easily trapped in local minima at low temperature. In the multicanonical method, we have to estimate the density of states $`g(E)`$ first, then perform a random walk with a probability of $`\frac{1}{g(E)}`$ to make the histogram flat in the desired region in the phase space, such as between two peaks of the canonical distribution at the first-order transition temperature. In a multicanonical simulation, the density of states need not necessarily be very accurate, as long as the simulation generates a relatively flat histogram and overcomes the tunneling barrier in energy space. This is because the subsequent re-weighting , does not depend on the accuracy of the density of the states as long as the histogram can cover all important energy levels with sufficient statistics. (If the density of states could be calculated very accurately, then the problem would have been solved in the first place and we needn’t perform any further simulation such as with the multicanonical simulational method.) Berg et al. proposed a recursive method to calculate the density of states by accumulating histogram entries and estimating the density of states iteratively with the histogram data . Their method works well for small systems since the number of all possible states is small; however, for large systems the number of possible states increases exponentially with the size of a system. A simple iteration method to construct a histogram gets trapped in a narrow energy range, and it needs an extremely long time and large number of iterations to get an accurate estimate for the density of states for the entire energy space, even for a small system. It is thus not practical to calculate the density of states for large systems with this approach. For a simple model such as the Potts model, Berg and Neuhaus used finite-size scaling theory to “guess” the density of states up to $`100\times 100`$ from the simulational results of the small systems . Such results are not as reliable as direct simulations, moreover for complex systems such as spin glass models, we can not simply apply such finite-size scaling at all. Berg and Celik applied their multicanonical method to a 2D spin-glass model, but they could test their method only up to a $`48\times 48`$ lattice . Recently Berg and Janke proposed a similar method (multioverlap simulational) for the 3D Ising spin-glass model, they obtained some reliable results for systems as large as $`12\times 12\times 12`$ .

Lee  independently proposed the entropic sampling method, which is basically equivalent to multicanonical ensemble sampling. He used an iteration process to calculate the microcanonical entropy at $`E`$ which is defined by $`S(E)=\mathrm{ln}[g(E)]`$ where $`g(E)`$ is the density of states. He also applied his method to the 2D ten-state ($`Q=10`$) Potts model and the 3D Ising model; however, just as for other simple iteration methods, it works well only for small systems. He obtained a good result with his method for the $`24\times 24`$ 2-dim $`Q=10`$ Potts model and the $`4\times 4\times 4`$ 3-dim Ising model.

Oliveira et al.  proposed the broad histogram method with which they calculated the density of states by estimating the probabilities of possible transitions between all possible states of a random walk in energy space. Using simple canonical average formulae in statistical physics, they then calculated thermodynamic quantities for any temperature. Though the authors believed that the broad histogram relation is exact, their simulational results have systematic errors even for the Ising model on a lattice as small as $`32\times 32`$ in references . They believed that the error was due to the particular dynamics adopted within the broad histogram method , but other work  argued that it violates the detailed balance condition. The algorithm was corrected in reference  and an approach to the broad-histogram method was proposed to calculate the density of states based on the number of potential moves during the random walk in energy space . They obtained more accurate results for a $`32\times 32`$ Ising model than the broad histogram did; but this method also suffers from the systematic errors and substantial deviations when system becomes larger than $`32\times 32`$ .

It is thus an extremely difficult task to calculate density of states directly with high accuracy for large systems. All methods based on accumulation of histogram entries, such as the histogram method of Ferrenberg and Swendsen , Lee’s version of multicanonical method (entropic sampling) , broad histogram method  and flat histogram method  have the problem of scalability for large systems. These methods suffer from systematic errors when systems are large, so we still need a superior algorithm to calculate the density of states for large systems.

Very recently, we introduced a new, general, efficient Monte Carlo algorithm that offers substantial advantages over existing approaches . In this paper, we will explain the algorithm in detail, including our implementation, and describe its application not only to 1st and 2nd order phase transitions, but also to a 3D spin glass model that has a rough energy landscape.

Unlike conventional Monte Carlo methods that directly generate a canonical distribution at a given temperature $`g(E)e^{E/K_\text{B}T}`$, our approach is to directly estimate the density of states $`g(E)`$ accurately via a random walk that produces a flat histogram in energy space. We modify density of states at each step of the random walk, and by carefully controlling the modification factor we can obtain a density of states that converges to the real value very quickly even for large systems. The resultant density of states is accurate enough to calculate thermodynamic quantities by applying canonical average formulas in statistical physics.

The remainder of this paper is arranged as follows. In section II, we present our general algorithm in detail. In section III, we apply our method to the 2D $`Q=10`$ Potts model which has a first-order phase transition. In section IV, we apply our method to a model with a second-order phase transition to test the accuracy of the algorithm. In section V, we consider the 3D $`\pm J`$ spin glass model, a system with rough landscapes. Discussion and the conclusion are presented in section VI.

## II A general and efficient algorithm to estimate the density of states with a flat histogram

Our algorithm is based on the observation that if we perform a random walk in energy space by flipping spins randomly for a spin system, and the probability to visit a given energy level $`E`$ is proportional to the reciprocal of the density of states $`\frac{1}{g(E)}`$, then a flat histogram is generated for the energy distribution. This is accomplished by modifying the estimated density of states in a systematic way to produce a “flat” histogram over the allowed range of energy and simultaneously making the density of states converge to the true value. We modify the density of states constantly during each step of the random walk and use the updated density of states to perform a further random walk in energy space. The modification factor of the density of states is controlled carefully, and at the end of simulation the modification factor should be very close to $`1`$ which is the ideal case of the random walk with the true density of states.

At the very beginning of our simulation, the density of states is a priori unknown, so we simply set all entries to $`g(E)=1`$ for all possible energies $`E`$. Then we begin our random walk in energy space by flipping spins randomly and the probability at a given energy level is proportional to $`\frac{1}{g(E)}`$. In general, if $`E_1`$ and $`E_2`$ are energies before and after a spin is flipped, the transition probability from energy level $`E_1`$ to $`E_2`$ is:

$$p(E_1E_2)=\mathrm{min}[\frac{g(E_1)}{g(E_2)},1].$$
(1)
Each time an energy level $`E`$ is visited, we modify the existing density of states by a modification factor $`f>1`$, i.e. $`g(E)g(E)f`$. (In practice, we use the formula $`\mathrm{ln}[g(E)]\mathrm{ln}[g(E)]+\mathrm{ln}(f)`$ in order to fit all possible $`g(E)`$ into double precision numbers for the systems we will discuss in this paper.) If the random walk rejects a possible move and stays at the same energy level, we also modify the existing density of states with the same modification factor. Throughout this study we have used an initial modification factor of $`f=f_0=e^12.71828\mathrm{}`$ which allows us to reach all possible energy levels very quickly even for a very large system. If $`f_0`$ is too small, the random walk will spend an extremely long time to reach all possible energies. However, too large a choice of $`f_0`$ will lead to large statistical errors. In our simulations, the histograms are generally checked about each 10000 MC sweeps. A reasonable choice is to make $`f_0^{10000}`$ have the same order of magnitude as the total number of states ($`Q^N`$ for a Potts model). During the random walk, we also accumulate the histogram $`H(E)`$ (the number of visits at each energy level $`E`$) in the energy space. When the histogram is “flat” in the energy range of the random walk, we know that the density of states converges to the true value with an accuracy proportional to that modification factor $`\mathrm{ln}(f)`$. Then we reduce the modification factor to a finer one using a function like $`f_1=\sqrt{f_0}`$, reset the histogram, and begin the next level random walk during which we modify the density of states with a finer modification factor $`f_1`$ during each step. We continue doing so until the histogram is “flat” again and then reduce the modification factor $`f_{i+1}=\sqrt{f_i}`$ and restart. We stop the random walk when the modification factor is smaller than a predefined value (such as $`f_{\text{final}}=\mathrm{exp}(10^8)1.00000001`$). It is very clear that the modification factor acts as a most important control parameter for the accuracy of the density of states during the simulation and also determines how many MC sweeps are necessary for the whole simulation. The accuracy of the density of states depends on not only $`f_{\text{final}}`$ but also many other factors, such as the complexity and size of the system, criterion of the flat histogram and other details of the implementation of the algorithm.

It is impossible to obtain a perfectly flat histogram and the phrase “flat histogram” in this paper means that histogram $`H(E)`$ for all possible $`E`$ is not less than $`x\%`$ of the average histogram $`H(E)`$, where $`x\%`$ is chosen according to the size and complexity of system and the desired accuracy of the density of states. For the $`L=32`$, 2D Ising model with only nearest-neighbor couplings, this percentage can be chosen as high as $`95\%`$, but for large systems the criterion for “flatness” may never be satisfied if we choose too high a percentage and the program may run forever.

One essential constraint on the implementation of the algorithm is that the density of states during the random walk converges to the true value. The algorithm proposed in this paper has this property. The accuracy of the density of states is proportional to $`\mathrm{ln}(f)`$ at that iteration; however, $`\mathrm{ln}(f_{\text{final}})`$ can not be chosen arbitrary small or the modified $`\mathrm{ln}[g(E)]`$ will not differ from the unmodified one to within the number of digits in the double precision numbers used in the calculation. If this happens, the algorithm no longer converges to the true value, and the program may run forever. Even if $`f_{\text{final}}`$ is within range but too small, the calculation might take excessively long to finish.

We have chosen to reduce the modification factor by a square root function, and $`f`$ approaches $`1`$ as the number of iterations approaches infinity. In fact, any function may be used as long as it decreases $`f`$ monotonically to 1. A simple and efficient formula is $`f_{i+1}=f_i^{1/n}`$, where $`n>1`$. The value of $`n`$ can be chosen according to the available CPU time and expected accuracy of the simulation. For the systems that we have studied the choice of $`n=2`$ yielded good accuracy in a relatively short time, even for large systems. When the modification factor is almost $`1`$ and the random walk generates a uniform distribution in energy space, the density of states should converge to the true value for the system.

The method can be further enhanced by performing multiple random walks, each for a different range of energy, either serially or in parallel fashion. We restrict the random walk to remain in the range by rejecting any move out of that range. The resultant pieces of the density of states can then be joined together and used to produce canonical averages for the calculation of thermodynamic quantities at essentially any temperature.

We should point out here that during the random walk (especially at the early stage of iteration) in the energy space, the algorithm does not satisfy the detailed balance condition exactly, since the density of states is modified constantly during the random walk. However, after many iterations, the density of states converges to the true value very quickly as the modification factor approaches $`1`$. If $`p(E_1E_2)`$ is the transition probability from the energy level $`E_1`$ to level $`E_2`$, from equation (1), the ratio of the transition probabilities from $`E_1`$ to $`E_2`$ and from $`E_2`$ to $`E_1`$ can be calculated very easily as:

$$\frac{p(E_1E_2)}{p(E_2E_1)}=\frac{g(E_1)}{g(E_2)}$$
(2)
where $`g(E)`$ is the density of states. In another words, our random walk algorithm satisfies the detailed balance condition:

$$\frac{1}{g(E_1)}p(E_1E_2)=\frac{1}{g(E_2)}p(E_2E_1)$$
(3)
where $`\frac{1}{g(E_1)}`$ is the probability at the energy level $`E_1`$ and $`p(E_1E_2)`$ is the transition probability from $`E_1`$ to $`E_2`$ for the random walk. We conclude that the detailed balance condition is satisfied with accuracy proportional to the modification factor $`\mathrm{ln}(f)`$.

Almost all recursive methods update the density of states by using the histogram data directly only after enough histogram entries are accumulated . Because of the exponential growth of the density of states in energy space, this process is not efficient because the histogram is accumulated linearly. In our algorithm, we modify the density of states at each step of the random walk, and this allows us to approach the true density of states much faster than conventional methods especially for large systems. (We also accumulate histogram entries during the random walk, but we only use it to check whether the histogram is flat enough to go to the next level random walk with a finer modification factor.)

We should point out here that the total number of configurations increases exponentially with the size of the system; however, the total number of possible energy levels increases linearly with the size of system. It is thus easy to calculate the density of states with a random walk in energy space for a large system. In this paper for an example, we consider the Potts model on a $`L\times L`$ lattice with nearest-neighbor interactions . For $`Q3`$, the number of possible energy levels is about $`2N`$, where $`N=L^2`$ is the total number of the lattice site. However, the average number of possible states (or configurations) on each energy level is as large as $`\frac{Q^N}{2N}`$, where $`Q`$ is the number of possible states of a Potts spin and $`Q^N`$ is the total number of possible configurations of the system. This is the reason why most models in statistical physics are well defined, but we can not simply use our computers to realize all possible states to calculate any thermodynamic quantities, this is also the reason why efficient and fast simulational algorithms are required in the numerical investigations.

By the end of simulation, we only obtain relative density, since the density of states can be modified at each time it is visited. We can apply the condition that the total number of possible states for the $`Q`$ state Potts model is $`\underset{E}{}g(E)=Q^N`$ or the number of ground state is $`Q`$ to get the absolute density of states.

## III Application to a first order phase transition

### A Potts model and its canonical distribution

In this section, we apply our algorithm to a model with a first-order phase transition . We choose the 2D ten state Potts model since it serves as an ideal laboratory for temperature-driven first-order phase transitions. Since some exact solutions and extensive simulational data are available, we have ample opportunity to compare our results with other values.

We consider the 2-dimensional $`Q=10`$ Potts model on $`L\times L`$ square lattice with nearest-neighbor interactions and periodic boundary conditions. The Hamiltonian for this model can be written as:

$$=\underset{<ij>}{}\delta (q_i,q_j)$$
(4)
and $`q=1,2,\mathrm{}Q`$. During the simulation, we select lattice sites randomly and choose integers between $`[1:Q]`$ randomly for new Potts spin values. The modification factor $`f_i`$ changes from $`f_0=e^1=2.71828`$ at the very beginning to $`f_{\text{final}}=\mathrm{exp}(10^8)1.00000001`$ by the end of the random walk. At the end of the simulations, our algorithm only provides a relative density of states for different energies, so to extract the correct density of states, we can either use the fact that the total number of possible states is $`Q^N`$ or that the number of ground states is $`Q`$, where $`N=L^2`$ is the total number of lattice sites. (Actually we can use one of these two conditions to get the absolute density of states, and use the other condition to check the accuracy of our result.) To guarantee the accuracy of thermodynamic quantities at low temperatures in further calculations, in this paper we use the condition that the number of the ground states is $`Q`$ to normalize the density of states. The densities of states for $`100\times 100`$, $`150\times 150`$ and $`200\times 200`$ are shown in Fig. 1. It is very clear from the figure that the maximum density of states for $`L=200`$ is very close to $`10^{40000}`$ which is actually about $`5.75\times 10^{39997}`$ from our simulational data.

Conventional Monte Carlo simulation (such as Metropolis sampling ) realizes a canonical distribution $`P(E,T)`$ by generating a random walk Markov chain at a given temperature:

$$P(E,T)=g(E)e^{E/k_BT}$$
(5)
From the simulational result for the density of states $`g(E)`$, we can calculate the canonical distribution by the above formula at essentially any temperature without performing multiple simulations. In Fig. 2(a), we show the resultant double peaked canonical distribution , at the transition temperature $`T_\text{c}`$ for the first-order transition of the $`Q=10`$ Potts model. The “transition temperatures” $`T_\text{c}(L)`$ are 0.70171 for $`L=60`$, 0.70143 for $`L=80`$ and 0.70135 for $`L=100`$ which are determined by the temperatures where the double peaks are of the same height. Note that the peaks of the distributions are normalized to 1 in this figure. The valley between two peaks is quite deep. e.g. is $`7\times 10^5`$ for $`L=100`$. The latent heat for this temperature driven first-order phase transition can be estimated from the energy difference between the double peaks. Our results for the locations of the peaks are listed in the Table 1. They are consistent with the results obtained by multicanonical method  and multibondic cluster algorithm  for those lattice sizes for which these other methods are able to generate estimates. As the table shows, our method produces results for substantially larger systems than have been studied by these other approaches.

Because of the double peak structure at a first-order phase transition, conventional Monte Carlo simulations are not efficient since an extremely long time is required for the system to travel from one peak to the other in energy space. With the algorithm proposed in this paper, all possible energy levels are visited with equal probability, so it overcomes the tunneling barrier between the coexisting phases in the conventional Monte Carlo simulations. The histograms for $`L`$ = 60, 80 and 100 are shown in the inset of the Fig. 2(a) and are very flat. The histogram in the figure is the overall histogram defined by the total number of visits to each energy level for the random walk. Here, too, we choose the initial modification factor $`f_0=e^1`$, and the final one as $`\mathrm{exp}(10^8)1.00000001`$; and the total number of iterations is 27. In our simulation, we do not set a predetermined number of MC sweeps for each iteration but rather give the criterion which the program checks periodically. Generally, the number of MC sweeps needed to satisfy the criterion increases as we reduce the modification factor to a finer one, but we cannot predict the exact number MC sweeps needed for each iteration before the simulation. We believe that it is preferable to allow the program to decide how great a simulational effort is needed for a given modification factor $`f_i`$. This also guarantees a sufficiently flat histogram resulting from a random walk which in turn determines the accuracy of the density of states at the end of the simulation. We nonetheless need to perform some test runs to make sure that the program will finish within a given time. The entire simulational effort used was about $`1.6\times 10^7`$ visits per energy level (or $`3.2\times 10^7`$ MC sweeps) for $`L=60`$, $`2.2\times 10^7`$ visits for $`L=80`$ and $`3.3\times 10^7`$ visits for $`L=100`$. With the program we implemented, the simulation for $`L=100`$ can be completed within two weeks in a single 600 MHz Pentium III processor.

To speed up the simulation, we needn’t constrain ourselves to performing a single random walk over the entire energy range with high accuracy. If we are only interested in a specific temperature range, such as near $`T_{\text{}}`$, we could first perform a low precision unrestricted random walk, i.e. over all energies, to estimate the required energy range, and then carry out a very accurate random walk for the corresponding energy region. The inset of Fig. 2(a) only shows the histograms for the extensive random walks in the energy range between $`E/N=1.90`$ and $`0.6`$. If we need to know the density of states more accurately for some energies, we also can perform separate simulations, one for low energy levels, one for high energy levels, the other for middle energy which includes double peaks of the canonical distribution at $`T_\text{c}`$. This scheme not only speeds up the simulation, but also increases the probability of accessing the energy levels for which both maximum and minimum values of the distributions occur by performing the random walk in a relatively small energy range. If we perform single random walk over all possible energies, it will take a long time to generate rare spin configurations. Such rare energy levels include the ground energy level or low energy levels with only few spins with different values and high energy levels where all, or most, adjacent Potts spins have different values.

With the algorithm in this paper, if the system is not larger than $`100\times 100`$, the random walk on important energy regions (such as that which includes the two peaks of the canonical distribution at $`T_\text{c}`$) can be carried out with a single processor and will give an accurate density of states within about $`10^7`$ visits per energy level. The results shown in Fig. 2(a) were obtained using a single processor. However, for a larger system, we can use a parallelized algorithm by performing random walks in different energy regions, each using a different processor. We have implemented this approach using PVM with a simple master-slave model and can then obtain an accurate estimate for the density of states with relatively short runs on each processor. The densities of states for $`150\times 150`$ and $`200\times 200`$, shown in Fig. 1, were obtained by joining together the estimates obtained from $`21`$ independent random walks, each constrained within a different regions of energy. The histograms from the individual random walks are shown in the inset of Fig. 2(b) both for $`150\times 150`$ and $`200\times 200`$ lattices. In this case, we only require that the histogram of the random walk in the corresponding energy segment is sufficiently flat without regard to the relative flatness over the entire energy range. In Fig. 2(b), the results for large lattices show clear double peaks for the canonical distributions at temperatures $`T_\text{c}(L)=0.70127`$ for $`L=150`$ and $`T_\text{c}(L)=0.701243`$ for $`L=200`$. The exact result is $`T_\text{c}=0.701232\mathrm{}.`$ for the infinite system. Considering the valley which we find for $`L=200`$ is as deep as $`9\times 10^{10}`$, we can understand why it is impossible for conventional Monte Carlo algorithms to overcome the tunneling barrier with available computational resources.

If we compare the histogram for $`L=60`$ in Fig. 2(a) with that for $`L=200`$ in Fig. 2(b), we see very clearly that the simulation effort for $`L=200`$ ($`9.8\times 10^6`$ visits per energy level) is even less then the effort for $`L=60`$ ($`1.6\times 10^7`$ visits per energy level.). It is more efficient to perform random walks in relatively small energy segments than a single random walk over all energies. The reason is very simple, the random walk is a local walk, which means for a given $`E_1`$, the energy level for the next step only can be one of 9 levels in the energy range $`[E_14,E_1+4]`$. (for the Potts model discussed in this section). The algorithm itself only requires that the histogram on such local transitions is flat. (A single random walk, subject to the requirement of a flat histogram for all energy levels, will take quite long.) For random walks in small energy segments, we should be very careful to make sure that all spin configurations with energies in the desired range can be equally accessed so we restart the random walk periodically from independent spin configurations.

An important question that must be addressed is the ultimate accuracy of the algorithm. One simple check is to estimate the transition temperature of the 2D $`Q=10`$ Potts model for $`L=\mathrm{}`$ since the exact solution is known. According to finite-size scaling theory, the “effective” transition temperature for finite systems behaves as:

$$T_\text{c}(L)=T_\text{c}(\mathrm{})+\frac{c}{L^d}$$
(6)
where $`T_\text{c}(L)`$ and $`T_\text{c}(\mathrm{})`$ are the transition temperatures for finite- and infinite-size systems, respectively, $`L`$ is the linear size of the system and $`d`$ is dimension of the lattice.

In Fig. 3, the transition temperature is plotted as a function of $`L^d`$. The data in the main portion of the figure are obtained from small systems($`L=1030`$), and the error bars are estimated by results from multiple independent runs. Clearly the transition temperature extrapolated from our simulational data is $`T_\text{c}(\mathrm{})=0.7014\pm 0.0004`$ which is consistent with the exact solution ($`T_\text{c}=0.701232\mathrm{}.`$) for the infinite system. To get an even more accurate estimate, and also test the accuracy of the density of states from single runs for large systems, we performed single, long random walk on large lattices($`L=50200`$). The results, plotted as a function of lattice size in the inset of the figure, show that the transition temperature extrapolated from the finite systems is $`T_\text{c}(\mathrm{})=0.701236\pm 0.000025`$ which is still consistent with the exact solution.

We also compare our simulational result for the $`Q=10`$ Potts model with the existing numerical data such as estimates of transition temperatures and double peak locations obtained with the multicanonical simulational method by Berg and Neuhaus  and the Multibondic cluster algorithm by Janke and Kappler . All results are shown in Table 1. With our random walk simulational algorithm, we can calculate the density of states up to $`200\times 200`$ within $`10^7`$ visits per energy level to obtain a good estimate of the transition temperature and locations of the double peaks. Using the multicanonical method and a finite scaling guess for the density of states, Berg et al. only obtained results for lattices as large as $`100\times 100`$ , and multibondic cluster algorithm data  were not given for systems larger than $`50\times 50`$.

In section IV, the accuracy of our algorithm will be further tested by comparing thermodynamic quantities obtained for 2-dim Ising model with exact solutions.

### B Thermodynamic properties of the $`Q=10`$ Potts model

One of advantages of our method is that the density of states does not depend on temperature; indeed with the density of states, we can calculate thermodynamic quantities at any temperature. For example, the internal energy can be calculated by:

$`U(T)={\displaystyle \frac{\underset{E}{}Eg(E)e^{\beta E}}{\underset{E}{}g(E)e^{\beta E}}}E_T`$ (7)
To study the behavior of the internal energy near $`T_{extc}`$ more carefully, we calculate the internal energy for $`L=60`$, 100 and 200 near $`T_\text{c}`$ as presented in Fig. 4(a). A very sharp “jump” in the internal energy at transition temperature $`T_\text{c}`$ is visible, and the magnitude of this jump is equal to the latent heat for the (first-order) phase transition. Such behavior is related to the double peak distribution of the first-order phase transition. When $`T`$ is slightly away from $`T_\text{c}`$, one of the double peaks increases dramatically in magnitude and the other decreases.

Since we only perform simulations on finite lattices, and use a continuum function to calculate thermodynamic quantities, all our quantities for finite-size systems will appear to be continuous if we use a very small scale. In the inset of Fig. 4(a) we use the same density of states again to calculate the internal energy for temperatures very close to $`T_\text{c}`$. On this scale the “discontinuity” at the first-order phase transition disappears and a smooth curve can be seen instead of a sharp “jump” in the main portion of Fig. 4(a). The discontinuity in Fig. 4(a) is simply due to the coarse scale, but when the system size goes to infinity, the discontinuity will be real.

From the density of states we can also estimate the specific heat from the fluctuations in the internal energy:

$$C(T)=\frac{U(T)}{T}=\frac{E^2_TE_{T}^{}{}_{}{}^{2}}{T^2}$$
(8)
In Fig. 4(b), the specific heat so obtained is shown as a function of temperature. We calculate the specific heat in the vicinity of the transition temperature $`T_\text{c}`$. The finite-size dependence of the specific heat is clearly evident. We find that specific heat has a finite maximum value for a given lattice size $`L`$ that, according to finite-size scaling theory for first-order transitions should vary as:

$$c(L,T)L^df((TT_\text{c}(\mathrm{}))L^d)$$
(9)
where $`c(L,T)=C(L,T)/N`$ is the specific heat per lattice site, $`L`$ is the linear lattice size, $`d=2`$ is the dimension of the lattice. $`T(L=\mathrm{})=0.70123\mathrm{}.`$ is the exact solution for the $`Q=10`$ Potts model . In the inset of Fig. 4(b), our simulational data for systems with $`L=60`$, $`100`$ and $`200`$ can be well fitted by a single scaling function, moreover this function is completely consistent with the one obtained from lattice sizes from $`L=18`$ to $`L=50`$ by standard Monte Carlo .

With the density of states, we not only can calculate most thermodynamic quantities for all temperatures without multiple simulations but can also access some quantities, such as the Gibbs free energy and entropy, which are not directly available from conventional Monte Carlo simulations. The free energy is calculated using

$`Z`$ $`=`$ $`{\displaystyle \underset{\{\text{configurations}\}}{}}e^{\beta E}={\displaystyle \underset{E}{}}g(E)e^{\beta E}`$ (10)
$`F`$ $`=`$ $`kT\mathrm{log}(Z)`$ (11)
Our results for the Gibbs free energy per lattice site is shown in Fig. 4(c) as a function of temperature. Since the transition is first-order the free energy appears to have a “discontinuity” in the first derivative at $`T_\text{c}`$. This is typical behavior for a first-order phase transition, and even with the fine scale used in the inset of Fig. 4(c), this property is still apparent even though the system is finite. The transition temperature $`T_\text{c}`$ is determined by the point where the first derivative appears to be discontinuous. With a coarse temperature scale we can not distinguish the finite-size behavior of our model; however, we can see a very clear size dependence when we view the free energy on a very fine scale as in the inset of Fig. 4(c).

The entropy is another very important thermodynamic quantity that cannot be calculated directly in conventional Monte Carlo simulations. It can be estimated by integrating over other thermodynamic quantities, such as specific heat, but the result is not always reliable since the specific heat itself is not easy to determine accurately, particularly considering the “divergence” at the first-order transition. With an accurate density of states estimated by our method, we already know the Gibbs free energy and internal energy for the system, so the entropy can be calculated easily:

$$S(T)=\frac{U(T)F(T)}{T}$$
(12)
It is very clear that the entropy is very small at low temperature and at $`T=0`$ is given by the density of states for the ground state. We show the entropy as a function of temperature in a wide region in Fig. 4.(d)

The entropy has a very sharp “jump” at $`T_\text{c}`$, just as does the internal energy and such behavior can be seen very clearly in the inset of Fig. 4(d), when we re-calculate the entropy near $`T_{\text{}}`$. The change of the entropy at $`T_\text{c}`$ shown in the figure can be obtained by the latent heat divided by the transition temperature, and the latent heat can be obtained by the jump in internal energy at $`T_\text{c}`$ in Fig. 4(a).

With the histogram method proposed by Ferrenberg and Swendsen , it is possible to use simulational data at specific temperatures to obtain complete thermodynamic information near, or between, those temperatures. Unfortunately it is usually quite hard to get accurate information in the region far away from the simulated temperature due to difficulties in obtaining good statistics, especially for large systems where the canonical distributions are very narrow. With the algorithm proposed in this paper, the histogram is “flat” for the random walk and we always have essentially the same statistics for all energy levels. Since the output of our simulation is the density of states, which does not depend on the temperature at all, we can then calculate most thermodynamic quantities at any temperature without repeating the simulation. We also believe the algorithm is especially useful for obtaining thermodynamic information at low temperature or at the transition temperature for the systems where the conventional Monte Carlo algorithm is not so efficient.

### C The tunneling time for the $`Q=10`$ Potts model at $`T_\text{c}`$

To study the efficiency of our algorithm, we measure the tunneling time $`\tau `$, defined as the average number of sweeps needed to travel from one peak to the other and return to the starting peak in energy space. Since the histogram that our random walk produces is flat in energy space, we expect the tunneling time will be the same as for the ideal case of a simple random walk in real space, i.e. $`\tau (N_E)N_E`$, where $`N_E`$ is the total number of energy levels. To compare our simulational results to those for the ideal case, we also perform a random walk in real space. We always use a fixed $`g(x_i)=1`$ in one-dimensional real space, where $`x_i`$ is a discrete coordinate of position that can be chosen simply as 1, 2, 3 ,4 …$`N_E`$. The random walk is a local random walk with transition probability $`p(x_ix_j)=1/2`$, where $`x_j=x_{i\pm 1}`$. We use the same definition of the tunneling time to measure the behavior of this quantity. The tunneling time for the ideal case satisfies the simple power law as $`\tau (N_E)N_{E}^{}{}_{}{}^{\alpha }`$ and the exponent $`\alpha `$ is equal to 1. ($`\tau `$ is defined using the unit of sweep of $`N_E`$ sites.) Our simulational data for random walks in energy space yield a tunneling time that is well described by the power law $`\tau N^\alpha `$ as shown in Fig. 5(a). The solid lines in the graph have the simple power law as $`\tau (L)N_E`$ , and we see that our simulation result is very close to the ideal case. Since our method needs an extra effort to update the density of states to produce a flat histogram during the random walk in energy space, the tunneling time is much longer than the real space case. Also because the tunneling time depends on the accuracy of the density of states, which is constantly modified during the random walk in energy space, it is not a well defined quantity in our algorithm. The tunneling time shown in Fig. 5(a) is the overall tunneling time which includes all iterations with the modification factors from $`f_0=e^12.71828\mathrm{}`$ to the final modification factor $`f_{\text{final}}=\mathrm{exp}(10^8)1.00000001`$.

We should point out that the two processes are not exactly the same, since the random walk in real space uses the exact density of states($`g(x_i)=1`$). However the random walk in energy space requires knowledge of the density of states, which is a priori unknown. The algorithm we propose in this paper is a random walk with the density of states which is modified at each step during the walk in energy space. At the end of our random walk, the modification factor approaches 1, and the estimated density of states approaches the true value. The two processes are then almost identical.

Conventional Monte Carlo algorithms (such as the heat-bath algorithm) have an exponentially fast growing tunneling time. According to Berg’s study in reference , the tunneling time obeys the exponential law $`\tau (L)=1.46L^{2.15}e^{0.080L}`$. The Multicanonical simulational method has reduced the tunneling time from an exponential law to a power law as $`\tau (N_E)N_{E}^{}{}_{}{}^{\alpha }`$. However, the exponent $`\alpha `$ is as large as $`1.33`$ , which is far away from the ideal case $`\alpha =1`$. Only very recently, Janke and Kappler introduced the multibondic cluster algorithm, the exponent $`\alpha `$ is reduced to as small as $`1.05`$ for 2-dim ten state Potts model . In Fig. 5(b), we show our result with the same quantities obtained with multicanonical method and heat bath algorithm in reference . We should point out that just like multicanonical simulational method, our algorithm has a power increasing tunneling time with a smaller exponent $`\alpha `$. For small systems, our algorithm offers less advantage because of the effort needed to modify the density of states during the random walk. Very recently, Neuhaus has generalized this algorithm to estimate the canonical distribution for $`T<T_\text{c}`$, in magnetization space for the Ising model . He found that for small systems the exponent for CPU time versus L for our algorithm and Multicanonical Ensemble simulations are almost identical. However we see that our algorithm is very efficient for large systems, especially for $`L100`$. Our results in Fig. 5(b) are only for single range random walks, and multiple range random walks have been proven more efficient for larger systems.

## IV Application to a second order phase transition

The algorithm we proposed in this paper is very efficient for the study of any order phase transitions. Since our method is independent of temperature, it reduces the critical slowing down at the second-order phase transition $`T_\text{c}`$ and slow dynamics at low temperature. We estimate the density of states very accurately with a flat histogram, the algorithm will be a very efficient for general simulational problems by avoiding the need for multiple simulations at multiple temperatures.

To check the accuracy and convergence of our method, we apply it to the 2D Ising model with nearest neighbor interactions on a $`L\times L`$ square lattice. This model provides an ideal benchmark for new algorithms  and is also an ideal laboratory for testing theory . This model can be solved exactly, therefore we can compare our simulational results with exact solutions.

With the exact solution for the partition function on finite-size systems , and expansion of the expression by Mathematica, the density of states for the 2D Ising model on a square lattice can be obtained exactly . Beale  obtained the exact density of states up to $`L=32`$; and because of the memory and speed limitation of present computers, we could only extend the calculation up to $`L=50`$ with the Mathematica program provided by Beale. With the algorithm proposed in this paper, the density of states for the Ising model is estimated for $`L=50`$. The final modification factor for our random walk was $`1.000000001`$ for $`L=50`$. In Fig. 6(a) the simulational densities of states are compared with exact results. With the logarithmic scale used in the figure simulational data and exact solution overlap perfectly with each other for $`L=50`$. In the inset of the figure, we show the relative error $`\epsilon `$, which is generally defined by the ratio between error of the simulational data and exact values for any quantity $`X`$:

$$\epsilon (X)\frac{|X_{\text{sim}}X_{\text{exact}}|}{X_{\text{exact}}}$$
(13)
$`\epsilon (\mathrm{log}(g))`$ for most of the region is smaller than 0.1%. Such errors for low energy levels are directly related to the errors for the thermodynamic quantities calculated from the density of states. The average relative error is 0.019% for $`L=50`$. It is possible to estimate the density of states for small systems with the broad histogram method . Recent broad histogram simulational data  for the 2D Ising model on a $`32\times 32`$ lattice with $`10^6`$ MC sweeps yielded an average deviation of the microcanonical entropy of about 0.08 % from the exact solution . With our algorithm we obtain an average error as small as 0.035 % on the $`32\times 32`$ lattice with $`7\times 10^5`$ MC sweeps. Procedures for allowing $`f1`$ have been examined by Hüller  who used data from two densities of states for two different values of $`f`$ to extrapolate to $`f=1`$. However, his data for a small Ising system yield larger errors than our direct approach. The applicability of his method to large systems also needs a more detailed study.

The absolute density of states in Fig. 6(a) is obtained by the condition that the number of ground states is 2 for the 2D Ising model (all up or down). This condition guarantees the accuracy of the density of states at low energy levels which are very important in the calculation of thermodynamic quantities at low temperature. With this condition, when $`T=0`$, we can get exact solutions for internal energy, entropy and free energy when we calculate such quantities from the density of states. If we apply the condition that the total number of states is $`2^N`$ for the ferromagnetic Ising model, we can not guarantee the accuracy of the energy levels at or near ground states because the rescaled factor is dominated by the maximum density of states.

In Fig. 6(b), we show our estimation of the density of states of Ising model on $`256\times 256`$ lattice. Since the density of states for $`E>0`$ has almost no contribution to the canonical average at finite positive temperature, we only estimate the density of states in the region $`E/N[2,0.2]`$ out of the whole energy $`[2,2]`$. To speed up our calculation, we divide the desired energy region \[-2, 0.2\] into 15 energy segments, and estimate the density of states for each segment with independent random walks. The modification factor changes from $`f_0=e^12.71828\mathrm{}`$ to $`f_{\text{final}}=\mathrm{exp}(10^7)1.0000001\mathrm{}`$. The resultant density of states can be joined from adjacent energy segments. To reduce the boundary effects of the random walk on each segments, we keep about several hundred overlapping energy levels for random walks on two adjacent energy segments. The histograms of random walks are shown in the inset of this figure. We only require a flat histogram for each energy segment. To reduce the error of the density of states relevant to the accuracy of the thermodynamic quantities near $`T_\text{c}`$ we optimize the parameter and perform additional multiple random walks for the energy range $`E/N[1.8,1]`$ with same number of processors. For this we use the density of states obtained from the first simulations as starting points and continue the random walk with modification factors changing from $`\mathrm{exp}(10^6)1.000001`$ to $`\mathrm{exp}(10^9)1.000000001`$. The total computational effort is about $`9.2\times 10^6`$ visits on each energy levels. Note that the total number of possible energy levels is $`N1`$ and we perform random walks only on \[-2, 0.2\] out of \[-2, 2\]. The real simulational effort is about $`6.1\times 10^6`$ MC sweeps for the Ising model with $`L=256`$. With the program we implemented, it took about 240 CPU hours on a single IBM SP Power3 processor.

For $`L=256`$, we perform multiple random walks on different energy ranges, and one problem arises, that is the error of the density of states due to the random walk in a restricted energy range. Since the exact density of states for large systems is not available, we use $`L=32`$ to study such effects. We perform three independent random walks in the ranges $`E/N=[1.7,1.2]`$, $`E/N=[1.8,1.1]`$ and $`E/N=[1.9,1.0]`$ to calculate the densities of states on these ranges. In Fig. 7(a), we show the errors of our simulation results from the exact values. We make our simulational densities of states match up with the exact results at the left edges. It is very clear that the width of the energy range of the random walks is almost not relevant to the errors of the density of states. The reason is that the random walks only require the local histogram to be flat as we discussed in the previous section. But the errors for the last two densities on the right edges are significantly larger than the others. The problem is the boundary effect because in this simulation we treat the densities of states at edges as same as those away from edges. That is not correct. Another reason is due to the flat histogram and high density at right side edges. Our simulational effort for the densities at the right edges is not enough compared to those at the left side. This is also the reason why we have not seen big errors at left edges.

To study the influence of the errors of the densities of states on the thermodynamic quantities calculated from them, in the energy range that we perform random walks, we replace the exact density of states with the simulational density of states. In Fig. 7(b), the specific heat calculated from such density of states is showed as a function of temperature. We also show the exact value with the simulational data, the difference is obvious.

To reduce the boundary effect, we delete the last two density entries, and insert them into the exact density of states again, then the results in Fig. 7(c) are much better.

With our test in the three different ranges of energy, it is quite safe to conclude that boundary effect will not be present in our multiply random walks if we have a couple of energy levels overlap for adjacent energy ranges. In our real simulations for large systems, we have hundreds of overlaping energy levels.

Since the exact density of states is only available on small systems ($`L50`$), it is not so interesting to compare the simulational density of states itself. The most important thing is the accuracy of estimations of thermodynamic quantities calculated from such density of states on large systems. With the density of states on large systems, we apply canonical average formulas to calculate internal energy, specific heat, Gibbs free energy and entropy. Ferdinand and Fisher  obtained the exact solutions of above quantities for 2D Ising model on finite-size lattices. Our simulational results on finite-size lattice can be compared with those exact solutions.

In Fig. 8(a), we show the internal energy as a function of temperature. The internal energy is estimated from the canonical average over energy of the system as the equation (7). We also draw the exact solution in the same graph. The exact and simulational data perfectly overlap with each other in a wide temperature region from $`T=0`$ to $`T=8`$. Since no difference is visible with the scale used in Fig. 8(a), a more stringent test of the accuracy is provided by the inset which shows the relative errors $`\epsilon (U)`$. With the density of states obtained with our algorithm, the relative error of simulational internal energy for $`L=256`$ is smaller than 0.09% for the temperature region from $`T=0`$ to $`8`$. From eqn. (7) it is very clear that the canonical distribution serves as a weighting factor, and since the distribution is very narrow, $`U(T)`$ is only determined by a small portion of the density of states. (For the $`L=50`$ 2D Ising model at $`T_\text{c}`$, only the density of states for $`E/N[1.6,1.2]`$ contributes in a major way to the calculation.) Therefore the error $`\epsilon (U)`$ is also determined by the errors of the density of states in the same narrow energy range.

With the density of states and equation (8), we also can calculate the specific heat per lattice site as a function of temperature. Both simulational data and exact results near $`T_\text{c}`$ are shown in Fig. 8(b) for $`L=64`$, 128 and 256 Ising model in the vicinity of $`T_\text{c}`$. Within the resolution of the figure, we can see the difference between the simulational data and exact solutions. In the inset of Fig. 8(b), we present relative errors for our simulational data as a function of temperature for $`L=256`$. The errors for the specific heat for the Ising model on a $`256\times 256`$ lattice are smaller than $`4.5\%`$ in all temperature region $`T<8`$. Very recently, Wang, Tay and Swendsen  estimated the specific heat of the same model on a $`64\times 64`$ lattice by the transition matrix Monte Carlo re-weighting method , and for a simulation with $`2.5\times 10^7`$ MC sweeps, the maximum error in temperature region $`T[0,8]`$ was about 1%. When we apply our algorithm to the same model on the $`64\times 64`$ lattice, with a final modification factor of $`1.000000001`$ and a total of $`2\times 10^7`$ MC sweeps on single processor, the errors in the specific heat are reduced below 0.7% for all temperatures. The relatively large errors at low temperature reflect the small values for the specific heat at low temperature. According to the simulational data by broad histogram method, the errors in specific heat are very large even for systems as small as $`32\times 32`$ . Only very recently, they have reduced the error near $`T_\text{c}`$ to a small value for $`L=32`$ .

In Fig. 8(c), the free energy for $`L=256`$ is plotted as a function of temperature. For comparison the exact solution  is shown in the same figure. As expected, simulational and exact data overlap perfectly within the resolution of the figure. Since the system has a second-order phase transition, unlike the $`Q=10`$ Potts model, the first derivative of the free energy is a continuous function of temperature. The result very close to $`T_\text{c}(L=256)`$ is shown in the inset of Fig. 8(c) for $`L=256`$. The behavior of the free energy near $`T_\text{c}`$ is quite different from that of the first-order phase transition in Fig. 4(c). According to our calculation, the relative errors for all temperatures from $`T=0`$ to $`8`$ are smaller than 0.0008%, which means that our simulational data agree almost perfectly with the exact solution . Since we use the condition that the number of ground states is $`g(E=2N)=2`$ to normalize the density of states, the errors at low temperatures are extremely small.

The entropy of the 2D Ising model can be calculated with the equation(12). In Fig. 8(d), the simulational data and exact results are presented in the same figure. With the scale in the figure, the difference between our simulational data and exact solutions are not visible. In the inset of Fig. 8(d), the relative errors of our simulational data are plotted as a function of temperature. For the Ising model on a $`256\times 256`$ lattice, the relative errors are smaller than 1.2% for all temperature range. We also notice that the errors near $`T=0`$ decrease dramatically. The reason is the condition we use to normalize the density of states. Very recently, with the flat histogram method  and the broad histogram method , the entropy was estimated with $`10^7`$ MC sweeps for the same model on $`32\times 32`$ lattice; however, the errors in reference  are even much bigger than our errors for $`256\times 256`$!

## V Application to 3D $`\pm J`$ EA model

Spin glasses  are magnetic systems in which the interactions between the magnetic moments produce frustration because of some structural disorder. One of the simplest theoretical models for such systems is the Edwards-Anderson model  (EA model) proposed twenty five years ago. For such disordered systems, analytical methods can provide only very limited information, so computer simulations play a particularly important role. However, because of the rough energy landscape of such disordered systems, the relaxation times of the conventional Monte Carlo simulations are very long. The dynamical critical exponent was estimated as large as $`z6`$ . Normally simulations can be performed only on rather small systems, and many properties concerning the spin glasses are still left unclarified .

In this paper, we consider the three-dimensional $`\pm J`$ Ising spin glass EA model. The model is defined by the Hamiltonian

$$=\underset{i,j}{}J_{ij}\sigma _i\sigma _j$$
(14)
where $`\sigma `$ is an Ising spin and the coupling $`J_{ij}`$ is quenched to $`\pm 1`$ randomly. The summation runs over the nearest-neighbors $`i,j`$ on a simple cubic lattice.

One of most important issues for a spin glass model is the low temperature behavior. Because of the slow dynamics and rough phase space landscape of this model, it is also one of most difficult problems in simulational physics. The algorithm proposed here is not only very efficient in estimating the density of states but also very aggressive in finding the ground states. From a random walk in energy space, we can estimate the ground state energy and the density of states very easily. For a spin glass system, after we finish the random walk, we can obtain the absolute density of states by the condition that total number of states is $`2^N`$. The entropy at zero temperature can be calculated from either $`S_0=\mathrm{ln}(g(E_0))`$ or $`\underset{T0}{lim}\frac{UF}{T}`$, where $`E_0`$ is the energy at ground states. Both relations will give the same result since $`U`$ and $`F`$ are calculated from the same density of states. Our estimates for $`s_0=S_0/N`$ and $`e_0=E_0/N`$ per lattice site, listed in Table 2, agree with the corresponding estimates made with the multicanonical method. With our algorithm, we can estimate the density of states up to $`L=20`$ by a random walk in energy space for few hours on a 400MHz processor.

If we are only interested in the quantities directly related to the energy, such as free energy, entropy, internal energy and specific heat, one dimensional random walk in energy space will allow us to calculate these quantities with a high accuracy as we did in the 2D Ising model. However for spin glass systems, one of the most important quantities is the order parameter which can be defined by

$$q^{\text{E}\text{A}}(T)\underset{t\mathrm{}}{lim}\underset{N\mathrm{}}{lim}q(T,t),q(T,t)\underset{i=1}{\overset{N}{}}\sigma _i(0)\sigma _i(t)/N.$$
(15)
Here, $`N=L^3`$ is the total number of the spins in the system, $`L`$ is the linear size of the system, $`q(T,t)`$ is the auto-correlation function, which depends on the temperature $`T`$ and the evolution time $`t`$, and $`q(T,0)=1`$. When $`t\mathrm{}`$, $`q(T,t)`$ becomes the order parameter of the spin glass. This parameter takes the following values

$$q^{\text{EA}}(T)\{\begin{array}{c}=1\text{if}T=0\\ =0\text{if}TT_\text{g}\\ 0\text{if}0<T<T_{\text{ g}}\end{array},$$
(16)
The value at $`T=0`$ can be different from 1 in the case where the ground state is highly degenerate.

In our simulation, there is no temperature introduced during the random walk. And it is more efficient to perform a random walk in single system than two replicas. So the order-parameter can be defined

$$q\underset{i=1}{\overset{N}{}}\sigma _i^0\sigma _i/N.$$
(17)
where $`\{\sigma _i^0\}`$ is one of spin configurations at ground states and $`\{\sigma _i\}`$ is any configuration during the random walk. The behavior of $`q`$ we defined above is basically the same as the order-parameter defined by the Edwards and Anderson . It is not exact same order-parameter defined by Edwards and Anderson, but was used in the early numerical simulations by Binder et al..

After first generating a bond configuration we perform a one-dimensional random walk in energy space to find a spin configuration $`\{\sigma _i^0\}`$ for the ground states. Since the order-parameter is not directly related to the energy, to get a good estimate of this quantity we have to perform a two-dimensional random walk to obtain the density of states $`G(E,q)`$ with a flat histogram in $`E`$-$`q`$ space. This also allows us to overcome the barriers in parameter space (or configuration space) for such a complex system. The rule for the 2D random walk is the same as 1D random walk in the energy space.

In Fig. 9, we show the histogram of the 2D random walk in energy-order parameter space, which is very flat. With the density of states $`G(E,q)`$, we can calculate any quantities as we did in the previous sections. It is very interesting to study the roughness of this model. First we study the canonical distribution as a function of the order-parameter:

$$P(q,T)=\underset{E}{}G(E,q)e^{E/k_\text{B}T}$$
(18)
In Fig. 10(a), we show a 3D plot for the canonical distribution at different temperatures for one bond configuration of $`L=6`$ EA model. At low temperatures, there are four peaks, and the depth of the valleys between peaks depends upon temperature. When the temperature is high, the multiple peaks converge to a single central peak. Because we use the linear scale to show our result in the Fig. 10(a). It is not clear how deep the dips among peaks are. In Fig. 10(b), we show the canonical distribution using logarithmic scale for the same distribution but only at $`T=0.5`$, and we find that the dips are as deep as $`10^4`$. We also noted there actually are six peaks, but the plot with linear scale does not show all of them because two are as small as $`10^3`$ compared to other four peaks.

In Fig. 11(a), we show the roughness of the canonical distribution for another realization on an $`8^3`$ lattice. Because of the wide variation in the distribution at low temperature we used a logarithmic scale: the relative size of dips are as deep as $`10^{30}`$ at $`T=0.1`$. There are several local minima even at high temperatures. With conventional Monte Carlo simulations, it is almost impossible to overcome the barriers at the low temperature, so the simulation will get trapped in one of the local minima as shown in the figure. With our algorithm all states will be visited with more or less the same probability and trapping is not a problem.

With the density of states $`G(E,q)`$, we also can calculate the energy landscape by:

$$U(q,T)=\frac{\underset{E,q}{}EG(E,q)e^{\beta E}}{\underset{E,q}{}G(E,q)e^{\beta E}}$$
(19)
In Fig. 11(b), we show the internal energy as a function of order-parameter for temperatures $`T=0.12.0`$. We find that the landscape is a very rough at low temperatures. The roughness of the energy landscape agrees with the one for canonical distribution. But the maxima in energy landscape are corresponding to the minima approximately in the canonical distribution.

As we already noted in the previous paragraph, the roughness of the landscape of the spin glass model makes the conventional Monte Carlo simulation extremely difficult to apply. Therefore, even a quarter of a century after the model was proposed, we even can not conclude whether there is a finite phase transition between glass phase and disordered phase. With Monte Carlo simulations on a large system ($`64^2\times 128`$) and a finite-size scaling analysis on a small lattice, Marinari et al.  expressed doubt about the existence of the “well-established” finite-temperature phase transition of the 3D Ising spin glass . Their simulational data can be described equally well by a finite-temperature transition or by a $`T=0`$ singularity of an unusual type. Kawashima and Young’s simulational data could also not rule out the possibility of $`T_\text{g}=0`$ . Thus even the existence of the finite-temperature phase transition is still controversial, and thus the nature of the spin glass state is uncertain. Although the best available computer simulation results  have been interpreted as a mean-field like behavior with replica-symmetry breaking(RSB) , Moore et al. showed evidence for the droplet picture  of spin glasses within the Migdal-Kadanoff approximation. They argued that the failure to see droplet model behavior in Monte Carlo simulations was due to the fact that all existing simulations are done at temperatures too close to transition temperature so that system sizes larger than the correlation length were not used. As discussed in the previous paragraph, the lower the temperature is, the rougher the canonical distribution and energy landscape are; hence, it is almost impossible for conventional Monte Carlo methods to overcome the barrier between local minima and globe minima. It is possible to heat the system up to increase the possibility of escape from local minima by simulated annealing and the more recent simulated tempering method  and parallel tempering method , but it is still very difficult to perform equilibrium simulations at low temperatures. Very recently Hatano and Gubernatis proposed a bivariate multicanonical Monte Carlo method for the 3D $`\pm J`$ spin glass model, and their result also favors the droplet picture . Marinari, Parisi et al. argued, however, that the data were not thermalized . The nature of spin glasses thus remains controversial .

The algorithm proposed in this paper provides an alternative for the study of complex systems. Because we need to calculate the order-parameter with high accuracy, and this quantity is not directly related to the energy, we need to perform a random walk in the two-dimensional energy - order parameter space. After we estimate the density of states in this $`2D`$ space, we can calculate the order-parameter at any temperature from the canonical average. In 12(a), we show our results for the $`3D`$ EA model for $`L=4`$, 6 and 8. Because we need to perform a 2D random walk with a total of about $`L^6`$ states, the simulation is only a practical for a small system($`L8`$). The results in the figure are the average over 100 realizations for $`L=4`$, 50 realizations for $`L=6`$ and 20 for $`L=8`$.

We notice that the behavior of $`q(T)`$ is very similar to the magnetization (the order-parameter for the Ising model), but the finite value at low temperature is not necessarily equal to 1 because of the high degeneracy of the ground state for the spin glass model. The fluctuation of the order-parameter at the different temperatures for $`L=4`$, 6 and 8 is shown in the inset of the figure.

To estimate the transition temperature of the spin glass system, we calculated the fourth order cumulant as a function of temperature. In Fig. 12(b), we show our simulational results for $`L=4`$, 6 and 8. All curves clearly cross around $`T_\text{g}=1.2`$. Below this temperature, the spin configurations are frozen into some disorder ground states and the order parameter assumes a finite value. Above this temperature $`T_\text{g}`$, the system is in a disordered states and the order parameter vanishes.

One complication for simulation of such random systems is the determination of the relative importance of the error due to the simulation algorithm and the error due to the finite sampling of bond distributions. From Fig. 12 (a) and (b), we can not tell what the origin of the error bars is so we also performed multiple independent simulations for the same bond configuration on a $`L=6`$ 3D EA model. We found that the statistical errors for the order parameter and the fourth order cumulant from these simulations were much smaller than the error bars shown in Fig. 12 (a) and (b) for all temperatures. We conclude that the error bars in the figure arise almost completely from the randomness of the system.

The computational resources devoted here to the EA model were not immense. All our simulations for one bond configuration ($`L=4`$, 6 and 8) were performed within two days on (multiple) Linux machines ($`200800`$MHz) in the Center for Simulational Physics. This effort should thus be viewed as a feasibility study, and substantially more effort would be required to determine the nature of the spin glass phase or to estimate the transition temperature with high accuracy. Nonetheless, we believe that these results show the applicability of our method to systems with a rough landscape. Because of the number of states is about $`N^2`$ for $`2D`$ random walks, such calculations not only require huge memory during the simulation but also substantial disk space to store the density of states for the later calculation of thermodynamic quantities.

## VI Discussion and Conclusion

In this paper, we proposed an efficient algorithm to calculate the density of states directly for large systems. By modifying the estimate at each step of the random walk in energy space and carefully controlling the modification factor, we can determine the density of states very accurately. Using the density of states, we can then calculate thermodynamic quantities at essentially any temperature by applying simple statistical physics formulas. An important advantage of this approach is that we can also calculate the Gibbs free energy and entropy, quantities that are not directly available from conventional Monte Carlo simulations.

We applied our method to the $`2D`$ $`Q=10`$ Potts model which demonstrates a typical first-order phase transition. By estimating the density of states with lattices as large as $`200\times 200`$, we calculated the internal energy, specific heat, free energy and entropy in a wide temperature region. We found a typical first-order phase transition with a “discontinuity” for the internal energy and entropy at $`T_\text{c}`$. The first derivative of the free energy also shows such a discontinuity at $`T_\text{c}`$. The transition temperature estimated from simulational data is consistent with the exact solution.

We also applied our algorithm to the $`2D`$ Ising model, which shows a second-order phase transition. The density of states obtained by the end of our simulations was compared directly with the exact solution on $`50\times 50`$ lattice. The relative errors for most important energy levels are less than 0.019%. It was also possible to calculate the density of states for a $`256\times 256`$ lattice with a computational effort of $`6.1\times 10^6`$ Monte Carlo sweeps. With the accurate density of states, we calculated the internal energy, specific heat, Gibbs free energy and entropy. For all temperatures between $`T=0`$ and $`T=8`$, the relative errors are smaller than $`0.09\%`$ for internal energy, $`0.0008\%`$ for free energy, $`1.2\%`$ for entropy and $`4.5\%`$ for specific heat.

We should point out that our simulational results for $`L=256`$ are close to the same accuracy as the specific heat of the $`2D`$ Ising model for $`L=64`$ with the transition matrix Monte Carlo re-weighting method with $`2.5\times 10^7`$ MC sweeps . Our estimate of the entropy for the 2D Ising model for the $`L=256`$ lattice is even more accurate than the results obtained for the same quantity for $`L=32`$ with the broad histogram method and flat histogram method , which needed about $`10^7`$ MC sweeps. Our simulational effort for $`L=256`$ was $`6.1\times 10^6`$ MC sweeps.

The algorithm was also applied with success to the $`3D`$ $`\pm J`$ EA spin glass model for which we could determine the roughness of the energy landscape and canonical distribution in order-parameter space. The internal energy and entropy at zero temperature were estimated up to a lattice size $`20^3`$, and the transition temperature was estimated as about $`T_\text{g}=1.2`$.

In this paper, we only concentrated the random walk in energy space (and order-parameter space); however, the idea is very general and we can apply this algorithm to any parameters . The energy levels of the models treated here are perfectly discrete and the total number of possible energy levels is known before simulation, but in a general model such information is not available. Since the histogram of the random walk with our algorithm tends to be flat, it is very easy to probe all possible energies and monitor the histogram entry at each energy level. For some models where all possible energy levels can not be fitted in the computer memory or the energy is continuous, e.g. the Heisenberg model, we may need to discretize the energy levels. According our experience on discrete and continuous models, if the total number of possible energies is around the number of lattice sites $`N`$, the algorithm is very efficient for studying both first- or second-order phase transitions.

In this paper, we only applied our algorithm to simple models, but since the algorithm is very efficient even for large systems it should be very useful in the studies of general, complex systems with rough landscapes. It is clear, however, that more investigation is needed to better determine under which circumstances our method offers substantial advantage over other approaches and we wish to encourage the application of this approach to other models.

## ACKNOWLEDGMENTS

We would like to thank S. P. Lewis, H-B Schuttler, T. Neuhaus and A. Hüller for comments and suggestions, K. Binder, N. Hatano, P. M. C. de Oliveira and C. K. Hu for helpful discussions. We also thank M. Caplinger for support on technical matters and P. D. Beale for providing his Mathematica program for the calculation of the exact density of states for the $`2D`$ Ising model. The research project is supported by the National Science Foundation under Grant No. DMR-0094422.