# Minimal energy ensemble Monte Carlo for the partition function of fermions coupled to classical fields

## Abstract

Models of non-interacting fermions coupled to auxilliary classical degrees of freedom are relevant to the understanding of a wide variety of problems in many body physics, e.g. the description of manganites, diluted magnetic semiconductors or strongly interacting electrons on lattices. Monte Carlo sampling over the classical fields is a powerful, yet notoriously challenging, method for this class of problems – it requires the solution of the fermion problem for each classical field configuration. Conventional Monte Carlo methods minimally utilize the information content of these solutions by extracting single temperature properties. We present a flat-histogram Monte Carlo algorithm that simulates a novel statistical ensemble which allows to acquire the full thermodynamic information, i.e. the partition function at all temperatures, of sampled classical configurations.

###### pacs:

05.10.Ln, 02.70.SsIntroduction - Bilinear Hamiltonians of lattice fermions coupled to classical degrees of freedom (continuous or discrete) are ubiquituous in contemporary many-body physics. These often arise as suitable approximations in the description of systems where many different degrees of freedom contrive to yield complex and interesting physics. In these cases, some subsystem is treated classically as e.g. localized spins in double-exchange models Zener (); Anderson (); deGennes (), models of Mn-doped (III, IV) semiconductors Schliemann (), the Ising t-J model IsingtJ (), adiabatic phonons in polaron models Holstein (); Millis (); Freericks1 (); Freericks2 (), or one species of fermions in the Falicov-Kimball model Freericks2 (); Hubbard (); FK (). Exactly solvable models can also take this form, such as the seminal honeycomb lattice Kitaev model Kitaev () where interacting spins are mapped onto Majorana fermions coupled to static gauge fields. More generally, auxiliary field methods, e.g. based on the Hubbard-Stratonovich transformation, allow to decouple fermion-fermion or fermion-boson interactions – the fields are then treated classically in conjunction with the application of the Suzuki-Trotter decomposition BSS (); HirschFye () (see Mukherjee14a () for a recent approximate scheme with static fields).

The simplicity of the form of such models belies their complexity. Although obtaining eigenstates for fixed configurations of classical fields is computationally easy, summing over the exponential number of such configurations to obtain thermodynamic properties is notoriously expensive. An obvious approach towards this problem is via Monte Carlo (MC) sampling Ulam (). Summing over the fermion states for fixed classical field configuration yields the conditional (grand) partition function. This quantity, at fixed temperature, serves as the Metropolis weight Metropolis () for performing a random walk in the space of classical field configurations Dagotto-review (). The serious bottleneck in these simulations is the repeated performance of the fermionic trace – and so exact diagonalization ED of the free fermion system – for executing the walk. Hence various improvements have been proposed to optimize the reevaluation of the weight - moment expansion of the fermion density of states by Chebyshev polynomials KPM (); TKPM (), low-rank matrix updates LRMUp () or Green’s functions GFKMP () Chebyshev expansion. On the other hand it seems naturally essential to optimize the extraction of information at each MC step. Indeed, while the expensive ED yields the conditional partition function at all temperatures, these are completely discarded by using only single temperature data in the above approaches. Here, we introduce an algorithm that fully exploits the thermodynamic information available at each diagonalization step in MC simulations of free fermions coupled to classical fields (FCCF).

The paradigmatic Metropolis algorithm Metropolis () suffers from critical slowing down at continuous phase transitions and prolonged trapping in metastable states at discontinuous transitions. These can be overcome to a large extent using cluster update schemes Swendsen (); Wolff (); loop (); worm (); loopoper (); directloop () or sampling extended ensembles singlehist (); multican (); broad (); paralel () such as Wang-Landau sampling of the density of states WL (); QWL (). While these approaches have been used in the study of classical and quantum systems, FCCF hold their own system-specific challenges rendering such applications difficult. In particular, the effective Hamiltonian (corresponding to the energies in the Metropolis weight) of the classical fields in general contains temperature dependent, long-range, multi-particle interactions. Often, molecular dynamics (MD) or hybrid Monte Carlo methods using Langevin’s equation hybrid1 (); hybrid2 (); hybrid3 (); molecular () are better suited than standard MC simulations of such Hamiltonians although the acceptance rate in these simulations crucially depends on the quality of the approximate action. Instead, we sample a novel extended ensemble bringing the advantages of Wang-Landau-like sampling to FCCF.

