text:
# Transition Matrix Monte Carlo Method

## 1 Introduction

The Monte Carlo technique  has served us well in the study of equilibrium statistical mechanics and other fields. The traditional local Monte Carlo method is simple, extremely general, and versatile. However, there are some intrinsic drawbacks. First, the convergence of the results to the exact values is slow. The basic probabilistic nature has limited the Monte Carlo error to decrease as $`1/\sqrt{t}`$, where $`t`$ is Monte Carlo steps or computer time. With exception of quasi Monte Carlo  for numerical integration, as long as we use a probabilistic approach, it does not appear possible to overcome this barrier. Most of the work to improve the efficiency of Monte Carlo method has been via variance reduction , which reduces the value of coefficient in front of the $`1/\sqrt{t}`$ law. Next, while the traditional Monte Carlo method is good for computing expectation values such as the internal energy and its derivatives, it is more difficult to compute the free energy or entropy .

Over the last few decades, a number of methods have been developed to compute the density of states. The histogram method  and the multiple histogram method  can be regarded from this point of view. The multi-canonical method  in some sense is also a computation of the density of states. Both of these methods involve the re-weighting of probabilities to construct the canonical distribution. Oliveira et al  proposed a broad histogram method, in which the density of states is also computed from simulation.

If the density of states can be computed with sufficient accuracy, then most thermodynamic quantities can be obtained with little further effort. This includes the moments of energy, entropy, and free energy. Moreover, the results are obtained as a continuous function of temperature from a single simulation. In this paper, we present such a method and study its efficiency. This method includes the use of a transition matrix , a stochastic matrix defined in the space of energy, and a class of related simulation algorithms . The present method has the elements of both the broad histogram method and multi-canonical method. The flat-histogram algorithm offers an effective way to compute density of states $`n(E)`$ for all energy $`E`$. With its multi-canonical element, it also offers fast dynamics for systems at first-order phase transitions. The use of transition matrix improves the efficiency of data analysis.

In the next section, we shall discuss the formalism and the essential aspects of the method. We also present the results of some numerical tests and discuss the connections of our method with previous methods. In the Appendices, we give details of a transition matrix dynamics.

## 2 Formalism

### 2.1 Markov Chain Monte Carlo

The Monte Carlo method aims at generating samples $`\sigma `$ with probability distribution $`P(\sigma )`$, where $`\sigma `$ is a particular state of the system. In the Ising model, which we shall use as a concrete example, $`\sigma `$ is a vector of all the spins $`\{\sigma _1,\sigma _2,\mathrm{},\sigma _N\}`$, where $`\sigma _i=\pm 1`$. In the usual application of the Monte Carlo method, the invariant distribution $`P(\sigma )`$ is given by the canonical distribution (Gibbs distribution) $`\mathrm{exp}\left(E(\sigma )/k_BT\right)`$. However, this need not be the case. In the equal-hit ensemble that we shall discuss later, $`P(\sigma )`$ is not known, and is not even unique. Nevertheless, it is still a valid Monte Carlo algorithm that can have significant advantages.

A sequence of states or samples is generated by a Markov chain  with transitions between states described by a matrix $`W(\sigma \sigma ^{})`$. This is the conditional probability that state moves to $`\sigma ^{}`$ given that the current state is $`\sigma `$. This matrix is known as a stochastic matrix and it must satisfy

$$\underset{\sigma ^{}}{}W(\sigma \sigma ^{})=1,W(\sigma \sigma ^{})0.$$
(1)
There is considerable freedom in choosing the matrix $`W`$, but the most important condition (or criterion) is detailed balance

$$P(\sigma )W(\sigma \sigma ^{})=P(\sigma ^{})W(\sigma ^{}\sigma ).$$
(2)
A Markov chain that satisfies the above condition is called a reversible Markov chain. This condition guarantees the invariance of the probability $`P(\sigma )`$ with respect to the transition matrix $`W`$, i.e.,

$$\underset{\sigma }{}P(\sigma )W(\sigma \sigma ^{})=P(\sigma ^{}).$$
(3)
Repeated application of $`W`$ to an arbitrary probability distribution makes the resulting probability distribution converge to a fixed point, $`P(\sigma )`$. We shall not elaborate on the condition that an invariant probability distribution exists and is unique. Roughly speaking, we must be able to make transitions in a finite numbers of steps from any initial state to any final state. This is known as ergodicity.

The standard Metropolis algorithm  is to take

$$W(\sigma \sigma ^{})=S(\sigma \sigma ^{})\mathrm{min}(1,\frac{P(\sigma ^{})}{P(\sigma )}),\sigma \sigma ^{},$$
(4)
where $`S(\sigma \sigma ^{})`$ is a selection function – a conditional probability of attempting to go to state $`\sigma ^{}`$ given that the current state is $`\sigma `$. Within the above formulation, it is required that the $`S`$ matrix is symmetric,

$$S(\sigma \sigma ^{})=S(\sigma ^{}\sigma ),$$
(5)
although this condition can also be relaxed . Following Oliveira , we call this condition microscopic reversibility. The diagonal elements of $`W`$ are fixed by the normalization condition, Eq. (1). Note that the diagonal elements are not needed explicitly in a computer simulation.

In a computer implementation, a move is selected according to $`S`$ (e.g., pick a site to flip a spin). The move is made if a random number $`\xi `$ between 0 and 1 is less than the flip rate $`\mathrm{min}(1,P(\sigma ^{})/P(\sigma ))`$; otherwise, it is rejected, and the original configuration $`\sigma `$ is counted once more as the next configuration in a Monte Carlo move.

Clearly the above formalism is very general. Although the procedure can be used to sample any distribution, it has its limitations. One drawback of standard algorithm is that the configurations generated are correlated. These correlations severely limit the efficiency of the method near phase transitions or for models with competing interactions and many local minima. A number of methods have been proposed to address this problem, such as the cluster algorithms , the multi-canonical methods , replica Monte Carlo , and simulated tempering . The flat histogram method presented in this paper is similar in some aspects to the multi-canonical method. The implementations of flat histogram method and transition matrix based methods are very simple and efficient.

### 2.2 Histogram

The concept of an energy histogram is essential to all of these methods. Other types of histogram of macroscopic quantities can be easily defined in analogy to the energy histogram and may be useful in some contexts, such as the joint histogram of energy and total magnetization. We define the energy histogram (in the case of a discrete energy spectrum) as the number of instances of each value of the energy $`E`$ generated during a Monte Carlo simulation; we denote the histogram by $`H(E)`$.

The histogram is important because of its direct relationship to the probability distribution of the energy in the system being simulated. If the probability of a state $`\sigma `$ is given by $`P(\sigma )`$ for a given simulation, then

$$h(E)=\underset{\sigma }{}\delta _{E(\sigma ),E}P(\sigma )=\underset{E(\sigma )=E}{}P(\sigma )$$
(6)
is the probability that the system has an energy $`E`$. If the Monte Carlo simulation in question generates $`m`$ configurations, then the expectation value of the histogram (average of the histogram over an infinite number of similar Monte Carlo runs) is given by

$$H(E)=mh(E).$$
(7)
For the canonical distribution, we have

$$P(\sigma )=f\left(E(\sigma )\right)/Z=Z^1\mathrm{exp}\left(E(\sigma )/k_BT\right),$$
(8)
where $`Z`$ is the partition function,

$$Z=\underset{\sigma }{}\mathrm{exp}\left(\frac{E(\sigma )}{k_BT}\right),$$
(9)
then $`h(E)=n(E)f(E)/Z`$. We define the density of states (for systems with discrete energy spectrum)

$$n(E)=\underset{E(\sigma )=E}{}1$$
(10)
as the number of states with energy $`E`$.

Note that the configuration dependence of the probability is only through energy implicitly. Thus, two configurations with the same energy will have the same probability. We shall call this the microcanonical property. The transition matrix Monte Carlo to be discussed below relies on this property crucially, while allowing the function $`f(E)`$ to be arbitrary.

The histogram $`H(E)`$ sampled during a Monte Carlo run (the number of visits to energy $`E`$) is an estimator to $`h(E)`$, i.e., $`H(E)h(E)`$. The usual canonical Monte Carlo method is equivalent to using the number of visits $`H(E)`$ to compute the moments of $`E`$ at the simulation temperature $`T_0`$. The histogram method of Ferrenberg and Swendsen  is based on the simple observation that density of states can be estimated (up to a proportionality constant) by $`n(E)H(E)/f(E)`$. With this information, the moments of $`E`$ can be extrapolated for nearby temperatures as well.

Clearly, if we can determine $`n(E)`$, then most of the energy related thermodynamic averages can be determined, such as internal energy, specific heat, free energy, and entropy. The free energy is given then by

$$F=k_BT\mathrm{ln}Z=k_BT\mathrm{ln}\underset{E}{}n(E)\mathrm{exp}(E/k_BT).$$
(11)
Other quantity of interest can also be computed if a “histogram” (as a function of $`E`$) of such quantity is also collected,

$$Q_T=\frac{{\displaystyle \underset{E}{}}Q_En(E)\mathrm{exp}(E/k_BT)}{{\displaystyle \underset{E}{}}n(E)\mathrm{exp}(E/k_BT)}\frac{{\displaystyle \underset{E}{}}Q_EH(E)}{{\displaystyle \underset{E}{}}H(E)}.$$
(12)
The main objective of this paper is to show that we can determine the density of states $`n(E)`$ for the whole range of energy with Monte Carlo sampling efficiently.

### 2.3 Detailed balance for histogram

The transition matrix defined below serves a dual purpose—for the computation of the density of states and for the construction of flat histogram algorithms. There are a number of ways to look at the transition matrix based methods. We shall take the detailed balance equation, (2), as a basic starting point. Consider all initial states $`\sigma `$ with energy $`E`$ and all final states $`\sigma ^{}`$ with energy $`E^{}`$. Each pair of states $`\{\sigma ,\sigma ^{}\}`$ has a detailed balance equation. Some of the equations may be the identity $`0=0`$ if the transition by a single-spin flip is not possible. Summing up the detailed balance equations for all the states $`\sigma `$ with a fixed energy $`E`$ and all the states $`\sigma ^{}`$ with a fixed energy $`E^{}`$, we have

$$\underset{E(\sigma )=E}{}\underset{E(\sigma ^{})=E^{}}{}P(\sigma )W(\sigma \sigma ^{})=\underset{E(\sigma )=E}{}\underset{E(\sigma ^{})=E^{}}{}P(\sigma ^{})W(\sigma ^{}\sigma ).$$
(13)
Assuming that the configuration probability distribution is a function of energy only, i.e., $`P(\sigma )f\left(E(\sigma )\right)`$, and defining the transition matrix in energy as

$$T(EE^{})=\frac{1}{n(E)}\underset{E(\sigma )=E}{}\underset{E(\sigma ^{})=E^{}}{}W(\sigma \sigma ^{}),$$
(14)
we have

$$n(E)f(E)T(EE^{})=n(E^{})f(E^{})T(E^{}E).$$
(15)
As a consequence of $`W`$ being a stochastic matrix, $`T(EE^{})`$ is also a stochastic matrix:

$$\underset{E^{}}{}T(EE^{})=1,T(EE^{})0,$$
(16)
with the histogram $`h(E)=n(E)f(E)/Z`$ being the invariant distribution:

$$\underset{E}{}h(E)T(EE^{})=h(E^{}).$$
(17)
Similar definition to Eq. (14) was introduced in ref.  in a different context.

### 2.4 Broad histogram equation

Because the matrix $`T(EE^{})`$ is composed of two factors, only the second of which depends on the specific ensemble under consideration, it is convenient to refer all calculations to the “infinite-temperature” transition matrix,

$$T_{\mathrm{}}(EE^{})=\frac{1}{N}N(\sigma ,E^{}E)_E,$$
(18)
where $`N`$ is the number of spins, or more generally, the number of allowed moves from a given state. If we define $`\mathrm{\Delta }E=E^{}E`$, we then have

