Histogram methods for quantum systems: from reweighting to Wang-Landau sampling

Histogram methods for quantum systems: from reweighting to Wang-Landau sampling

Matthias Troyer^{I, II}; Fabien Alet^{I, II}; Stefan Wessel^I

^ITheoretische Physik, ETH Zürich, CH-8093 Zürich, Switzerland

^IIComputational Laboratory, ETH Zürich, CH-8092 Zürich, Switzerland

ABSTRACT

We discuss how reweighting and histogram methods for classical systems can be generalized to quantum models for discrete and continuous time world line simulations, and the stochastic series expansion (SSE) method. Our approach allows to apply all classical reweighting and histogram techniques for classical systems, as well as multicanonical or Wang-Landau sampling to the quantum case.

1 Introduction

Quantum Monte Carlo simulations have become almost as powerful as classical simulations, despite their usually more complex formulation and implementation.

Not only local update algorithms,[1] but also generalizations of the classical cluster Monte Carlo algorithms [2] have been developed for quantum systems, solving the problem of critical slowing down at some important second order phase transitions.[3,4,5] Some cluster update algorithms, such as the worm algorithm [6] have even originally been developed for quantum models before being applied to classical ones.[7]

Here we show that also classical reweighting, histogram sampling methods, and the Wang-Landau algorithm[8] can be adapted to quantum systems, expanding on technical issues omitted in Ref. [9]. Histogram reweighting [10] allows simulation results to be obtained in a range of temperatures (or couplings) from simulations performed at a single temperature. It is useful in accurately finding critical points, the Wang-Landau algorithm [8], which is related to multicanonical sampling [11] eases tunneling through free energy barriers at first order phase transitions.

We start with a short review of the classical histogram reweighting before discussing how it can be generalized to quantum systems, and finish with applications of a quantum version of Wang-Landau sampling.

2 Classical Monte Carlo simulations

2.1 Histogram reweighting

In a classsical system the thermal average of a quantity A at an inverse termperature b is defined as a sum over all configurations

where A_c is the measurement of the observable A in the configuration c and E_c the energy of that configuration. In a Monte Carlo simulations, a subset of M configurations {c_i}, drawn with probability p_c = exp(–bE_c) is sampled, and the thermal average estimated by the sample mean

This sampling scheme gives results only for the inverse temperature b, but actually there is much more information available than just the simple average Eq. (2). Indeed, in the search for a phase transition a range of temperatures needs to be explored and information at a nearby inverse temperature b¢ » b can be obtained from a simulation performed at b. This is done by reweighting the configurations sampled with the Boltzmann weight p_c = exp(–bE_c) to obtain averages for the Boltzmann weight = exp(–b¢E_c):

where Db = b¢ – b.

Instead of storing the full time series of measurements {} and energies {}, we now store a histogram H(E), counting how often the energy level E occurs in the time series {}, and the average A(E) of all the measurements performed on configurations with energy E:[10]

With these two histograms the average áA(b¢)ñ can be calculated as a sum over all energies

Histogram reweighting works well only if the configurations sampled at the inverse temperature b are also relevant at b¢, requiring that Db is small. Otherwise the errors become too large since there will not be sufficient entries in H(E) for the energies E important at b¢. Multiple histograms obtained at different temperatures can be used to broaden the accessible temperature range.[10]

2.2 Wang-Landau sampling

While cluster updates, optionally combined with histogram reweighting, are an efficient solution to the problem of critical slowing down at a second order phase transition they do not address the problem of tunneling between local free energy minima at a first order phase transition. There, the multicanonical algorithm[11] or Wang-Landau sampling can help.[8]

The key idea of the multicanonical method and Wang-Landau sampling is to calculate the density of states r(E) directly by a random walk in energy space instead of performing a canonical simulation at a fixed temperature. If, instead of the Boltzmann weight p(E) µ exp(–bE) one uses the inverse density of states p(E) µ 1/r(E), the probability of encountering a state with energy E changes from

to a constant probability

This algorithm should be much more efficient and solve the exponential tunneling problem, since one now performs a perfect random walk in energy space with a flat histogram. Naively one expects the scaling to be the O(N²) scaling of a random walk over N energy levels. In a detailed scaling analysis[12] we found that due to memory effects this is not a Markovian random walk. The scaling is slightly worse, but remains polynomial, for homogeneous non-frustrated and frustrated systems, and becomes exponential only for spin glasses.