The problem and method - We first generalize our considerations: a system is bi-separable if it may be separated into two subsystems and such that for a given state of all states of the system can be efficiently summed over to obtain the conditional (grand) partition function or equivalently the conditional free energy (grand potential) . This definition covers both FCCF (where subsystem is classical) and many classical models (as e.g. the Ising model on a bipartite lattice where and are the spins on the two sublattices respectively). For bi-separable systems, the partition function is decomposed as:

(1) |

i.e. the partition function is obtained by averaging the conditional partition functions over all configurations attainable in the system . Notice that once the conditional energy spectrum is obtained, complete thermodynamic information associated with the exponentially large number of configurations of encoded in for all inverse temperatures becomes potentially available at each simulation step.

Our problem can be stated as follows. Standard MC sampling for FCCF systems consists of walking, at fixed temperature, between different configurations of with with Metropolis weights . Here, we aim to obtain the entire partition function (1) by acquiring, at a given simulation step, the conditional partition function for all arguments from the information (full spectrum) abundantly available for each fixed configuration on . In principle this may also be achieved by parallel tempering (to efficiently sample configurational space of ) and a reweighting procedure, which however requires well chosen set of temperature intervals. Instead, in our method, we perform a random walk directly in the configurational space of without referring to any specific temperature. The basic challenge here is to obtain an appropriate importance sampling scheme over the (exponential number of) configurations of .

The key behind any thermodynamic Monte Carlo simulation lies in sweeping through energy space efficiently. A simple observation reveals that, for most systems, only a few configurations of will lead to the spectrum of system energies containing the ground state energy. These configurations of obviously must be effectively found by the importance sampling scheme. On the other hand the conditional DOSes of energy spectra associated with typical configurations of are expected to differ most also in their lower range. Therefore a key discriminant for configurations of is the minimal energy attainable by system at a given configuration which we will denote as . We consider two configurations of as belonging to the same class when they have equal values of .

Formally, any importance sampling scheme can be represented by assigning a weight function to the set of configurations of subsystem , such that

(2) |

The principle that we identify and implement for acquiring the partition function is that all minimal energy classes be visited by the algorithm. For this, notice that the minimal attainable energies can be associated to a density of states – , which enumerates the number of configurations of subsystem attaining a given minimal energy. We emphasize that this ”auxiliary” DOS is distinct from the true DOS of the system, and does not determine the latter.

Our algorithm consists of two separate sampling stages associated with first generating the weight distribution for importance sampling and then utilizing it to acquire thermodynamic information about the system. (i) The auxiliary DOS is readily obtained by performing a Wang-Landau (WL) algorithm WL () in the space of minimal energies, where plays the role of ”energy” of a given configuration of . (ii) Next, one performs a random walk in the space of configurations of with weight function

(3) |

to sample (Eq.(2)) for all . This choice of weight function not only allows to visit all classes of configurations but visits each class approximately the same number of times yielding a flat histogram.

The above principle can be viewed as a minimal necessary requirement for sampling energy space. It yields a coarse-grained view of energy space disregarding any differentiation between configurations belonging to a given class. Below, we show that the application of this principle leads to remarkably accurate results.

Testing the algorithm - The convergence properties of the WL algorithm which constitutes the first stage of our algorithm have been widely studied. We note, that the random walk in the second stage of our algorithm (with fixed auxiliary DOS) fulfills, by fiat, detailed balance. We test the convergence of our algorithm on the (a) Ising and (b) Potts model. Finally, we present continuous temperature results for the Falicov-Kimball model as a prominent example of FCCF.