$$\frac{N(\sigma ,\mathrm{\Delta }E)_E}{N}=\frac{1}{n(E)}\underset{E(\sigma )=E}{}\frac{N(\sigma ,\mathrm{\Delta }E)}{N}=\underset{E(\sigma )=E}{}\underset{E(\sigma ^{})=E^{}}{}\frac{S(\sigma \sigma ^{})}{n(E)}.$$
(19)
In a random single-spin-flip dynamics, $`S(\sigma \sigma ^{})`$ equals $`1/N`$ if the two configurations $`\sigma `$ and $`\sigma ^{}`$ differ by one spin, and zero otherwise. Thus, the second summation over $`\sigma ^{}`$ gives the number $`N(\sigma ,E^{}E)`$ of configurations of energy $`E^{}`$ that can be reached from $`\sigma `$ of energy $`E`$ by a spin flip. The first summation is over the configurations with energy $`E`$, i.e., a microcanonical average of the quantity $`N(\sigma ,\mathrm{\Delta }E)`$. The constancy of $`S(\sigma \sigma ^{})`$ for the nonzero matrix elements is important for this interpretation of $`N(\sigma ,\mathrm{\Delta }E)_E`$. The quantity $`N(\sigma ,\mathrm{\Delta }E)`$ is central to the current method, as well as to the broad histogram method .

Within the single-spin-flip dynamics, the matrix $`T`$ is then given by

$$T(EE^{})=T_{\mathrm{}}(EE^{})a(EE^{}),$$
(20)
where any flip rate $`a(EE^{})`$ can be inserted once $`T_{\mathrm{}}(EE^{})`$ has been determined. Substituting Eq. (20) into the energy detailed balance equation (15), we can cancel $`f(E)`$ and $`a(EE^{})`$ as for a valid dynamics which generates distribution $`P(\sigma )=f\left(E(\sigma )\right)/Z`$, we have the usual detailed balance,

$$\frac{a(EE^{})}{a(E^{}E)}=\frac{f(E^{})}{f(E)}.$$
(21)
The final equation is known as broad histogram equation, initially presented by Oliveira et al ,

$$n(E)N(\sigma ,E^{}E)_E=n(E^{})N(\sigma ^{},EE^{})_E^{}.$$
(22)
In terms of the transition matrix notation, this becomes

$$n(E)T_{\mathrm{}}(EE^{})=n(E^{})T_{\mathrm{}}(E^{}E).$$
(23)
The name “broad histogram” equation is historical and clearly a misnomer. The above equation has a very simple interpretation. Consider all pairs of states $`\sigma `$ with energy $`E`$ and states $`\sigma ^{}`$ with energy $`E^{}`$ such that the moves (or transitions) between $`\sigma `$ and $`\sigma ^{}`$ are allowed. These states correspond to states for which the matrix elements in $`S`$ are non-zero. Due to the microscopic reversibility, if $`\sigma `$ to $`\sigma ^{}`$ is allowed, so is the reverse move from $`\sigma ^{}`$ to $`\sigma `$. There are two ways to count the total number of moves, summing up from states with energy $`E`$ or summing up from states with energy $`E^{}`$. The state $`\sigma `$ has $`N(\sigma ,E^{}E)`$ ways to move into energy $`E^{}`$. The total number of moves to energy $`E^{}`$ from all states with energy $`E`$ is $`_{E(\sigma )=E}N(\sigma ,E^{}E)`$. By the reversibility requirement, we must have

$$\underset{E(\sigma )=E}{}N(\sigma ,E^{}E)=\underset{E(\sigma ^{})=E^{}}{}N(\sigma ^{},EE^{}).$$
(24)
A microcanonical average of any quantity is defined by

$$Q_E=\frac{{\displaystyle \underset{E(\sigma )=E}{}}Q(\sigma )}{{\displaystyle \underset{E(\sigma )=E}{}}1}=\frac{1}{n(E)}\underset{E(\sigma )=E}{}Q(\sigma ).$$
(25)
Using this definition, the previous equation (24) is reduced to Eq. (22). This argument is first put forth by Oliveira  and by Berg and Hansmann . Clearly, the result does not depend on what type of moves we use, so long as it satisfies the reversibility condition.

### 2.5 $`\mathrm{𝐓𝐓𝐓}`$ identity

The detailed balance condition imposes a restriction on the transition matrix, which we call the $`TTT`$ identity. Consider three distinct energy levels, $`E`$, $`E^{}`$, and $`E^{\prime \prime }`$, for which energy transition matrix elements among them are nonzero. Let us write down three equations associated with the transitions among them:

$`h(E)T(EE^{})`$ $`=`$ $`h(E^{})T(E^{}E),`$ (26)
$`h(E^{})T(E^{}E^{\prime \prime })`$ $`=`$ $`h(E^{\prime \prime })T(E^{\prime \prime }E^{}),`$ (27)
$`h(E^{\prime \prime })T(E^{\prime \prime }E)`$ $`=`$ $`h(E)T(EE^{\prime \prime }).`$ (28)
Multiplying the left and right sides of the three equations together and assuming that the product $`h(E)h(E^{})h(E^{\prime \prime })`$ is nonzero, we can cancel this factor from the equation and obtain:

$$T(EE^{})T(E^{}E^{\prime \prime })T(E^{\prime \prime }E)=T(E^{}E)T(EE^{\prime \prime })T(E^{\prime \prime }E^{}).$$
(29)
This is the $`TTT`$ identity. The importance of this equation is that it does not require the knowledge of the stationary distribution to check for agreement of the data with the condition of detailed balance. While for normal Monte Carlo simulation, the detailed balance is built-in directly to the transition matrix $`W`$, this is not the case for some of the transition matrix Monte Carlo schemes that have been proposed. One implication of detailed balance violation is that the microcanonical property that all states with the same energy have the same probability is violated. This detailed balance violation for the initial choice of Oliveira’s broad histogram dynamics, a particular choice of the transition rate $`W`$, has been demonstrated explicitly for small systems .

The significance of the $`TTT`$ identity is that given the probability $`h`$ of energy having value $`E`$, if we can predict $`h^{\prime \prime }`$ at energy $`E^{\prime \prime }`$ by the detailed balance equations in two ways, one directly from $`E`$ to $`E^{\prime \prime }`$, one by two hops, from $`E`$ to $`E^{}`$, and then $`E^{}`$ to $`E^{\prime \prime }`$, then the $`TTT`$ identity guarantees that the results are exactly the same. That is,

$$h^{\prime \prime }=h\frac{T(EE^{\prime \prime })}{T(E^{\prime \prime }E)},$$
(30)
and

$$\stackrel{~}{h}^{\prime \prime }=h^{}\frac{T(E^{}E^{\prime \prime })}{T(E^{\prime \prime }E^{})},h^{}=h\frac{T(EE^{})}{T(E^{}E)}.$$
(31)
The $`TTT`$ identity says that the two predictions based on the detailed balance are equal, $`h^{\prime \prime }=\stackrel{~}{h}^{\prime \prime }`$.

Detailed balance implies $`TTT`$ identity. Is the reverse also true? I.e., given a complete set of $`TTT`$ identities, do they imply detailed balance in the sense of equation $`h(E)T(EE^{})=h(E^{})T(E^{}E)`$ for all $`E`$ and $`E^{}`$? The answer is yes. The $`TTT`$ identity turns out to guarantee a consistent (detailed balance) solution involving three jumps, say, $`E`$ to $`E^{}`$, to $`E^{\prime \prime }`$, to $`E^{\prime \prime \prime }`$, versus $`E`$ to $`E^{\prime \prime \prime }`$ directly when such jumps are allowed. Therefore, $`TTTT`$ identities that follow from detailed balance are automatically fulfilled and neither provide further information, nor require separate proof. Naturally, when we consider adding more complex Monte Carlo moves, either to improve the efficiency of a calculation or to reflect the nature of a more complex model, more $`TTT`$ identities are generated and identities with more factors of $`T`$ are automatically satisfied.

We define a quantity to measure the detailed balance violation for three energy levels for which the transitions among them are nonzero as

$$v=\left|1\frac{\widehat{T}(EE^{})\widehat{T}(E^{}E^{\prime \prime })\widehat{T}(E^{\prime \prime }E)}{\widehat{T}(E^{}E)\widehat{T}(EE^{\prime \prime })\widehat{T}(E^{\prime \prime }E^{})}\right|,$$
(32)
where $`\widehat{T}(\mathrm{})`$ is Monte Carlo estimate of $`T(\mathrm{})`$. For a single-spin-flip dynamics with the Metropolis rate, the energy transition matrix is given by Eq. (20), with $`a(EE^{})=\mathrm{min}(1,f(E^{})/f(E))`$, thus the above equation is equivalent to

$$v=\left|1\frac{N(\sigma ,E^{}E)_EN(\sigma ^{},E^{\prime \prime }E^{})_E^{}N(\sigma ^{\prime \prime },EE^{\prime \prime })_{E^{\prime \prime }}}{N(\sigma ^{},EE^{})_E^{}N(\sigma ,E^{\prime \prime }E)_EN(\sigma ^{\prime \prime },E^{}E^{\prime \prime })_{E^{\prime \prime }}}\right|.$$
(33)
For Ising model where energies are equally spaced, we consider three levels at $`E`$, $`E+4J`$ and $`E+8J`$, where $`J`$ is the coupling constant. A plot of $`v`$ for Ising model is presented in ref. .

### 2.6 Flat histogram dynamics

If a configuration has probability $`P(\sigma )=f\left(E(\sigma )\right)/Z`$, then the histogram is $`H(E)n(E)f(E)`$. A flat histogram is obtained if we take $`f(E)1/n(E)`$. A single spin flip with the flip rate $`\mathrm{min}(1,n(E)/n(E^{}))`$ can be used to do the simulation. Lee’s reformulation of multi-canonical method  is essentially this. The trick there is to determine $`n(E)`$ efficiently .

From the equation describing detailed balance for the transition matrix, (23), we can write the acceptance rate as

$$a(EE^{})=\mathrm{min}(1,\frac{T_{\mathrm{}}(E^{}E)}{T_{\mathrm{}}(EE^{})}).$$
(34)
This is the first equation derived for flat histogram dynamics, although we will show below that it is not unique. This rate is first proposed in ref. , and is independently discovered by Li . Unlike the quantity $`n(E)`$, a good approximation is already available in the very beginning, since we can use the instantaneous value $`N(\sigma ,E^{}E)`$ as a preliminary estimate for $`T_{\mathrm{}}(EE^{})`$. A cumulative average of contributions to $`T_{\mathrm{}}(EE^{})`$ can be used as a convenient, and remarkably accurate, approximation for the microcanonical average. We shall discuss how good this scheme is in a later section.

There is another equivalent way of looking at the problem. Consider the energy detailed balance equation in the form

$$h(E)T_{\mathrm{}}(EE^{})a(EE^{})=h(E^{})T_{\mathrm{}}(E^{}E)a(E^{}E).$$
(35)
If we require that the histogram is a constant, $`h(E)=h(E^{})=\mathrm{const}`$, then the spin-flip rate must satisfy the following equation,

$$\frac{a(EE^{})}{a(E^{}E)}=\frac{T_{\mathrm{}}(E^{}E)}{T_{\mathrm{}}(EE^{})}.$$
(36)
Clearly, Eq. (34) satisfies the above equation. Moreover, there is a whole family of choices of the transition rates. Some of the choices are given in Table 1.

### 2.7 $`𝐍`$-fold way

The standard Metropolis algorithm contains two steps. First, a move is proposed. Next, this move is accepted with probability $`a(EE^{})`$ or rejected with probability $`1a(EE^{})`$, where $`0a(EE^{})1`$. Is it possible to always accept a move without sampling the same configuration repeatedly? The answer is yes, if we are willing to keep extra information during the simulation. In the flat histogram method, this extra information is already there for free. It is precisely $`N(\sigma ,\mathrm{\Delta }E)`$.

In an $`N`$-fold way simulation (also known as event-driven simulation) , we do not change the dynamics; it is fully equivalent to the usual single-spin-flip dynamics. However, there is a substantial improvement in efficiency in those situations for which the rejection rate is high. The method begins by computing the total probability that a move would be accepted with the standard approach. For the Ising model, this probability is