Besides solving the tunneling problem at first order phase transitions, the Wang-Landau algorithm calculates the density of states r(E) in an iterative procedure. Knowledge of r(E) then allows the direct calculation of the free energy from the partition function

The internal energy, entropy, specific heat and other thermal properties are easily obtained as well, by differentiating the free energy.

3 Quantum Monte Carlo

3.1 Generalizing histogram reweighting and Wang-Landau sampling

Since simulations of quantum systems suffer from the same problems as classical simulations, in particular from long tunneling times at first order phase transitions and the inability to calculate the free energy directly, an extension of this algorithm to quantum systems is highly desired. The generalization is not immediately obvious since the partition function of a quantum system

cannot be cast in the form of Eq. (9) unless the complete spectrum of the Hamilton operator H is known.

Instead of a representation like Eq. (9) we will aim for a generalized representation of the form

where the tuple (n₁, ¼, n_d) describes properties of the configuration that are sampled, p(b, n₁,¼, n_d) is the weight of that configuration, and the histogram g(n₁,¼, n_d) counts how many configurations there are with the same properties (n₁,¼, n_d). In the classical simulations discussed above we had one-dimensional histiograms d = 1, the desired property was the energy n₁ = E, the weight the Boltzmann weight p(b, E) = exp(–bE) and the histogram g(E) the density of states.

Both histogram reweighting and Wang-Landau sampling are straightforwardly generalized to a representation (11). E.g., to employ Wang-Landau sampling, the weight p(b, n₁,¼, n_d) of a configuration in standard sampling should be replaced by 1/g(n₁,¼, n_d) and the algorithm will perform a random walk with a flat histogram in the space of the (n₁,¼, n_d).

We will now describe how the mapping to Eq. (11) can be done for quantum systems in three different representations: in discrete[13] and continuous time[14, 4] path integrals and in the stochastic series expansion (SSE) representation.[15]

3.2 Discrete Time Path Integrals

The historically first approach to quantum Monte Carlo simulations of lattice models was based on the Trotter-Suzuki decomposition, by splitting the Hamilton operator H into two or more terms H = H_i so that the exponentials of each of the terms exp(–bH_i) are easy to calculate. Although the H_i do not commute, the error in estimating the exponential

is small for small prefactors e. Applying this approximation to the partition function we get Suzuki's famous mapping, here shown for the simplest case of a chain with nearest neighbor interactions with two terms H₁ collecting terms on even bonds and H₂ terms on odd bonds:

where the time step is Dt = b/M, the |i_kñ each are complete orthonormal sets of basis states, and the transfer matrices are U_i = exp(–DtH_i). The evaluation of the matrix elements ái|U₁|i¢ñ is straightforward since the H_i are chosen to be easily diagonalized.

To see how this can be mapped to Eq. (11) let us consider a spin-1/2 model with Hamilton operator

Since each of the H_i is a sum of commuting bond terms, we can focus on a single bond term:

the exponential of which is

with

The product of matrix elements in Eq. (13) can be expressed as integer powers of these matrix elements:

with

Using a three-dimensional histogram instead of a one-dimensional one all classical reweighting approches, such as histogram reweighting, [10] multicanonical sampling,[11] and the Wang-Landau algorithm[8] can now be applied to a quantum system.

3.3 Continuous Time Path Integrals

Simulations in discrete imaginary time with finite time steps Dt require extrapolations to continuous time Dt ® 0 from calculations performed at several discrete time steps Dt. For models with a finite number of states per lattice site, the limit Dt ® 0 can already be taken in the construction of the algorithm and simulations can be performed directly in the continuous time limit Dt ® 0. Following Ref. [14], the continuous time limit can be formulated through time-dependent perturbation theory in imaginary time:

where the symbol denotes time-ordering of the exponential. The Hamiltonian H = H₀ + V is split into a diagonal term H₀ which in our example of a spin-1/2 model is

and an offdiagonal perturbation V:

The time-dependent perturbation in the interaction representation is V(t) = exp(tH₀)Vexp(–tH₀). In order to gain an expression in the form of Eq. (11) we insert sums over complete sets of basis states and make the b dependence explicit by redefining t _i ® bt_i:

Making use of the observation that all matrix elements of V are J_xy/2 we end up with an expression

where E(c) is the averaged diagonal energy of the configuration c = (i₁,¼, i_n):

and for the special case n = 0

We end up with Eq. (25), which is of the desired form of Eq. (11) with a two-dimensional histogram g(E, n) of the diagonal energies E and expansion orders n, and a weight function

At this point we need to address two technical issues before we are able to implement histogram reweighting or Wang-Landau sampling. The first is that the energy E is a continuous variable and we need to discretize it by introducing narrow energy bins of finite width DE. The second is that the sum over order n runs to infinity and needs to be truncated at some upper limit L. This truncation restricts the validity of the histogram reweighting or Wang-Landau sampling to inverse temperatures b < O(L/N) where N is the number of lattice sites. By choosing L large enough we can still reach down to any desired temperature.

3.4 Stochastic Series Expansion

In going from discrete to continuous time we not only got rid of the time discretization error but also simplified the histogram from a three-dimensional to a two-dimensional one. We can formally do this by setting H₀ = 0 and V = H in the above continuous time representation, which will result in a constant E = 0 but at the cost of now perturbing in the diagonal part of the Hamiltonian (the J_z term) which was previously treated exactly. After integrating out the time integrals we obtain the stochastic series expansion (SSE) representation:[15] We start by expressing the partition function as a high temperature expansion

where Z(0) = g(0) is just the number of states in the Hilbert space and the weight function is (assuming for simplicity a Heisenberg model with J_xy = J_z = J on a bipartite lattice)

The n-th order series coefficient g(n) now plays the role of the density of states in the classical algorithm and we can apply histogram reweighting or Wang-Landau sampling.

4 Application examples

As a first example of the efficiency of quantum Monte Carlo algorithms, we show in Fig. 1 Binder cumulant ratios of a three-dimensional (3D) quantum Heisenberg ferromagnet obtained from histogram reweighting of multiple histograms obtained in an SSE simulation. This demonstrates that quantum Monte Carlo simulations can give results with accuracy comparable to classical simulations and allow high precision determination of critical behavior.

To demonstrate Wang-Landau sampling for quantum systems, we first consider a second order thermal phase transition in the Heisenberg antiferromagnet on a simple cubic lattice. From simulations of systems with L³ sites, L = 4,6,8,12,16, we can calculate the staggered structure factor S(p, p, p) for any value of the temperature using the measured histograms. Fig. 2 shows the scaling plot of S(p, p, p)/L^2-h with h = 0.034. The estimate for the critical temperature T_c = 0.947 J, obtained in only a couple of days on an 800 MHz Pentium-III CPU, compares well with earlier estimates [16]. In Fig. 2, we show error bars as obtained from ten independent runs for each system size. Note that, compared to the histogram reweighting results in Fig. 1 we show results for smaller systems but over a much wider temperature range.

Next we demonstrate the efficiency of the algorithm at a first order phase transition by considering two-dimensional hardcore bosons with next-nearest neighbour interactions [17] described by:

where (a_i) creates (annihilates) a hardcore boson at site i, n_i is the density at this site, t the hopping amplitude between nearest neighbour sites, V₂ the next-nearest neighbour repulsion, and m the chemical potential. At low temperature and half filling this model is in an insulating phase with striped charge order and provides the simplest quantum mechanical model with stripes. Simulations with conventional update schemes suffer from exponentially increasing tunneling times needed to change the stripe orientation from a horizontal to a vertical arrangement. These tunneling times quantify the autocorrelation effects within a single simulation. The flat histogram in the order n in our algorithm reduces the tunneling times by many orders of magnitude already on small lattices (c.f. Fig. 3) which demonstrates the efficiency of our algorithm at first order phase transitions.

As a last example we show in Fig. 4 magnetization curves as a function of temperature for a spin ladder model.[18]

5 Coupling expansion

Histogram reweighting and Wang-Landau sampling can not only be applied to reweighting in the temperature but also to reweighting in the strength of a coupling constant. We will show as an example the SSE version.[9] Defining the Hamiltonian as H = H₀ + lV we can rewrite the partition function Eq. (29) as