The Ising model with nearest neighbour interactions on a square lattice is an ideal benchmark for testing new algorithms since both an exact solution in the thermodynamic limit as well as large system high precision Monte Carlo results are available. This model is bi-separable which allows for direct application of our algorithm. Indeed for given configuration of spins on one of the sublattices (subsystem ) the conditional partition function corresponding to all configurations of spins from the second sublattice (subsystem ) is simply:

(4) |

where is the coupling between spins, is the number of sites on sublattice . and are the number of spins on sublattice subject to the nearest neighbour fields with absolute value and respectively. Of course, and are dictated by the configuration of spins on sublattice . Fig. 1 shows the continuous temperature dependence of the specific heat, which is the temperature derivative of the free energy, obtained with our and Wang-Landau algorithms. These are generated for a 60 x 60 site lattice (with periodic boundary conditions PBC). The two curves show good agreement, indicating sufficiency of our prescribed importance sampling principle. Any essential problems regarding the effectiveness of this sampling scheme would have been expected to be apparent for this moderately large lattice size.

We now highlight some important distinctions between the WL and our algorithms. A single step in the second stage of our algorithm entails accumulation of full thermodynamic information from all configurations of subsystem (a remarkable configurations in the presented simulation) during each update move on subsystem . In contrast the standard WL algorithm updates information from one configuration per move. However this ”exponential update” comes at the initial cost of first obtaining the ”auxiliary” DOS via a WL procedure. Interestingly, the ”auxiliary” DOS has to be determined with essentially a higher histogram flatness requirement than the system DOS in the direct WL simulation of the Ising model. This has two sources: (i) The ”subsystem Hamiltonian” with energies contains multi-spin and long range interactions unlike the original Ising interaction Hamiltonian. (ii) Sensitivity to the precision of the ”auxiliary” DOS. We have found that adding small fluctuations to rapidly deteriorates results. Hence, our algorithm, should not be viewed as an alternative to the standard WL algorithm for classical models. Finally while the WL algorithm directly outputs the system DOS, from which the partition function may be easily calculated, our algorithm yields the partition function. This difference has important practical consequences. In the WL algorithm the DOS values (which generally grow exponentially with system size) are accumulated in a given run by multiplying them with the so-called modification factor every time a given energy is visited. In subsequent runs is gradually reduced to 1. This ingenious trick allows to build up very big DOS values in a reasonable number of steps). In our procedure we sum the accumulated quantities , and the values of both factors are already very big or very small. Therefore care must be taken to avoid roundoff error acumulation during summation.

To illustrate that our algorithm succesfully simulates discontinuous phase transitions, we consider the 10-state Potts model on the square lattice with nearest neighbour interaction. Bi-separability here may be shown in similar fashion as in the Ising model. However the form of the conditional partition function is more complicated than (4) due to the multitude of values of local fields set by configurations of nearest neighbour Potts spins. Comparison in Fig. 2 of the specific heat, obtained by the WL and our algorithm, for a 30 x 30 – site lattice with PBC reveals the effectiveness of our sampling scheme in this case as well.

Now, consider the Falicov-Kimball Hamiltonian (FKH) of FCCF:

(5) |

where and are creation operators of mobile and immobile fermions respectively, , , is the nearest neighbour hopping integral, is the Hubbard on-site interaction and is the chemical potential for fermions (the number of immobile fermions is fixed). For a given configuration of immobile fermions , the set of single particle eigenenergies of mobile fermions is readily obtained, rendering the model bi-separable.