$$A=A(\sigma )=\underset{i=1}{\overset{N}{}}\frac{1}{N}a\left(E(\sigma )E(\overline{\sigma }^i)\right),$$
(37)
where the factor $`1/N`$ is due to the fact that each spin located at site $`i`$ is picked up with probability $`1/N`$. The notation $`\overline{\sigma }^i`$ refers to a configuration with the spin at site $`i`$ reversed in sign. The quantity $`A`$ is the probability that any spin is flipped. Since the flip rate depends on the initial and final energies only, we can simplify the above equation as

$$A=\underset{\mathrm{\Delta }E}{}\frac{N(\sigma ,\mathrm{\Delta }E)}{N}a(EE+\mathrm{\Delta }E).$$
(38)
We divide the possible moves into classes according to their energy increment $`\mathrm{\Delta }E`$. Within a given class, each spin has the same flipping probability. We now set the probability of the class with energy change $`\mathrm{\Delta }E`$ being chosen as

$$P(\mathrm{\Delta }E)=\frac{1}{AN}N(\sigma ,\mathrm{\Delta }E)a(EE+\mathrm{\Delta }E).$$
(39)
A spin in this class is then chosen at random and flipped with probability one. As a practical consideration in designing algorithms, we note that the condition that $`a(EE^{})`$ must be between 0 and 1 can now be relaxed because of the normalization by $`1/A`$ in this equation.

In the original algorithm , Monte Carlo time is rescaled to make the dynamics equivalent to that of the usual algorithm. For the purpose of calculating averages, we will reweight each configuration to achieve the same effect. Each configuration $`\sigma `$ in the original single spin flip has equal weight. The Monte Carlo average of a quantity $`Q`$ is computed as

$$Q=\frac{1}{m}\underset{t=1}{\overset{m}{}}Q_t.$$
(40)
Here $`Q_t`$ at step $`t`$ and the subsequent steps could be the same, due to the possibility of rejecting a move. The probability that a move is rejected is $`R=1A`$. Thus the life-time of a configuration has the probability distribution

$$P_t=(1R)R^{t1},t=1,2,3,\mathrm{}.$$
(41)
The average life-time for configuration $`\sigma `$ is then

$$\overline{t}=\underset{t=1}{\overset{\mathrm{}}{}}tP_t=\frac{1}{1R}=\frac{1}{A}.$$
(42)
In the $`N`$-fold-way simulation, each step generates a distinct state from the immediate preceding one; each of these states is supposed to last for a during of $`\overline{t}`$, on average. Thus in replacing Eq. (40) the correct formula for statistical average with $`N`$-fold-way simulation is

$$Q=\underset{t=1}{\overset{m}{}}\frac{Q_t}{A}/\underset{t=1}{\overset{m}{}}\frac{1}{A}.$$
(43)
Naively the computation of $`N(\sigma ,\mathrm{\Delta }E)`$ seems to require $`O(N)`$ basic steps at each single-spin-flip move. But the effect of changing configuration by a flip is local, involving only the site in question and its neighboring sites. We only need to compute the changes in $`N(\sigma ,\mathrm{\Delta }E)`$. Thus each move takes $`O(1)`$ in computer time. It is few times slower than a corresponding straightforward single-spin-flip program. The $`N`$-fold way does require extra memory, as a list of spins for each class is required in order to be able to pick a spin from the class with a computer time of $`O(1)`$.

### 2.8 Equal-hit algorithms

The flat histogram ensemble in some sense is the best ensemble to evaluate the density of states, for each energy level is sampled with the same frequency. However, as we have seen in the $`N`$-fold way, this is not entirely true, as some configurations are weighted more than others. The equal-hit algorithms  generate fresh configurations with equal probability at each energy. There is a very interesting aspect of these algorithms—the histograms in such algorithms are not unique and depend on the details of the dynamics.

The energy histogram in the normal single-spin-flip dynamics (that is, not in the $`N`$-fold way) is computed as

$$h(E)=\delta _{E(\sigma ),E},$$
(44)
i.e., the contribution to the histogram from the configuration $`\sigma `$ of energy $`E(\sigma )`$ is 1 for $`E=E(\sigma )`$ and zero for other energies. The angular brackets refer to Monte Carlo average. In the $`N`$-fold way, the statistics are weighted with $`1/A`$ to get the equivalent result in the original dynamics, so we have

$$h(E)=\delta _{E(\sigma ),E}\frac{1}{A(\sigma )}_\mathrm{N}/\frac{1}{A(\sigma )}_\mathrm{N}.$$
(45)
The angular brackets indicate simple arithmetic average over the samples generated in an $`N`$-fold way, c.f. Eq. (43). We put a subscript N to emphasize the fact that the average is over $`N`$-fold way samples. We define the hits as

$$u(E)=\delta _{E(\sigma ),E}_\mathrm{N}.$$
(46)
This quantity measures the average number of fresh configurations generated at each energy, since in the $`N`$-fold way, each configuration $`\sigma `$ in the sample is distinct from the previous sample due to the fact that there is no rejection in $`N`$-fold way. We can relate the histogram to the hits by

$$h(E)=u(E)1/A_{E,\mathrm{N}}/1/A_\mathrm{N},$$
(47)
where

$$1/A_{E,\mathrm{N}}=\frac{\delta _{E(\sigma ),E}\frac{1}{A(\sigma )}_\mathrm{N}}{\delta _{E(\sigma ),E}_\mathrm{N}}=A_E^1$$
(48)
is the average of inverse acceptance rate for states of fixed energy in the $`N`$-fold way samples. Note that this is not the same as the microcanonical average of $`1/A`$, which in general is

$$Q(\sigma )_E=\frac{Q(\sigma )/A_{E,\mathrm{N}}}{1/A_{E,\mathrm{N}}}.$$
(49)
Substituting Eq. (47) into energy detailed balance equation, (35), we have

$`u(E)1/A_{E,\mathrm{N}}T_{\mathrm{}}(EE^{})a(EE^{})=`$
$`u(E^{})1/A_{E^{},\mathrm{N}}T_{\mathrm{}}(E^{}E)a(E^{}E).`$ (50)
Equal-hit is realized if the flip rate satisfies the above equation with $`u(E)=u(E^{})=\mathrm{const}`$. For example,

$$a(EE^{})=\mathrm{min}(1,\frac{1/A_{E^{},\mathrm{N}}T_{\mathrm{}}(E^{}E)}{1/A_{E,\mathrm{N}}T_{\mathrm{}}(EE^{})}).$$
(51)
Other choices are also possible; some of them are given in Table 1.

We note that since the transition rate depends on the underlying dynamics through $`A`$, there is no guarantee that the histogram $`h(E)`$ is unique. In fact, in the equal-hit dynamics, $`h(E)`$ is not known, and has to be determined self-consistently through the equal-hit dynamics. We note that $`A`$ is defined by Eq.(37) using $`a`$, which uses $`A`$ in Eq. (51). There is also a peculiarity that $`h(E)`$ diverges for some choices of the transition rates for the ground states.

## 3 Numerical Tests

In this section, we evaluate the performance of various proposed algorithms. In Table 1, we list a dozen possible choices of flip rates, $`a(EE^{})`$. For each rate, the simulation can be implemented with or without the $`N`$-fold way.

As the flip rates require the knowledge of exact microcanonical average, which is not available, we use the crucial approximation

$$T_{\mathrm{}}(EE+\mathrm{\Delta }E)\frac{1}{H(E)N}\underset{t=1}{\overset{m}{}}\delta _{E(\sigma ^t),E}N(\sigma ^t,\mathrm{\Delta }E),$$
(52)
where energy histogram $`H(E)`$ is the number of samples generated for a given energy, $`H(E)=_{t=1}^m\delta _{E(\sigma ^t),E}`$, $`m`$ is the total number of samples generated so far, and $`\sigma ^t`$ is the configuration at step $`t`$ in the algorithm. We collect a sample after every attempt of moves. The above expression is suitable for the normal single-spin-flip algorithm. In $`N`$-fold way, the statistics have to be weighted by $`1/A(\sigma )`$. Whenever information is not available, we set the flip rate to 1. This biases the system to visit unexplored energy levels.

In the above method, the simulation is started automatically, without an iteration process. This bootstrap is efficient and is also correct in the sense that it converges. We shall call the above method iteration 0. Strictly speaking, iteration 0 need not be a valid Monte Carlo algorithm as the transition rates are fluctuating quantities, thus the normal Markov chain theory for convergence does not apply. However, numerical results do support convergence, although a rigorous proof is lacking.

### 3.1 Histograms

Figure 1 presents energy and hit histograms for a $`16\times 16`$ square lattice for an 8-state ferromagnetic Potts model. As expected, the energy histogram is flat for the flat histogram algorithm. The hits are roughly proportional to the number of different configurations generated. This is not exactly true since in the $`N`$-fold way, we can only guarantee that the next configuration is different from the immediate preceding one. The new configuration can be the same with configurations in earlier steps. By requiring that the hits are equal for all energies, we obtain the equal-hit algorithm with the corresponding energy and hit distribution shown in the lower part of the figure.

Due to the statistical nature of the histogram and also due to the fact that the energy range is explored similarly to a random walk, the histogram is not exactly flat, but fluctuates around a mean. During a simulation of length $`t`$ of Monte Carlo steps (sweeps), we have generated $`t/\tau `$ independent samples, where $`\tau `$ is correlation time related to the histogram. These samples are distributed to order $`N`$ bins of different energies (for the two-dimensional Ising model, it is exactly $`N1`$ bins). Thus, each bin has about $`t/(\tau N)`$ independent samples. The relative fluctuation of the histogram is (asymptotically for large $`t`$)

$$\frac{\delta H(E)}{H(E)}\sqrt{\frac{\tau N}{t}}.$$
(53)
Although the above argument applies for the fluctuation between different runs, it is also reasonable to apply it to the fluctuation among different energies, since the samples between different energies are assumed independent. The same analysis also applies to the hits in equal-hit algorithm.

Equation (53) can be used as a definition for the correlation time $`\tau `$. A perfect Poisson process has a constant correlation time ($`\tau =1`$); an ideal random walk in energy is $`\tau N`$. We compute $`\tau `$ for the nearest neighbor Ising model on one-dimensional chain, two-dimensional square lattice, and three-dimensional cubic lattice. We found that numerically Eq. (53) is approximately satisfied, with the correlation time growing with size, $`\tau L`$ in one dimension, $`\tau L^{0.7}`$ in two dimensions, and $`\tau L^{1.2}`$ or $`L^2`$ in three dimensions, see Fig. (2). The algorithm 0 becomes inefficient (slow convergence) for large lattices in two and three dimensions. We shall discuss this further later on.

### 3.2 Rate of convergence

We tested the zeroth iteration algorithms for convergence to the exact values for the infinite temperature transition matrix, $`T_{\mathrm{}}(EE^{})=\frac{1}{N}N(\sigma ,E^{}E)_E`$, on a $`5\times 5`$ Ising square lattice model. The exact result $`T_{\mathrm{}}`$ is obtained by an exhaustive enumeration. Figure 3 is a plot of overall error in the $`N`$-fold-way simulation data $`\widehat{T}`$, defined by

$$ϵ_T=\underset{E,E^{}}{}\left|\widehat{T}_{\mathrm{}}(EE^{})T_{\mathrm{}}(EE^{})\right|,$$
(54)
versus Monte Carlo simulation length $`t`$ (averaged over many runs). The Monte Carlo times are in units of sweeps ($`N`$ moves). The asymptotic value $`ϵ_T\sqrt{t}`$ for large $`t`$ is listed in Table 1, which characterizes the rate of convergence.

As expected, the best algorithm is algorithm number $`1`$, when the input $`n(E)`$ is exact. The algorithms (or rates) number 0, 1, and 7 are only slightly worse than the best. Numbers 5 and 6 are the second best. The rate numbers 3, 4, and 11 may still converge extremely slowly or not converge at all. We do not understand why this is so. Number 9 does not converge, which we know since it does not satisfy detailed balance equation. Using instantaneous values rather than average values is clearly wrong.