where on the right hand side we have collected all terms associated with into (n_l).

We consider as a first example the quantum phase transition in a bilayer Heisenberg antiferromagnet whose ground state changes from quantum disordered to Néel ordered as the ratio l = J/J¢ of intra-plane (J) to inter-plane (J¢) coupling is increased [19]. From the histograms generated within one simulation we can calculate the staggered structure factor S(p, p) of the system at any value of l. In Fig. 5 we show a scaling plot of S(p, p) / L^2–z–h as a function of l. In short simulations, taking only a few days on an 800 MHz Pentium-III CPU, we find the quantum critical point at l = 0.396, which again compares well with earlier results [19]. Statistical error bars are obtained from ten independent runs, and shown in Fig. 5 at selected values of J¢/J.

As a final example we show magnetization curves for a spin flop transition in a two-dimensional quantum Heisenberg antiferromagnet in Fig. 6.

6 Conclusions

To summarize, we have shown how classical histogram reweighting and Wang-Landau sampling can be applied to quantum systems, allowing us to simulate quantum systems with similar accuracy as is possible for classical ones. This is possible in any representation, such as discrete or continuous time path integrals or in a stochastic series expansion (SSE) representation. At first order phase transitions we see, like in the classical case, a speedup of simulations by many orders of magnitude. These algorithms open up new possibilities for quantum Monte Carlo simulations.

We thank D.P. Landau for discussions and acknowledge support of the Swiss National Science Foundation. MT is grateful to the organizers of the third Brazilian Meeting on Simulational Physics for the invitation to visit their country.

References

[1] M. Suzuki et al., Prog. Theor. Phys. 58, 1377 (1977); J.E. Hirsch et al., Phys. Rev. Lett. 47, 1628 (1981).

[2] R.H. Swendsen and J.-S. Wang, Phys. Rev. Lett. 58, 86 (1987); U. Wolff, Phys. Rev. Lett. 62, 361 (1989).

[3] H.G. Evertz et al., Phys. Rev. Lett. 70, 875 (1993).

[4] B. B. Beard and U. J. Wiese, Phys. Rev. Lett. 77, 5130 (1996).

[5] A.W. Sandvik, Phys. Rev. B 59, R14157 (1999); O. F. Syljuåsen and A. W. Sandvik, Phys. Rev. E 66, 046701 (2002).

[6] N. V. Prokof'ev et al., Sov. Phys. - JETP 87, 310 (1998).

[7] N. Prokof'ev and B. Svistunov, Phys. Rev. Lett. 87, 160601 (2001).

[8] F. Wang and D.P. Landau, Phys. Rev. Lett. 86, 2050 (2001); Phys. Rev. E 64, 056101 (2001).

[9] M. Troyer, S. Wessel and F. Alet, Phys. Rev. Lett. 90, 120201 (2003).

[10] A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 61, 2635 (1988); ibid. 63, 1195 (1989).

[11] B.A. Berg and T. Neuhaus, Phys. Rev. Lett. 68, 9 (1992).

[12] P. Dayal et al., Preprint cond-mat/0306108.

[13] M. Suzuki, Prog. Theor. Phys. 65, 1454 (1976).

[14] N.V. Prokof'ev et al., Pis'ma v Zh.Eks. Teor. Fiz., 64, 853 (1996) [English translation is Report cond-mat/9612091].

[15]A.W. Sandvik and J. Kurkijärvi, Phys. Rev. B 43, 5950 (1991); A. W. Sandvik, J. Phys. A 25, 3667 (1992).

[16] A.W. Sandvik, Phys. Rev. Lett. 80, 5196 (1998).

[17] F. Hebert et al., Phys. Rev. B 65, 014513 (2002).

[18] T. Barnes et al., Phys. Rev. B 47, 3196 (1993); M. Troyer et al., Phys. Rev. B 50, 13515 (1994).

[19] A. W. Sandvik and D.J. Scalapino, Phys. Rev. Lett. 72, 2777 (1994).

Received on 2 September, 2003

Publication Dates

• Publication in this collection
  01 Sept 2004
• Date of issue
  June 2004

History

• Received
  02 Sept 2003
• Accepted
  02 Sept 2003