Unlike the classical models above, an efficient standard WL algorithm is not directly applicable to the simulation of the FKH. We consider the FKH on a square lattice under PBC with , at half filling (), for different values of . Under these conditions, the FK model undergoes, at low temperatures, a transition to the charge density wave (CDW) ordered state, with ordering wave-vector. The transition is continuous for large with some evidence pointing to a change to discontinous transitions for small Maska (); Zonda (). The application of our algorithm yields all-temperature results unlike standard Metropolis sampling over the immobile fermions. Moreover, since the Metropolis algorithm suffers from slow kinetics at discontinous phase transitions our method should be particulary useful in further investigations of the small regime. Here, we restrict to proving that our sampling principle works. In Fig. 3, a comparison of results for the specific heat and internal energy obtained by Metropolis sampling and our algorithm is presented. On all diagrams the results are in good agreement within the accuracy of the local update Metropolis algorithm, which again indicates the effectiveness of our sampling. We mention here, that in obtaining the ”auxiliary” DOS in the first stage of our algorithm, we discretized values of minimal energies of the total system (but we have not discretized single particle energies anywhere else) by assigning them to (energy) bins of width or and checked convergence.

Summary and conclusions - In this Letter, we have presented a new Monte Carlo algorithm for problems of fermions coupled to classical degrees of freedom. This algorithm is based on Wang-Landau-like sampling in contradistinction to the commonly used Metropolis sampling. This allows in principle to overcome drawbacks of Metropolis schemes and importantly to fully exploit all available thermodynamic information at each diagonalization step – information that is mostly wasted in other MC schemes. The scheme is based on the notion of minimal energy attainable for a given classical field configuration. As a minimal requirement, we devised a rule that all such minimal energies be visited by the algorithm. This was achieved by first obtaining an auxiliary DOS for minimal energies via Wang-Landau sampling, and then using its inverse as the weight function to perform a random walk in classical field configurations accumulating the full partition functions. We benchmarked our recipe for several paradigmatic models of statistical physics achieving excellent agreement with known results. We mention here that while our principle of sampling minimal energies yields satisfactory results in the presented examples, more generally, supplementary conditions may need to be identified in other cases. However, our results show that in principle accumulation of conditional partition functions at all temperatures at once using simple, temperature independent importance sampling is possible. An intriguing possibility is the application of this method to problems of quenched disorder, where free energies need to be efficiently averaged over the disorder realizations.

Finally, we comment on the most distinctive feature of our algorithm – i.e. accumulation of full conditional partition functions per update. Standard algorithms update information from only single configurations per move. However, this is by no means the only possibility. N-fold way nfold () algorithms may be seen as updating information from configurations during each move (where is the number of system sites). Our algorithm is remarkable in that it uniquely allows the update of information from an exponential number of configurations during each move ( for some constant ). Clearly, no algorithms exist that can update more information from a complexity point of view.

###### Acknowledgements.

Acknowledgements - We thank Mathias Troyer for interesting and useful comments. The authors acknowledge support from the (Polish) National Science Center Grant No DEC-2011/03/B/ST2/01903.### References