Errors for single-spin-flip algorithm are about 1.5 larger than with the $`N`$-fold-way, or about $`1.5^22.3`$ less efficient in Monte Carlo steps. However, since the $`N`$-fold way is slower in CPU time than single-spin flips by a factor of 2 or more, the two methods are comparable in overall efficiency in this particular instance.

From the above results we conclude that even the zeroth iteration converges to the correct results. This does not mean that the rate of convergence is uniform in $`E`$. In fact, we found for large lattices, the violation of detailed balance, $`v`$, is large at the two ends of the energy spectrum .

### 3.3 Determining the density of states

There are a number of different ways of determining the density of states. The matrix $`T_{\mathrm{}}(EE^{})`$ has eigenvalue 1 with corresponding left eigenvector $`n(E)`$. However, the solutions of the eigenvalue problem are numerically unstable. The broad histogram equation, (22), can be used. In the simplest application, we can ignore the extra equations, and consider only these with smallest $`\mathrm{\Delta }E=E^{}E`$, and obtain solution from iteration,

$$\mathrm{ln}n(E^{})=\mathrm{ln}n(E)+\mathrm{ln}\frac{T_{\mathrm{}}(EE^{})}{T_{\mathrm{}}(E^{}E)}.$$
(55)
Since there are more equations of this type than the unknown $`n(E)`$, we can use least-squares method. Our experience suggests that we should view the problem as an optimization with nonlinear constraints.

There are two possible models for an optimization solution. Introducing the optimization variable $`S(E)=\mathrm{ln}n(E)`$, consider,

$$\mathrm{minimize}\underset{E,E^{}}{}\frac{1}{\sigma _{E,E^{}}^2}\left(S(E^{})S(E)\mathrm{ln}\frac{T_{\mathrm{}}(EE^{})}{T_{\mathrm{}}(E^{}E)}\right)^2,$$
(56)
subject to any known constraints. For the $`d`$-dimensional Ising model, we have

$$S(E)=S(E),S(dJN)=\mathrm{ln}2,\underset{E}{}\mathrm{exp}\left(S(E)\right)=2^N.$$
(57)
The three conditions are the symmetry between low and high energies, the degeneracy for the ground states, and the total number of states. The weight $`\sigma ^2(E,E^{})`$ is the variance of the Monte Carlo estimates of the quantity $`\mathrm{ln}\left(T_{\mathrm{}}(EE^{})/T_{\mathrm{}}(E^{}E)\right)`$. The above minimization problem is essentially a linear problem. We solve it by an iterative method.

The second, different formulation of the optimization is expressed in variable $`T_{\mathrm{}}(EE^{})`$,

$$\mathrm{minimize}\underset{E,E^{}}{}\frac{1}{\sigma _{E,E^{}}^2}\left(T_{\mathrm{}}(EE^{})\widehat{T}_{\mathrm{}}(EE^{})\right)^2,$$
(58)
where $`\widehat{T}_{\mathrm{}}(EE^{})`$ is Monte Carlo estimate with error $`\sigma _{E,E^{}}`$, and $`T_{\mathrm{}}(EE^{})`$ is unknown. The minimization is subject to the conditions

$`0T_{\mathrm{}}(EE^{})1,`$ $`{\displaystyle \underset{E^{}}{}}T_{\mathrm{}}(EE^{})=1,`$
$`{\displaystyle T_{\mathrm{}}(EE^{})}`$ $`=`$ $`{\displaystyle T_{\mathrm{}}(E^{}E)}.`$ (59)
The last one is symbolically a $`TTT`$ identity, see Eq. (29) for a concrete example. For the Ising model, there is also an additional symmetry relation

$$T_{\mathrm{}}(EE+\mathrm{\Delta }E)=T_{\mathrm{}}(EE\mathrm{\Delta }E).$$
(60)
This set of constraints is much more difficult to handle. Intuitively, this second optimization problem should give better result than the first one. However, this is not the case, at least for the two-dimensional Ising model. For this model, optimization in $`T_{\mathrm{}}(\mathrm{})`$, Eq. (58), gives twice the error ($`ϵ_n`$ defined below) of the first optimization method. A simple iteration with Eq. (55) and $`\mathrm{\Delta }E=\pm 4J`$ has 4 times the error comparing to optimizing solution of Eq. (56) and (57).

Figure 4 is another convergence test plot for algorithm 0 and a two-stage simulation, both with $`N`$-fold way, on the two-dimensional Ising square lattices. In this plot, we consider the relative error per energy level for density of states,

$$ϵ_n=\frac{1}{N1}\underset{E}{}\left|\frac{\widehat{n}(E)}{n(E)}1\right|\left|\widehat{S}(E)S(E)\right|.$$
(61)
The normalization by $`N1`$ is to exclude ground state energy and its symmetric state energy, for which the exact values are imposed. The exact value $`n(E)`$ is obtained according to ref. ; $`\widehat{n}(E)`$ is Monte Carlo estimate obtained from solving Eq. (56) with $`\sigma _{E,E^{}}=1`$.

We see signs of slower convergence for large lattice sizes for algorithm 0. A two-stage simulation improves the efficiency to that of using the exact density of states in flip rates. We first apply the algorithm 0 using cumulative average in the flip rates; we then apply algorithm $`1`$ (Lee’s method of multi-canonical algorithm), using the density of states obtained in the first step. Both steps use the same number of Monte Carlo sweeps. With the two-stage algorithm, the slowing down is roughly $`\tau _nL^{1.5}`$, using a similar definition as given in Eq. (53).

### 3.4 Dynamic characteristics of the algorithms

We have already presented two correlation times — $`\tau `$ characterizes the convergence of the histogram and $`\tau _n`$ characterizes the convergence of the density of states. Both of them showed the effect of reduced performance when the size is increased. By definition of the correlation times $`\tau `$ and $`\tau _n`$, they also measure roughly how many Monte Carlo steps are needed to generate independent samples. Of course, in $`\tau _n`$, it also reflects the effect of data analysis methods.

Another measure of the performance of algorithms is the tunneling time. In our study, we define the tunneling time as the average Monte Carlo steps that the system in the lowest energy state goes to the highest energy state, or vice versa. More precisely, as soon as a ground state is reached, we record the current Monte Carlo move $`t_1`$, and then look for the highest energy, which may happen at $`t_2`$. We then look for the ground state again. The difference $`(t_2t_1)/N`$ consists of one sample for the tunneling time.

The tunneling time for the two-dimensional Ising model is very well described by $`\tau _t0.4L^{2.8}`$ (in units of sweeps). An ideal random walk in the space of energy would have an exponent of 2 ($`\tau _tN`$ in general). The dynamics is close but not quite random walk in energy.

The spin glass is one of the most difficult systems to simulate. The performance of the flat histogram algorithm for the two-dimensional $`\pm J`$ spin glass is presented in ref. . In Fig. 5, we show the tunneling time as a function of system linear size $`L`$ for three algorithms (algorithm 0 and 1 with $`N`$-fold way, and algorithm 0 without $`N`$-fold way) on the three-dimensional Ising spin glass. The gain from $`N`$-fold way as comparing to standard single-spin-flip is by a constant factor. The equal-hit algorithm is about a factor of 6 faster in tunneling times than the algorithm 0 without $`N`$-fold way. Unfortunately, this gain is not very significant as the $`N`$-fold-way program runs few times slower than standard single-spin flips. The slowing down exponent is about 8, this is comparable to that of multi-canonical method .

## 4 Energy Transition Matrix Dynamics

The stochastic matrix $`W`$ describes Monte Carlo dynamics in the space of spin configurations. Such state space is extremely large, containing $`2^N`$ states, from which Monte Carlo moves sample only a very small fraction. On the other hand, we introduced a new stochastic matrix $`T`$ in a coarse-grained space of energy. The matrix $`W`$ and $`T`$ are related by Eq. (14). We shall call the dynamics induced by this stochastic matrix energy transition matrix dynamics. In discrete time step $`i`$, the dynamics describes the evolution of the histogram $`h_i(E)`$ as

$$h_{i+1}(E)=\underset{E^{}}{}h_i(E^{})T(E^{}E).$$
(62)
What is the significance of this dynamics? The dynamics describes the change of energy distribution through the following single-spin-flip moves: given the current state $`\sigma `$ with energy $`E`$, pick a new state $`\sigma ^{}`$ with the same energy $`E`$ among all the $`n(E)`$ degenerate spin states with equal probability, flip a spin according to the usual spin flip rate as embedded in $`W`$. As we can see, since the state changes at random to a completely new state of the same energy to try another flip, its dynamics is substantially faster than single spin flip or even cluster flip. Unfortunately, such dynamics is not realizable on a computer, but it is of interest for comparison with realizable algorithms.

We can say more about the dynamics given by Eq. (62). Let us first convert the equation into continuous time which is more convenient for analytic treatment, and which is a valid description for moderately large system. Introducing $`t=i/N`$ and $`\mathrm{\Delta }t=1/N`$ and define $`h(E,i/N)=h_i(E)`$, we have

$$h(E,\frac{i+1}{N})h(E,\frac{i}{N})=\underset{E^{}}{}h(E^{},i/N)\left(T(E^{}E)\delta _{E^{},E}\right).$$
(63)
Dividing both side by $`1/N`$, taking the limit of large $`N`$, we have

$$\frac{h(E,t)}{t}=\underset{E^{}}{}h(E^{},t)\stackrel{~}{T}(E^{}E),$$
(64)
where the continuous time transition matrix $`\stackrel{~}{T}`$ is related to the discrete step matrix by $`\stackrel{~}{T}(EE^{})=\left(T(EE^{})\delta _{E,E^{}}\right)N`$.

Two results were initially reported in ref. . Detailed derivations will be given in Appendices. Firstly, an explicit form of $`\stackrel{~}{T}`$ can be given for the one-dimensional Ising model with Glauber flip rate, as

$`\stackrel{~}{T}(kk1)`$ $`=`$ $`{\displaystyle \frac{k(2k1)}{N1}}(1+\gamma ),`$ (65)
$`\stackrel{~}{T}(kk+1)`$ $`=`$ $`{\displaystyle \frac{(N2k)(N2k1)}{2(N1)}}(1\gamma ),`$ (66)
where $`k=(E/J+N)/4`$, $`\gamma =\mathrm{tanh}\left(2J/(k_BT)\right)`$, and $`N`$ is the chain length. The diagonal term is computed from the relation

$$\stackrel{~}{T}(kk1)+\stackrel{~}{T}(kk)+\stackrel{~}{T}(kk+1)=0,$$
(67)
and the rest of the elements $`\stackrel{~}{T}(kk^{})=0`$ if $`|kk^{}|>1`$. Secondly, equation (64) is continuous in time but still discrete in energy. We can go one step further to consider the limit of both time and energy to be continuous. For transition matrix associated with canonical ensemble, we found a partial differential equation in such limit as

$$\frac{h(x^{},t^{})}{t^{}}=\frac{}{x^{}}\left(\frac{h(x^{},t^{})}{x^{}}+x^{}h(x^{},t^{})\right),$$
(68)
where $`x^{}`$ and $`t^{}`$ are properly scaled energy fluctuation and time:

$$x^{}=\frac{Eu_0N}{(N\overline{c})^{1/2}},u_0N=\overline{E},$$
(69)
and $`t^{}=At`$ with

$$A=\underset{N\mathrm{}}{lim}\underset{E}{}\stackrel{~}{T}(E\overline{E}),$$
(70)
where $`u_0`$ is the average energy per spin and $`\overline{c}=k_BT^2c`$ is the reduced specific heat per spin. This equation, (68), is equivalent to the one-dimensional quantum harmonic oscillator equation, thus the analytic solutions are readily obtained.

The most important consequence of this equation is that the relaxation time is proportional to the specific heat of the system. This result can also be seen from a less rigorous point of view. Since this artificial dynamics involves a random walk on the probability distribution of the energy, the characteristic relaxation time will be proportional to the square of the energy fluctuation, $`\left(EE\right)^2`$, which is in turn proportional to the specific heat.

## 5 Connections with Other Methods

### 5.1 Single histogram method

In the single histogram method , one performs a canonical ensemble simulation at a fixed temperature $`T_0`$, and collects the histogram $`H(E)`$. The histogram is proportional to $`n(E)\mathrm{exp}(E/k_BT_0)`$, so an estimate to the density of states is obtained from

$$n(E)H(E)\mathrm{exp}\left(\frac{E}{k_BT_0}\right).$$
(71)
Once we have the density of states, we can use it to evaluate thermodynamic quantities at any temperature.

Unfortunately, since $`H(E)`$ is approximately a Gaussian function with mean $`E_{T_0}`$ and variance $`k_BT_0^2cN`$, where $`c`$ is the specific heat per spin, the accuracy of the estimates deteriorates exponentially outside the energy window of order $`\sqrt{N}`$, c.f. Fig. 6. Detailed error analyses for energy and specific heat are given in . The region of good accuracy coincides with the critical region at a second order phase transition, so the single histogram method is still an extremely valuable tool to study phase transitions.

The transition matrix approach can also be used in a way similar to single histogram method, i.e., collecting the statistics of the transition matrix in a canonical simulation. Numerical comparison suggests that the two methods are comparable in accuracy. In fact, as we can see in Fig. 6, for a certain interval, the transition matrix gives results which are up to 10 times better, but become comparable or worse outside the limited range of $`E`$. If we use the two results to compute the average energy or heat capacity, we found that the errors are about the same . The reason is that the contributions to errors are dominated by the tails of the histogram distribution $`H(E)`$, at these ranges, the density of states is of comparable accuracy.

It is somewhat disappointing that the single histogram method and transition matrix Monte Carlo analysis are of the same accuracy. Some improvement can be made by a careful analysis using Baysian method . But it is unlikely that we can bring about an improvement of order $`\sqrt{N}`$ for the accuracy.

### 5.2 Multiple histogram and multi-canonical methods

Both of the multiple histogram method  and the multi-canonical method  give density of states over a wide range. While the multiple histogram method uses a collection of standard canonical simulations, the multi-canonical method uses only one simulation. In reality, a multi-canonical simulation needs to be iterated few times to converge to the desired ensemble. In this respect, the flat histogram method takes at most three iterations. The first iteration already gives excellent results, although there are noticeable biases for large systems; the second iteration with fixed flip rates greatly improves the accuracy; the third iteration would give correct sample average for the transition matrix as well as correct multi-canonical ensemble.

The additional benefit of using the transition matrix is improved accuracy comparing to other methods, within the same simulation runs. Figure 7 shows the accuracy of the density of states for two-dimensional Ising model on a $`16\times 16`$ lattice. We note that the accuracy is sensitive to the constraints imposed with the optimization. This extra accuracy comes about due to the nature of the samples that are taken. In histogram or multi-canonical methods, each new state gives one count in the histogram, while $`N`$ counts are collected from each state for the matrix. Naively, we expect an improvement by a factor of $`N`$ in terms of the variance, since each state contributes 1 for the histogram, and each state contributes a number of $`O(N)`$ for the transition matrix. While the accuracy of the transition matrix elements does improve as the system size increases (may be by $`1/\sqrt{N}`$ for the error), as has been pointed out by others , this accuracy is lost in the density of states. This is due to accumulation of errors, as the transition matrix elements only determine the ratio of the density of states, c.f. Eq. (55). If we use a simple iteration method starting from the ground state, we see that this extra accuracy in the matrix elements gets canceled exactly by the accumulation of errors. However, the optimization methods of determining $`n(E)`$ make the error analysis difficult. Here are some quantitative comparisons. We take the exact multi-canonical rate (algorithm $`1`$) in simulation, and collect both the transition matrix and the histogram. With transition matrix, the average relative errors of the density of states defined by Eq. (61) for the two-dimensional Ising model of size $`L=4`$, 8, 16, 32, and 50 with $`10^6`$ Monte Carlo steps in each run are 0.0003, 0.0012, 0.0037, 0.011, 0.024, respectively. The corresponding results computed by histograms are 0.0033, 0.010, 0.027, 0.058, 0.10, respectively. In general, transition matrix method performs better than histogram based methods, as we have shown numerically.

### 5.3 Transition probability methods

A proposal to use the transition matrix was given by Smith and Bruce  in 1995 in connection with multi-canonical simulation. This is further developed by Fitzgerald et al . The canonical transition probability (CTP) method  also estimates the transition matrix and uses energy detailed balance equations to estimate the canonical distribution. In the simplest version, instead of collecting the histogram $`H(E)`$, a matrix $`H(EE^{})`$ is incremented by 1 for every Monte Carlo move from state with energy $`E`$ to state with energy $`E^{}`$. Clearly, this quantity is an estimator of $`h(E)T(EE^{})`$. The transition matrix is obtained by $`H(EE^{})/H(E)`$. Both of the above methods and the present method are similar in the way that transition matrix is used. However, there are two important differences: (1) in CTP method, only the current move is used for statistics, not all possible moves of the state; (2) the simulation is performed in canonical ensemble at a given temperature.

### 5.4 Wang-Landau method

Very recently, Wang and Landau  proposed an intriguing method to determine the density of states. The dynamics follows the usual multi-canonical simulation or entropic sampling, by the single-spin-flip rate $`\mathrm{min}(1,n(E)/n(E^{}))`$. However, $`n(E)`$ is not a constant, but is updated at each step of trial move with

$$n(E)n(E)f$$
(72)
for the current energy $`E`$. This is somewhat like the Lee’s method of entropic sampling, but to some extent the updating of the weights are done at every move. If $`f`$ were greater than 1, the algorithm would never converge, so the idea is to reduce $`f`$ after some Monte Carlo steps, by $`ff^{1/2}`$, for example. A criterion of flatness of the histogram was used to determine if $`f`$ should be reduced.

Wang and Landau’s idea can be adopted in the context of transition matrix. For example, we can consider updating the logarithm of density of states, $`S(E)=\mathrm{ln}n(E)`$, using the information from the transition matrix by

$$S(E)S(E)+\eta \left(S^{\mathrm{pred}}(E)S(E)\right),$$
(73)
where $`0<\eta <1`$ is some small parameter and

$$S^{\mathrm{pred}}(E)=\frac{1}{M}\underset{E^{}}{}\left(S(E^{})+\mathrm{ln}\frac{T_{\mathrm{}}(E^{}E)}{T_{\mathrm{}}(EE^{})}\right)$$
(74)
is the predicted logarithmic density of states, based on $`M`$ possible hops from $`E^{}`$ to the current $`E`$. If we already know the ground state degeneracy, we can fix it to the constant. Unlike the updating rule $`n(E)n(E)f`$ which makes $`n(E)`$ grow indefinitely, $`S(E)`$ will converge to the exact value. This generalization does give more accurate results than algorithm 0 if it converges to flat histogram. However, it appears to have the problem of sticking to a Gaussian-like distribution for the histogram for large systems.

We made an extensive test of the accuracy of the random walk algorithm of Wang and Landau. In general, the random walk algorithm is far inferior to the flat histogram algorithm in terms of accuracy and rate of convergence to flatness, particularly for small systems. For the original implementation of the method, once the system passes the transient period, the error becomes independent of the total Monte Carlo steps used and the system sizes, and is primarily determined by how slowly $`f`$ is reduced. This feature is useful for its robustness, particularly for large systems.

In Fig. 8 we plot average errors in the density of states as defined by Eq. (61) for fixed Monte Carlo steps of $`10^6`$ ($`2\times 10^6`$ for the two-stage algorithm). It is clear that all methods have bigger errors for larger systems. But the random walk algorithm is generally order of magnitudes worse than the best transition matrix based methods. The sharp increase of the error with system sizes from 16 to 32, and to 50 for the random walk algorithm is an indication that with $`10^6`$ steps, the system is still in transients for these sizes. If $`10^7`$ Monte Carlo steps are used, the random walk algorithm comes closer to the accuracy of algorithm 0 at size $`L=50`$. The two-stage method (algorithm 0 followed by algorithm $`1`$, both using $`N`$-fold way) gives the best performance.

## 6 Generalization and Some Comments

The transition matrix approach can be used for continuous degrees of freedom. In such applications, both the dynamic variables and energy spectrum have to be discretized. Muñoz et al have applied the broad histogram method to the XY model . The important issue here is how to control the numerical error caused by discretization.

It is straightforward to generalize the transition matrix to more than one macroscopic variable, such as energy and magnetization; in fact, this has already been done  with the broad histogram method. This approach may have problems, particularly when the Hamiltonian is complicated. First, the matrix may be too large to handle in general. Second, with more elements to fill, the statistics for individual entry is poor. This makes the method less robust.

The transition matrix simulation can be parallelized very efficiently, where each processor works on separate configurations, using and updating a common transition matrix. Each processor can be limited to work on a range of energies. The advantage of this is clearly a fast way of approaching the flat distribution for the histogram.

## 7 Conclusion

Starting from the detailed balance equation, we have formulated the transition matrix in energy. The infinite-temperature version of this matrix serves as the basic data from which we determine the density of states and at the same time is used for construction of flat-histogram and equal-hit algorithms. This method of simulation together with optimization method to determine the density of states offers a better way of computing thermodynamic quantities by Monte Carlo simulation. In such an approach, a single simulation produces the whole function of temperature (or other parameter of the model) through re-weighting. It is efficient and easy to implement. As the use of accumulated average for the transition rates causes slow convergence for large systems, a two-stage iteration is recommended and is enough to get the best convergence. Dynamically, the flat-histogram algorithm for long simulations is equivalent to the multi-canonical method. Using the equal-hit algorithm together with $`N`$-fold way offers additional benefits.

## Acknowledgements

We thank Fugao Wang for providing us his programs for the random walk simulation. We also thank Lik Wee Lee, Bernd A. Berg, and Fei Fan for comments on the manuscript. J.-S. W. is supported by NUS and Singapore-MIT Alliance research grants.

## Appendices

## Appendix A Exact expression for the energy transition matrix in one dimension

We consider the single-spin-flip dynamics with a random pick of spins and the Glauber flip rate

$$a(EE+\mathrm{\Delta }E)=\frac{1}{2}\left[1\mathrm{tanh}\frac{\mathrm{\Delta }E}{2k_BT}\right].$$
(A.1)
In one dimension for the Ising model, the quantity $`N(\sigma ,\mathrm{\Delta }E)_E`$ can be evaluated exactly. Let us first define a set of new variables $`n_i=(\sigma _i\sigma _{i+1}+1)/2`$. $`n_i`$ is 1 for a satisfied bond, and 0 for an unsatisfied bond. The mapping from $`\sigma `$ to $`n`$ is unique modulo an overall spin up-down symmetry. We assume periodic boundary condition and lattice size $`N`$ to be even. There are three possible energy changes, $`4J`$, 0, $`+4J`$. If the original spin of a site and the spins of two neighboring sites all have the same sign, it contributes one count for the total $`N(\sigma ,4J)`$. In terms of $`n_i`$, it requires two consecutive satisfied bonds, i.e., $`n_in_{i+1}=1`$. Thus we can express $`N(\sigma ,+4J)`$ in terms of $`n_i`$ as

$$N(\sigma ,+4J)=\underset{i=1}{\overset{N}{}}n_in_{i+1}.$$
(A.2)
Note that only $`N1`$ variables $`n_i`$, $`i=1,2,\mathrm{},N1`$, are independent (since $`_{i=1}^Nn_i`$ must be even).

A microcanonical average at energy $`E`$ needs to be carried out. Let us use $`k`$ to label the equally spaced energy levels, $`k=0,1,2,\mathrm{},N/2`$. Then $`E=NJ+4Jk`$, and $`_in_i=n=N2k`$. The microcanonical average can be expressed as a summation over all $`n_i`$ subject to $`k=(E/J+N)/4`$ being an integer constant. Thus we have