- C. Zener, Phys. Rev. 82, 403 (1951).
- P. W. Anderson and H. Hasegawa, Phys. Rev. 100, 675 (1955).
- P. G. de Gennes, Phys. Rev. 118, 141 (1960).
- J. Schliemann, J. Konig, and A. H. MacDonald, Phys. Rev. B 64, 165201 (2001).
- M. M. Maśka, M. Mierzejewski, A Ferraz and E A Kochetov, J. Phys.: Condens. Matter 21 045703 (2009).
- T. Holstein, Ann. Phys., NY 8 325, (1959).
- A. J. Millis, R. Mueller and B. I. Shraiman, Phys. Rev. B 54 5389 (1996).
- J. K. Freericks, M. Jarrell, and G. D. Mahan, Phys. Rev. Lett. 77, 4588 (1996).
- J. K. Freericks and V. Zlatić, Rev. Mod. Phys. 75, 1333 (2003).
- J. Hubbard, Proc. R. Soc. London, Ser. A 276, 238 (1963); 281, 401 (1964).
- L. M. Falicov and J. C. Kimball, Phys. Rev. Lett. 22, 997 (1969).
- A. Kitaev, Annals of Physics 321,2 (2006)
- R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Phys. Rev. D 24, 2278 (1981).
- J. E. Hirsch and R. M. Fye, Phys. Rev. Lett. 56, 2521 (1986).
- A. Mukherjee, N. D. Patel, S. Dong, S. Johnston, A. Moreo, and E. Dagotto, Phys. Rev. B 90, 205133 (2014).
- S. Ulam, R. D. Richtmyer, and J. von Neumann. 1947. Statistical methods in neutron diffusion. Los Alamos Scientific Laboratory report LAMS-551.
- N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, J. Chem. Phys. 21, 1087 (1953).
- E. Dagotto, Nanoscale Phase Separation and Colossal Magnetoresistance: The Physics of Manganites and Related Compounds (Springer-Verlag, Berlin, 2003).
- Y. Motome and N. Furukawa, J. Phys. Soc. Jpn. 68, 3853 (1999); ibid. 69, 3785 (2000); ibid. 70, 3186 (2001).
- Y. Motome and N. Furukawa, J. Phys. Soc. Jpn. 72, 2126 (2003); N. Furukawa and Y. Motome, ibid. 73, 1482 (2004).; G. Alvarez, C. Sen, N. Furukawa, Y. Motome, and E. Dagotto, Comput. Phys. Commun. 168, 32 (2005).
- G. Alvarez, P. K.V.V. Nukala, and E. D’Azevedo, J. Stat. Mech. P08007 (2007).
- A. Weiße, Phys. Rev. Lett. 102, 150604 (2009).
- R.H. Swendsen and J.-S. Wang, Phys. Rev. Lett. 58, 86 (1987).
- U. Wolff, Phys. Rev. Lett. 60, 1461 (1988).
- H.G. Evertz, G. Lana and M. Marcu, Phys. Rev. Lett. 70, 875 (1993).
- N.V. Prokof’ev, B.V. Svistunov and I.S. Tupitsyn, Sov. Phys. - JETP 87, 310 (1998).
- A.W. Sandvik, Phys. Rev. B 59, R14157 (1999); A. Dorneich and M. Troyer, Phys. Rev. E 64, 066701 (2001).
- O.F. SyljuĂĽsen and A.W. Sandvik, Phys. Rev. E 66, 046701 (2002); O.F. SyljuĂĽsen, Phys. Rev. E 67, 046701 (2003).
- A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 61, 2635 (1988).
- B. A. Berg and T. Neuhaus, Phys. Rev. Lett. 68, 9 (1992); Phys. Lett. B. 267, 249 (1991).
- P. M. C. de Oliveira et al., Braz. J. Phys. 26, 677 (1996).
- K. Hukushima and K. Nemoto, J. Phys. Soc. Japan 65, 1604 (1996); E. Marinari and G. Parisi, Europhys. Lett. 19, 451 (1992); A. P. Lyubartsev et al., J. Chem. Phys. 96, 1776 (1992).
- F. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 (2001); Phys. Rev. E 64, 056101 (2001).
- M. Troyer, S. Wessel and F. Alet, Phys. Rev. Lett. 90, 120201 (2003).
- R. T. Scalettar, D. J. Scalapino, and R. L. Sugar, Phys. Rev. B 34, 7911 (1986).
- S. Duane, A. Kennedy, B. J. Pendleton, and D. Roweth, Phys. Lett. B 195, 216 (1987).
- J. L. Alonso, L. A. Fernández, F. Guinea, V. Laliena, and V. Martín-Mayor, Nucl. Phys. B596, 587 (2001).
- K. Barros and Y. Kato, Phys. Rev. B 88, 235101 (2013).
- M.M. Maśka and K. Czajka, Phys. Rev. B 74, 035109 (2006).
- M. Žonda, P. Farkašovský, H. Čenčariková, Solid State Commun. 149, 45 (2009).
- A. B. Bortz, M. H. Kalos, and J. L. Lebowitz, J. Comput. Phys. 17, 10 (1975).