$$n(E)N(\sigma ,+4J)_E=2\underset{_in_i=n}{}\underset{i=1}{\overset{N}{}}n_in_{i+1},$$
(A.3)
and similarly

$$n(E)N(\sigma ,4J)_E=2\underset{_in_i=n}{}\underset{i=1}{\overset{N}{}}(1n_i)(1n_{i+1}).$$
(A.4)
The factor 2 is due to the two-to-one mapping from $`\sigma `$ to $`n`$. In order to compute the above sums, we consider the statistical mechanics problem of a one-dimensional lattice gas with the Hamiltonian,

$$H=ϵ\underset{i=1}{\overset{N}{}}n_in_{i+1}h\underset{i=1}{\overset{N}{}}n_i.$$
(A.5)
The partition function of this system (at $`1/k_BT=1`$) is

$$Z=\underset{\{n_i\}}{}\mathrm{exp}\left(ϵ\underset{i=1}{\overset{N}{}}n_in_{i+1}+h\underset{i=1}{\overset{N}{}}n_i\right).$$
(A.6)
Taking the derivative with respect to $`ϵ`$, we have

$$\frac{Z}{ϵ}|_{ϵ=0}=\underset{n=0}{\overset{N}{}}\left(\underset{{\scriptscriptstyle n_i}=n}{}\underset{i=1}{\overset{N}{}}n_in_{i+1}\right)e^{hn}=\underset{n=0}{\overset{N}{}}\mathrm{\Delta }_{\frac{Nn}{2}}^+\mu ^n,$$
(A.7)
where $`\mu =e^h`$. The desired quantity is obtained from the generating function $`Z`$ as $`N(\sigma ,+4J)_E=2\mathrm{\Delta }_{(E/J+N)/4}^+/n(E)`$. The partition function is obtained by the standard trick of transfer matrix. We find

$$Z=\lambda _+^N+\lambda _{}^N,$$
(A.8)
where

$$\lambda _\pm =\frac{1}{2}\left(1+\mu e^ϵ\pm \sqrt{1+4\mu 2\mu e^ϵ+\mu ^2e^{2ϵ}}\right)$$
(A.9)
are the eigenvalues of the matrix

$$\left[\begin{array}{cc}1& \sqrt{\mu }\\ \sqrt{\mu }& e^ϵ\mu \end{array}\right].$$
(A.10)
After some algebra, we find

$$\frac{Z}{ϵ}|_{ϵ=0}=N\mu ^2(1+\mu )^{N2}=N\underset{n=0}{\overset{N2}{}}\frac{(N2)!}{n!(N2n)!}\mu ^{n+2}.$$
(A.11)
Thus

$$n(E)N(\sigma ,+4J)_{(N+4k)J}=2\mathrm{\Delta }_k^+=\frac{2N(N2)!}{(2k)!(N2k2)!}.$$
(A.12)
A similar derivation from a slightly different Hamiltonian gives

$$n(E)N(\sigma ,4J)_{(N+4k)J}=2\mathrm{\Delta }_k^{}=\frac{2N(N2)!}{(2k2)!(N2k)!}.$$
(A.13)
Combining the above results with the density of states for the Ising model, $`n(E)=2N!/[(2k)!(N2k)!]`$, which is readily obtained by the combinatorial problem of putting $`2k`$ unsatisfied bonds in $`N`$ places, we obtain the expressions given in Eq. (65) and (66).

## Appendix B Transition matrix dynamics in the continuum limit

We start from the dynamical equation

$$\frac{h(E,t)}{t}=\underset{E^{}}{}h(E^{},t)\stackrel{~}{T}(E^{}E),$$
(A.14)
where time $`t`$ is continuous but energy $`E`$ is discrete. The aim is to consider the continuous energy limit. This limit is natural and is a very good approximation for large systems. We follow the general method known as $`\mathrm{\Omega }`$ expansion . Let us introduce a new variable,

$$x=\frac{ENu_0}{\sqrt{N}},$$
(A.15)
where $`u_0`$ will be determined later. Since $`E`$ is of order $`N`$, naively $`x`$ is of order $`\sqrt{N}`$. However, by choosing $`u_0`$ to be the average of $`E/N`$, we cancel the leading $`N`$ dependence; $`x`$ is actually of order 1, measuring the relative fluctuation around mean. We look for nontrivial solution in variable $`x`$ in the scaling limit of $`N\mathrm{}`$, keeping $`x`$ fixed. More precisely, we find equation in $`x`$ such that the coefficients of the differential equation are independent of $`N`$. Consider the function in terms of $`x`$ as $`\stackrel{~}{h}(x)=h(Nu_0+x\sqrt{N})`$. We also write $`\stackrel{~}{T}`$ in $`x`$ as

$$\stackrel{~}{T}_i(x)=\stackrel{~}{T}(Nu_0+x\sqrt{N}+iaNu_0+x\sqrt{N}),$$
(A.16)
where $`i=0,\pm 1,\pm 2,\mathrm{},\pm d`$ is the change of energy associated with $`E^{}`$, and $`E^{}=E+ia`$, $`a=4J`$. We assume a $`d`$-dimensional Ising model in the derivation. Replacing $`E`$ by $`Nu_0+x\sqrt{N}`$, $`E^{}`$ by $`Nu_0+x\sqrt{N}+ia`$, Equation (A.14) becomes

$$\frac{\stackrel{~}{h}(x,t)}{t}=\underset{i=0,\pm 1,\mathrm{},\pm d}{}\stackrel{~}{T}_i(x)\stackrel{~}{h}(x+\frac{ia}{\sqrt{N}},t).$$
(A.17)
The crucial step now is to take Taylor expansion assuming $`a/\sqrt{N}`$ small, and to find equation that is leading order in the large $`N`$ limit. For notational convenient, we shall drop the tildes on $`\stackrel{~}{T}`$ and $`\stackrel{~}{h}`$, which actually denote different functions.

We know in the limit $`N\mathrm{}`$, $`h(x,t)`$ and $`T_i(x)`$ are smooth functions in $`x`$. Then (omitting the variable $`t`$ for clarity)

$$h\left(x+\frac{ia}{\sqrt{N}}\right)=h(x)+\frac{ia}{\sqrt{N}}h^{}(x)+\frac{1}{2}\left(\frac{ia}{\sqrt{N}}\right)^2h^{\prime \prime }(x)+O(1/N^{3/2})h^{\prime \prime \prime }(x).$$
(A.18)
The primes denote derivatives with respect to $`x`$. We should note that in the large $`N`$ limit with $`x`$ fixed, $`h(x)`$ and its derivatives do not contain $`N`$. The $`N`$ dependence is made explicit by the above expansion. Substituting Eq. (A.18) into Eq. (A.17), we have

$$\frac{h(x)}{t}=A(x)h(x)+B(x)h^{}(x)+C(x)h^{\prime \prime }(x)+O(1/\sqrt{N})h^{\prime \prime \prime }(x),$$
(A.19)
where

$`A(x)`$ $`=`$ $`{\displaystyle \underset{i}{}}T_i(x),`$ (A.20)
$`B(x)`$ $`=`$ $`{\displaystyle \underset{i}{}}T_i(x){\displaystyle \frac{ia}{\sqrt{N}}},`$ (A.21)
$`C(x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}T_i(x){\displaystyle \frac{(ia)^2}{N}}.`$ (A.22)
Naively, since $`T`$ is $`O(N)`$, we expect $`A(x)O(N)`$, $`B(x)O(\sqrt{N})`$, and $`C(x)O(1)`$ and the equation does not have a well-defined large $`N`$ limit. But this is not true, due to special relations among $`T_i`$. Two relations are important in the derivation below to show that both $`A(x)`$, $`B(x)`$, and $`C(x)`$ are of $`O(1)`$ and the third derivative can be dropped in the large $`N`$ limit.

The existence of an equilibrium implies

$$\underset{E}{}\stackrel{~}{T}(E^{}E)=0.$$
(A.23)
Expressed in $`T_i(x)`$, it is

$$\underset{i=0,\pm 1,\mathrm{},\pm d}{}T_i\left(x\frac{ia}{\sqrt{N}}\right)=0.$$
(A.24)
Thus, using $`T_0(x)`$ from the above

$`A(x)`$ $`=`$ $`T_0(x)+{\displaystyle \underset{i0}{}}T_i(x)`$ (A.25)
$`=`$ $`{\displaystyle \underset{i0}{}}\left[T_i(x)T_i\left(x{\displaystyle \frac{ia}{\sqrt{N}}}\right)\right]`$ (A.26)
$`=`$ $`{\displaystyle \underset{i0}{}}\left[T_i^{}(x){\displaystyle \frac{ia}{\sqrt{N}}}{\displaystyle \frac{1}{2}}T_i^{\prime \prime }(x){\displaystyle \frac{(ia)^2}{N}}+\mathrm{}\right].`$ (A.27)
The last equation above used a Taylor expansion. Since $`T_i(x)`$ has a scaling form $`T_i(x)Ng(x/\sqrt{N})`$ in the large $`N`$ limit, we find that the $`k`$-th derivative of $`T_i(x)`$ at $`x=0`$ is of order $`T_i^k(0)O(N^{1k/2})`$. Thus $`T_i^{}(x)`$ is of order $`\sqrt{N}`$, and we can safely replace $`x`$ by 0, and find

$$A(x)=\left\{\underset{i}{}T_i^{}(0)\frac{ia}{\sqrt{N}}\right\}+O\left(\frac{1}{\sqrt{N}}\right)=A+O\left(\frac{1}{\sqrt{N}}\right).$$
(A.28)
For $`B(x)`$, we make an expansion in $`x`$ for $`T_i(x)`$ and find

$`B(x)`$ $`=`$ $`{\displaystyle \underset{i}{}}T_i(x){\displaystyle \frac{ia}{\sqrt{N}}}`$ (A.29)
$`=`$ $`{\displaystyle \underset{i}{}}\left\{T_i(0){\displaystyle \frac{ia}{\sqrt{N}}}+T_i^{}(0){\displaystyle \frac{ia}{\sqrt{N}}}x+{\displaystyle \frac{1}{2}}T_i^{\prime \prime }(0){\displaystyle \frac{ia}{\sqrt{N}}}x^2+\mathrm{}\right\},`$ (A.30)
where in the last formula, the first term is of order $`\sqrt{N}`$, the second term is of order 1, and last term is of order $`1/\sqrt{N}`$. If the first term were there, we would have an ill-defined limit. So we must require that

$$D=\underset{i}{}T_i(0)\frac{ia}{\sqrt{N}}=0.$$
(A.31)
This is in fact a condition to fix $`u_0`$. We shall show that this condition requires $`u_0=E/N_T`$, the canonical average of energy per spin.

We evoke the energy detailed balance equation, (15). In terms of the new variables $`x`$ and $`i`$, it is

$$T_i(x)h_{eq}\left(x+\frac{ia}{\sqrt{N}}\right)=T_i\left(x+\frac{ia}{\sqrt{N}}\right)h_{eq}(x).$$
(A.32)
Let $`\delta =ia/\sqrt{N}`$, Taylor expanding the terms involving $`\delta `$, we find:

$$T_i(0)T_i(0)=T_i(0)\left[\mathrm{ln}h_{eq}(0)\right]^{}\delta +T_i^{}(0)\delta +O\left(\frac{1}{N}\right).$$
(A.33)
Note that the first term is of order $`O(\sqrt{N})`$, the second term of order $`O(1)`$. It is important to realize that we are looking for the scaling limit of $`N\mathrm{}`$, fixing $`x`$.

Substituting this equation into the expression for $`D`$, we find

$$D=\underset{k=1,2,\mathrm{},d}{}\left[T_k(0)T_k(0)\right]\frac{ka}{\sqrt{N}}\underset{k=1,2,\mathrm{},d}{}k^2T_k(0)[\mathrm{ln}h_{eq}(0)]^{}\frac{a^2}{N}.$$
(A.34)
The requirement that $`D=0`$ is equivalent to say that $`x=0`$ is at the extreme of equilibrium probability distribution.

When the first term in $`B(x)`$ is set to 0, we have

$$B(x)=\left\{\underset{i}{}T_i^{}(0)\frac{ia}{\sqrt{N}}\right\}x+O\left(\frac{1}{\sqrt{N}}\right)=Ax+O\left(\frac{1}{\sqrt{N}}\right).$$
(A.35)
The coefficient $`C(x)`$ is a constant to leading order in $`N`$:

$`C(x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}T_i(x){\displaystyle \frac{(ia)^2}{N}}`$ (A.36)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}\left\{T_i(0){\displaystyle \frac{(ia)^2}{N}}+T^{}(0){\displaystyle \frac{(ia)^2}{N}}x+\mathrm{}\right\}`$ (A.37)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}T_i(0){\displaystyle \frac{(ia)^2}{N}}+O({\displaystyle \frac{1}{\sqrt{N}}})=C+O({\displaystyle \frac{1}{\sqrt{N}}}).`$ (A.38)
The equilibrium distribution of energy for large system is a Gaussian distribution with mean $`E_T=Nu_0`$, and variance $`N\overline{c}`$ where $`\overline{c}`$ is reduced specific heat per spin, thus

$$h_{eq}(x)=\frac{1}{\sqrt{2\pi \overline{c}}}\mathrm{exp}\left(\frac{x^2}{2\overline{c}}\right).$$
(A.39)
Substituting this result into the partial differential equation in equilibrium

$$\frac{h}{t}=Ah+Axh^{}+Ch^{\prime \prime }=0,$$
(A.40)
we find $`C=A\overline{c}`$. This same relation can also be obtained using detailed balance equation. Changing variables from $`t`$ to $`t^{}=At`$ and from $`x`$ to $`x^{}=x/\sqrt{\overline{c}}`$, we obtain Eq. (68).

## Appendix C Data file for errors

The file below is the raw data for various errors. Formally this is not part of the paper. We included here which could be useful for future benchmarking use.

```
Convergence Test Data

MCDIS = -1, SIZE L = 5, N-Fold-Way, 2D Ising, algo -1 to 12
e_T = sum | T - T_ex |

mcs  -1     0      1      2      3      4      5      6      7      8      9      10    11     12

1    14.961 15.059 15.043 13.995 14.299 15.357 13.330 14.241 14.455 15.883 13.859 16.06 13.776 15.237
10   6.480  8.107  7.795  6.509  6.33   8.83   5.52   5.40   6.79   9.91   6.32   9.68  6.41   8.09
100  0.681  0.819  0.736  0.759  2.28   1.09   0.710  0.713  0.667  1.29   1.47   1.00  1.171  0.825
1000 0.188  0.204  0.210  0.234  1.50   0.320  0.221  0.221  0.199  0.226  0.581  0.227 0.370  0.256
1e4  0.0608 0.0644 0.0663 0.080  1.00   0.14   0.070  0.070  0.064  0.076  0.539  0.079 0.256  0.094
1e5  0.0181 0.0206 0.0215 0.030  0.698  0.070  0.0219 0.0224 0.0201 0.0246 0.5353 0.026 0.2136 0.0349
1e6  0.0060 0.0064 0.0070 0.0098 0.43   0.044  0.0067 0.0077 0.0064 0.0085 0.533  0.008 0.131  0.0144
1e7 0.00186 0.0020 0.0021 0.0040 0.357  0.0172 0.0024 0.0023 0.0021 0.0030 0.5342 0.0026 0.0926 0.0047
1e8         0.00062 0.00066      0.29   0.012                                            0.0842

same as above parameters, but with single-spin-flip (i.e., no N-fold)

mcs  -1     0      1      2      3      4      5      6      7      8      9      10    11     12
1   16.448  15.26  15.263 15.178 15.307 15.31  15.18  15.31  15.41  15.30  15.31  15.31 15.54  15.30
10   9.329  9.162  9.165
100  1.286  1.583  1.599
1000 0.287  0.312  0.310
1e4  0.093  0.097  0.095
1e5  0.0287 0.0311 0.0314

           algo no (error at N-fold mcs=1e4)
best        -1     (0.061)
very good   0,1,7  (0.065)
good        5,6,8  (0.07)
OK          2,10   (0.08)
not so good 4,12   (0.1)
bad (converge) 3,11 (>0.25)
don’t converge 9   (0.5)

All the rest has MCDIS=0 (no MCS discarded)

2D Ising model, e_n = (1/(N-1)) sum | n(E)/n_ex(E) - 1 |
                  algo 0 (N-fold-way)
mcs      L=4          L=8           L=16       L=32        L=50
1        1.200
10       0.152        1.38^
100      0.0340       0.212         1.3^
1000     0.0109       0.0492        0.20       1^
1e4      0.00342      0.0154        0.0568     0.28        0.31
1e5      0.00108      0.0054        0.0189     0.084       0.23
1e6      0.000338     0.00154       0.0066     0.07        0.22+/-0.02
1e7      0.000110     0.00051       0.0024     0.042       0.191+/-0.0138
1e8      0.000035     0.000190      0.0009     0.022       0.16

^wide distribution, occasionally not converging

tunnel   18.983      144.68         983         6448
time

2D Ising model, e_n = (1/(N-1)) sum | n(E)/n_ex(E) - 1 |
            algo 1 (N-fold-way, equal-hit)
mcs      L=4          L=8           L=16       L=32        L=50
1        1.20
10       0.142        1.27
100      0.0410       0.178
1000     0.0134       0.051         0.26
1e4      0.0042       0.0152        0.048      0.18^       0.4
1e5      0.00138      0.0050        0.015      0.04        0.12^
1e6      0.00044
1e7      0.00015

2D Ising model, e_n = (1/(N-1)) sum | n(E)/n_ex(E) - 1 |
         algo -1 (N-fold-way) (exact rate, best we can do)
mcs      L=4          L=8           L=16       L=32        L=50
1       1.32
10      0.131       14
100     0.031        0.18
1000    0.0097       0.038         0.16
1e4     0.0030       0.0125        0.041       0.13        0.23
1e5     0.00096      0.00383       0.0121      0.037       0.07
1e6     0.00031      0.00119       0.0037      0.0114      0.024
1e7     0.000095     0.00040       0.00122     0.00337     0.007
1e8     0.000032     0.000103      0.00045     0.00091

Ising model, same as above algo -1 N-fold-way, but use histogram to
compute n(E), i.e. Berg/Lee method
mcs      L=4          L=8           L=16       L=32        L=50
1e5     0.010         0.033        0.0733      0.32        0.5
1e6     0.0033        0.0099       0.027       0.058       0.10

2D Ising model, e_n = (1/(N-1)) sum | ...|
two stage simulation, stage 1, using algo 0,
                      stage 2, using algo -1 of n(E) generated from stage 0.
same run length for the two stages
mcs       L=4         L=8       L=16           L=32        L=50
100      0.0311       0.26
1000     0.0096       0.0404    0.10
10000    0.0030       0.0122    0.0401         0.2         0.5
1e5      0.00100      0.00349   0.0131         0.038       0.075
1e6      0.00036      0.00123   0.0034         0.0133      0.017
1e7                   0.00038   0.0016         0.0029      0.0074
1e8                             0.0004         0.0010

Modified Wang FG/Wang JS method (algo 13’), e_n = (1/(N-1) sum | ...|
using S(E) <- S(E) + eta (S^pred - S(E))
(eta = 0.1)  (two numbers - using Tmatrix/using S(E) directly)
mcs  L=4            L=8            L=16          L=32
100
1000 0.00936/0.358
1e4  0.00306/0.0087 0.39/0.96
1e5  0.00096/0.0028 0.0039/0.015  0.99/0.99
1e6                 0.00115/0.004 0.71/0.99
1e7                 0.00049/0.0014

Wang Fugao’s program (fixed MCS, f_max = 2.718, f_min=1, 80%H, check every 1k steps)
e_n = < | n(E)/n_ex(E) - 1 | >
mcs       L=4         L=8       L=16           L=32        L=50
1e4      0.172      10^
1e5      0.065       0.073     1
1e6      0.066       0.049     0.066           0.36        9(+/-4)
1e7      0.055       0.073     0.043           0.05        0.18
1e8      0.069       0.047     0.07            0.025       0.14

Wang FG method e_n (my program algo 13, 80%H, checks 30 times, single-spin-flip)
mcs       L=4         L=8       L=16           L=32        L=50
1e4       0.079       1.62
1e5       0.034       0.071     0.7
1e6       0.012       0.031     0.056         0.42         5.5
1e7       0.0044      0.010     0.025         0.06         0.14
1e8       0.0013                                           0.02

```


format: text

text:
# Chromaticity effects in microlensing by wormholes

## 1 Introduction

A wormhole is a region of space-time with non-trivial topology. It has two mouths connected by a throat. The mouths are not hidden by event horizons, as in the case of black holes, and, in addition, there is no singularity that could avoid the passage of particles from one side to the other. After the pioneering paper by Morris & Thorne , wormhole solutions to the Einstein field equations has been extensively studied in the literature (see Refs.  and references cited therein).

Static wormhole structures require the violation of the average null energy condition in the wormhole throat in order to exist. Plainly stated, this means that the matter threading the wormhole must exert gravitational repulsion to stay stable against collapse. Although there are known violations to the energy conditions (e.g. the Casimir effect), it is currently far from clear whether large macroscopic amounts of “exotic matter” exist in the nature. If natural wormholes actually exist in the universe (e.g. if the original topology after the Big-Bang was multiply connected), then there should be observable electromagnetic signatures that could lead to their identification.

Cramer et al.  and Torres et al.  have recently studied the microlensing effects of a negative mass (e.g. a natural stellar-size wormhole) on the light from a background point-like source \[for a macrolensing study see \]. These authors have shown that the typical lightcurves expected from microlensing events produced by wormholes should be very different from the lightcurves in ordinary (positive mass) microlensing of point sources. In the standard case, a time symmetric burst in the flux density from the background source occurs. If the lens, instead, is a wormhole, the gravitational repulsion creates an obscure umbra region, deflecting light rays that should, otherwise, reach the observer. The wormhole basically acts as a divergent lens. However, light is concentrated on the border of the umbra region producing two separated intensity enhancements that are observed before and after the occultation event. These two burst are individually asymmetric under time reversal, although they are mirror images of each other. In the first burst, the observed flux increases towards a divergency and then drops to zero. The infinity, or “caustic”, and the vertical drop to zero are consequences of the ideal model adopted in the calculations, with a point source at the background. A more realistic model should consider an extended source with a given intensity distribution over it.

Recently, Anchordoqui et al.  searched in existent gamma-ray bursts databases for signatures of wormhole microlensing. Although they detected some interesting candidates, no conclusive results could be obtained. Peculiarly asymmetric gamma-ray bursts , although highly uncommon, might be probably explained by more conventional hypothesis, like precessing fireballs (see, for instance, Ref. ). Even in the case of galactic microlensing, highly asymmetric curves can be the effect of microlensing by binary stellar systems . Hence, although binary lenses produce asymmetric lightcurves with different shapes than those that would be produced by wormholes, discriminators other than the lightcurves are needed in order to identify wormhole microlensing events.

In this paper, we expand the analysis of Cramer et al.  and Torres et al.  to the extended source case. This allow us to present more realistic lightcurves for wormhole microlensing events. Then, using this formalism, we compute the effects of a finite source extent on the spectral features of microlensing. We show that limb darkening of the intensity distribution on a stellar source induces specific chromaticity effects that are very different from what is expected in the positive mass lens case. Thus, multi-color optical observations can be used to search for galactic natural wormholes. In addition, we study the extragalactic wormhole microlensing of background AGNs taking into account the variations of the source size with the observing wavelengths. This extends the previous work by Torres et al. , who first considered these kind of events at gamma-rays adopting a point-like source. We show that also in this case the spectral evolution during the wormhole microlensing event presents peculiar signatures that allow its identification.

## 2 Microlensing by wormholes in the extended source case

The total amplification $`A=I_{\mathrm{obs}}/I_0`$ –where $`I_{\mathrm{obs}}(I_0)`$ is the observed (intrinsic) luminosity– of a background point source due to gravitational lensing by a single point mass is given by

$$A_0=\frac{B_0^2\pm 2}{B_0\sqrt{B_0^2\pm 4}},$$
(1)
where the plus (minus) sign corresponds to the positive (negative) mass case, and

$$B_0=\frac{b_0}{R_\mathrm{E}}$$
(2)
is the shortest lens-source separation (on the lens plane) in units of the Einstein ring radius, $`R_\mathrm{E}`$, defined as

$$R_\mathrm{E}=\sqrt{\frac{4G|M|}{c^2}\frac{D_{\mathrm{ol}}D_{\mathrm{ls}}}{D_{\mathrm{os}}}}.$$
(3)
In this last expression, $`D_{\mathrm{os}}`$ is the observer-source distance, $`D_{\mathrm{ol}}`$ is the observer-lens distance, $`D_{\mathrm{ls}}`$ is the lens-source distance, and $`M`$ is the mass of the gravitational lens.111It should be kept in mind that any form of negative masses, not only wormholes, would produce the same observational effects.

When the source is extended, we must take into account the contributions coming from the different parts of it. The amplification observed at a frequency $`\nu `$, then, results:

$$A_\nu =\frac{_0^{2\pi }_0^r_{}_\nu (r,\phi )A_0(r,\phi )r𝑑r𝑑\phi }{_0^{2\pi }_0^r_{}_\nu (r,\phi )r𝑑r𝑑\phi },$$
(4)
where $`(r,\phi )`$ are polar coordinates in a reference frame centered in the source (e.g. a star), $`r_{}`$ is the radius of the source and $`_\nu (r,\phi )`$ is its surface intensity distribution. For a radially symmetric distribution we have:

$$A_\nu =\frac{_0^{2\pi }_0^r_{}_\nu (r)A_0(r,\phi )r𝑑r𝑑\phi }{2\pi _0^r_{}_\nu (r)r𝑑r}.$$
(5)
If the lens is moving with constant velocity $`v`$, the lens-source separation evolves in time as

$$b(t)=\sqrt{\left(b_0+r\mathrm{sin}\phi \right)^2+\left(vt+r\mathrm{cos}\phi \right)^2}.$$
(6)
In the published version of this paper, this is illustrated in the first figure. In units of the Einstein radius, we have:

$$B(T)=\frac{b(t)}{R_E}=\sqrt{\left(B_0+R\mathrm{sin}\phi \right)^2+\left(T+R\mathrm{cos}\phi \right)^2},$$
(7)
where $`T=vt/R_E`$. Replacing now $`B_0`$ by $`B(T)`$ in Eq. (1), we obtain the time dependent amplification.

In order to compute the expected lightcurves for the case of negative mass lenses, we take the negative sign in Eq. (1) and replace it in Eq. (5). In Figure 2 we show different lightcurves for different impact parameters, calculated assuming a circular source of uniform brightness and radius $`R_{}=r_{}/R_E=0.1`$. These curves should be compared with the ideal curves presented by Torres et al.  in their Figure 1. The main difference is that in the extended source case the divergencies are eliminated, and there is not a discontinuous transition to or from the umbra region. The basic features of the microlensing event, however, remain. For $`B_0<2`$, we have two successive intensity enhancements, separated by a period of absence of radiation. The individual bursts are quite asymmetric. In the case $`B_0>2`$, there is only one event, symmetric under time reversal and similar to what is observed in the standard case.

Real background sources in microlensing events can present non-uniform brightness distributions on their surfaces and a dependency of their emission with the observing frequency. These complications can result in chromaticity effects, i.e. in spectral changes induced by differential lensing during the event. In the next sections we shall compute such changes in order to establish whether they can provide a specific signature of wormhole microlensing.

## 3 Chromaticity effects in microlensing of stars

Stars are brighter in their center. The obscuration of the intensity profile of a star towards its border is known as “limb darkening”. They also radiate as a blackbody with an effective temperature $`T_{\mathrm{eff}}`$, and consequently their emission is frequency dependent. To compute the chromaticity effects due to these characteristics, we represent the intensity profile of typical star by :

$$_\nu (r)=1C_\nu \left(1\sqrt{1\left(\frac{r}{r_{}}\right)^2}\right),$$
(8)
where the limb-darkening coefficients in the I and U bands can be taken as $`C_{\nu _1}=0.503`$ and $`C_{\nu _2}=1.050`$, respectively, for a K-giant star with $`T_{\mathrm{eff}}=4750`$K.

Defining the dimensionless radius as $`R=r/R_E`$, we obtain

$$A_\nu =\frac{_0^{2\pi }_0^R_{}_\nu (R)A_0(R,\phi )R𝑑R𝑑\phi }{2\pi _0^R_{}_\nu (R)R𝑑R}$$
(9)
and

$$_\nu (R)=1C_\nu \left(1\sqrt{1\left(\frac{R}{R_{}}\right)^2}\right),$$
(10)
where $`R_{}=r_{}/R_E`$ is the dimensionless radius of the star. The colour change caused by limb darkening can then be computed by

$$\mathrm{\Delta }(m_{\nu _1}m_{\nu _2})=2.5\mathrm{log}\left(\frac{A_{\nu _1}}{A_{\nu _2}}\right).$$
(11)
In Figure 3 we show several color curves for both positive (left panel) and negative (right panel) mass lenses. The different curves correspond to different impact parameters $`k=B_0/R_{}`$ specified in the caption. The star radius was considered, as before, equal to $`R_{}=0.1`$. In the case of an ordinary lens, we see that as it gets closer to the star, the color of the observed source becomes redder due to the differential amplification of the coldest regions. When the lens transits towards the star interior, the hot center starts to dominate the amplification, producing a dramatic change in the slope of the color curve. This is basically the result recently found by Han et al. .

If we now look at the color curves for the negative mass lens case, we notice a clear contrast in the behaviour. The spectral changes start long before than in the standard situation. Initially, the source also becomes redder and then experiences a switch when shorter wavelengths begin to dominate. Contrary to what happens with positive masses, the spectral trend changes again, with the source appearing colder and colder until it vanishes in the umbra during the transit. When the source is seen again, the inverse behaviour is observed.

## 4 Chromaticity effects in extragalactic microlensing

Let us now turn to the extragalactic microlensing case, which was previously studied by Torres et al.  using point sources. We shall study the lightcurves produced by a stellar-mass extragalactic wormhole when it crosses the line of sight to a background compact source, namely and Active Galactic Nuclei (AGN). As discussed by Torres et al. , a critical requirement to observe a significant light magnification during the microlensing event is that the size of the source, when projected onto the lens plane, is not larger than the Einstein ring of the lensing mass. This is necessary because otherwise light from outside the ring could become dominant and smooth out the gravitationally induced variability. In the case of an AGN, since the size of the emitting regions varies with the wavelength, the gravitational amplification will be more effective at those wavelengths at which the source is more compact. The wavelengths where the amplifications are maximized are located at gamma-ray energies .

The compactness of the gamma-ray central region of an AGN is determined by the opacity to the propagation of gamma-ray photons in the X-ray field produced by the object. Since the opacity is a function of the photon energy, the size of the observable gamma-ray photosphere will also depend on the observing energy. Differential amplification during the microlensing event, then, should lead to chromaticity effects. In order to compute these effects, we shall adopt a radius of the photosphere given by the maximum height at which photons of energy $`E`$ are absorbed by pair creation in the X-ray radiation field of the inner accretion disk. According to Becker & Kafatos , this size is

$$r_\gamma E^{\alpha /(2\alpha +3)},$$
(12)
where $`\alpha `$ is the X-ray spectral index of the accretion disk radiation field. In our calculations, we shall adopt an average value $`\alpha =1.1`$ . The larger photospheres, then, are those observed at the higher energies. To compute the microlensing effects on the AGN, we define a reference source with radius $`r_{\mathrm{ref}}`$ and gamma-ray energy $`E_{\mathrm{ref}}`$, such that

$$R_\gamma (E)=R_{\mathrm{ref}}\left(\frac{E}{E_{\mathrm{ref}}}\right)^{\alpha /(2\alpha +3)},$$
(13)
where we have written $`r_\gamma `$ and $`r_{\mathrm{ref}}`$ in units of the Einstein radius $`R_\mathrm{E}`$.

The intensity of the source (without lensing) is uniform and its spectrum follows approximately a power law:

$$I_0(E)E^\xi $$
(14)
with $`\xi `$ typically in the range $`\xi (1.7,2.7)`$ . We can write this latter Eq. in the form:

$$I_0(E)=I_{\mathrm{ref}}\left(\frac{E}{E_{\mathrm{ref}}}\right)^\xi ,$$
(15)
where $`I_{\mathrm{ref}}`$ is the intensity of the reference source. The surface intensity distribution is assumed to be constant at a given energy

$$_0(E)=\frac{I_0(E)}{\pi \left(R_\gamma (E)\right)^2}.$$
(16)
Then, the total amplification of the source due to microlensing can be computed as

$$A=\frac{I}{I_0}=\frac{_0^{2\pi }_0^{R_\gamma (E)}_0(E)A_0(R,\phi )R𝑑R𝑑\phi }{_0^{2\pi }_0^{R_\gamma (E)}_0(E)R𝑑R𝑑\phi }.$$
(17)
Since $`_0(E)`$ does not depend on $`(R,\phi )`$, we have

$$A=\frac{_0^{2\pi }_0^{R_\gamma (E)}A_0(R,\phi )R𝑑R𝑑\phi }{\pi \left(R_\gamma (E)\right)^2}.$$
(18)
From Eq. (15) and Eq. (17) we obtain:

$$I=AI_0=AI_{\mathrm{ref}}\left(\frac{E}{E_{\mathrm{ref}}}\right)^\xi .$$
(19)
Let us now introduce the dimensionless minimum impact parameter $`k`$ as:

$$k=\frac{B_0}{R_\gamma (E)}.$$
(20)
Adopting $`\xi =2`$, $`R_{\mathrm{ref}}=0.2`$, and $`E_{\mathrm{ref}}=100`$ MeV, we can use Eq. (18) and Eq. (19) to compute the expected lightcurves and the spectral evolution during the microlensing event. The results are shown, for a variety of impact parameters, in Figures 4-7, where we define $`J`$ as the ratio between $`I`$ and $`I_{\mathrm{ref}}`$.

The main point to be emphasized in the spectral evolution of the lensed source is that, when the source is a wormhole, the spectrum emerges from the umbra at its lowest level (the bottom curve in the spectral plots) and then it jumps to its maximum, from where it decreases slowly keeping its slope in the log-log diagram. This occurs because after the transit the intensity have to recover from the umbra null values. The effects of the differential obscuration also appear as a clear break in the spectrum at early times since the transit. This break, with a hardening at low energies, results because the smaller and less energetic gamma-spheres emerge from the umbra after than the more energetic (and bigger) ones, and then are initially less amplified. The break in the spectrum is best appreciated in Figure 6, where the impact parameter is $`k=1`$.

In the case of positive lenses (upper set of panels in the figures) the evolution is smoother, decreasing the emission at all wavelengths after the transit. The chromaticity effects in this case are limited to a small slope change at medium energies.

## 5 Concluding remarks

In this paper we have developed the formalism for gravitational microlensing of extended sources when the lens is a wormhole-like object. We have studied the chromaticity effects introduced by the finite source size and the non-uniform intensity distribution on it, finding a set of very peculiar features that allow a clear identification of any lens having a negative energy density. These chromaticity effects, when combined with the observed lightcurves, provide a unique tool to study the possible existence of natural wormholes in the universe.

## Acknowledgments

This work has been supported by Universidad de Buenos Aires (UBACYT 01TW77, EE), CONICET (DFT, and PIP 0430/98, GER), ANPCT (PICT 98 No. 03-04881, GER), and Fundación Antorchas (GER and DFT).