Numerical Simulation of Correlated Electron Systems
Correlated electron systems belong to the general class of complex systems whose defining property is that the whole is more than the sum of the individual parts. This emergent complexity gives rise to a rich variety of fascinating phenomena in fields ranging from biological systems to materials. At the same time, the complexity poses a challenging scientific problem: Even if we find a way to solve the corresponding equations of motion of the constituents of a many-body system, it is by no means guaranteed that we can understand the emergent collective phenomena.
In the case of crystalline materials, we can start with a quantum-mechanical description of the electronic motion in the periodic potential of the ions in terms of Bloch states. From these states, we can construct the many-body wavefunction by taking into account the Pauli exclusion principle that demands the wave function to be antisymmetric under exchange of identical fermions. These two ingredients, Bloch states and Fermi statistics, underlie band structure theory and produce a host of very interesting states, most notably insulators and metals with completely or partially filled highest occupied band, respectively. As discovered theoretically by Kane and Mele , band theory also permits more exotic insulating states referred to as topological insulators. These states, which arise from strong spin-orbit coupling and time-reversal symmetry, are characterized by insulating behaviour in the interior but metallic behaviour at the edge or surface of the sample. Importantly, the characteristic surface states are robust with respect to disorder, and can therefore be observed in experiments . In two dimensions, the edge states are helical, that is the direction of motion of electrons is coupled to the direction of their spin . In principle, this property and the topological protection make topological insulators a promising candidate for applications in the field of spintronics, which offers the vision of overcoming the problem of heating inherent to current CPU technology.
Our discussion so far has neglected electron-electron interactions. However, as students learn already in introductory courses on solid state physics, the energy scale associated with such correlations is very large, typically several electron volts. Nevertheless, in many cases, electronic correlations do not significantly alter the properties of the system – they are irrelevant in the language of Wilson’s renormalization group. This remarkable fact hinges on the screening of the Coulomb repulsion, and Fermi statistics which drastically limits the phase space available for scattering due to the Coulomb repulsion, and is known as Fermi liquid theory.
Materials in which Fermi liquid theory does not hold are of particular interest. In such systems, correlations lead to exotic and/or collective states of matter. In one dimension, the concept of the Fermi liquid breaks down and correlations cause a fractionalization of electrons into independent, collective spin and charge excitations. More generally, Mott insulators are correlated materials which are incorrectly predicted by band theory to be metallic. The insulating behaviour is often attributed to the observed ordering of, for example, the spin degrees of freedom (that is, the formation of a magnetic state). However, a central question in this context is if these two phenomena are always coupled. Remarkably, the answer to this question is no. Mott insulating states with no associated symmetry breaking can be realized in a number of models of correlated electron systems. In two dimensions, these states are yet another example of exotic states of matter which support anyonic excitations with fractional statistics. Mott insulators have also received considerable attention because they produce high-temperature superconductivity upon doping with electrons or holes. Recently, the possibility of a quantum spin liquid phase has been hotly debated for the case of the Hubbard model on the honeycomb lattice [3–5]. Motivated by the discovery of topological insulators, there have been numerous studies of the honeycomb Hubbard model with additional spin-orbit coupling, see  for a review. The Mott transition of fermions on the honeycomb lattice has recently been realized experimentally using cold atoms in an optical lattice .
The Numerical Challenge
For quantum many-body problems, the Hilbert space grows exponentially with the number of electrons, or the volume of the system. A brute force diagonalization of the full problem hence requires an exponential effort. The question is whether or not we have to consider all possible states to understand the problem. Here, we choose the stochastic approach, and use the auxiliary-field quantum Monte Carlo method . The method is based on a path integral formulation where the interaction is decomposed with the help of a HubbardStratonovitch field that mediates the electronic correlations. The integration over the field configurations is carried out with the Monte Carlo method. A generic problem of this approach is the sign problem (the appearance of configurations with negative weights), causing an exponential increase of computer time with volume and inverse temperature to reach a given accuracy. However, for a rather large class of non-trivial models [9–12], symmetries permit us to avoid the sign problem. For such models, the numerical effort scales as N 3 β where β is the inverse temperature and N the number of electrons. The method used here can be efficiently parallelized and run on modern architectures. Recently, various optimization schemes, in particular cache optimization, have allowed us to gain up to an order of magnitude in performance.
|Figure 1: The Kane-Mele-Hubbard model describes electrons on the honeycomb lattice with nearest-neighbour hopping t, spin-orbit coupling λ, and an onsite repulsion U between electrons occupying the same lattice site.
The Kane-Mele-Hubbard Model
Motivated by the experimental progress with graphene and related systems (such as silicene and artificial graphene), we consider the paradigmatic Hubbard model for strongly correlated electrons on the honeycomb lattice  illustrated in Fig. 1. Electrons can hop from site to site with hopping amplitude t, and two electrons at the same lattice site experience a repulsion U that captures the essence of Coulomb repulsion. This model has time reversal, SU(2) spin, and sublattice symmetry, and the noninteracting band structure is that of massless Dirac fermions. Adding a spin-orbit term in the form of a complex second-neighbour hopping λ  leads to the Kane-Mele-Hubbard model for correlated topological insulators . The spin-orbit term breaks sublattice symmetry and thereby opens a topological mass gap resulting in a quantum spin Hall insulator (a widely used name for two-dimensional topological insulators). At half filling (one electron per lattice site), the above symmetries permit us to apply the auxiliary-field quantum Monte Carlo method without a sign problem, and hence to investigate the role of electronic correlations accurately and without uncontrolled approximations.
|Figure 2: The zero-temperature phase diagram of the Kane-Mele-Hubbard model as obtained from quantum Monte Carlo simulations with up to 648 electrons [5, 11]. Here AFM indicates an antiferromagnetic Mott insulator with magnetic order either in the xyz or in the xy directions, respectively. The xyz AFM and semimetal phases exist only at λ = 0.
In the limits of weak and strong coupling, the low-energy physics of the Kane-Mele-Hubbard model can be under- stood quite easily. The semimetallic state of Dirac fermions at λ = U = 0 is stable with respect to correlations and hence expected to survive at small U > 0. Similarly, at U = 0, the spin-orbit coupling establishes a gapped quantum spin Hall state, and the bulk physics remains unaffected by small interactions U. For strong coupling U/t » 1, the charge degrees of freedom become frozen (any deviations from the state with one electron per site are suppressed by U) but the spin degrees of freedom remain active. Using perturbation theory in the small parameter t/U, one can derive effective spin models. In the absence of spin orbit coupling, the SU(2) spin symmetry leads to a Heisenberg model with antiferromagnetic nearest-neighbour exchange. Since the lattice is non-frustrated, we expect long-range antiferromagnetic order at zero temperature, and a gapless Goldstone mode corresponding to spin-wave excitations. This antiferromagnetic Mott insulator (AFM) state is labelled xyz AFM in Fig. 2. Spin-orbit coupling reduces the spin symmetry from SU(2) to U(1). The resulting spin model is more complex, and includes frustrated interactions in the z direction of spin . Long-range magnetic order can develop in the transverse direction, as is the case in the xy AFM phase in Fig. 2.
The validity of the above considerations can be verified using numerical simulations. Importantly, such methods also permit to study the intermediate-coupling regime and hence the evolution from small to large U/t. It is of particular interest whether or not there exists a direct transition between the weak- and strong-coupling phases, or if instead exotic intermediate phases appear. Moreover, if continuous quantum phase transitions occur, it is of great interest to identify their universality class. In the following, we summarize our present understanding of correlation effects in the Kane-Mele-Hubbard model, as it emerges from large-scale quantum Monte Carlo simulations [3, 5,10,11,14].
Impact of U for λ = 0
The following two scenarios can be envisaged for the transition from the semimetal at weak coupling to the antiferromagnetic Mott insulator at strong coupling. Starting from strong coupling, and noting that the insulator to semimetal transition is numerically found to occur at values of U smaller than the electronic bandwidth W = 3t, one can construct an effective spin model for the region close to the transition where higher-order ring-exchange terms proliferate with decreasing U/t . This point of view suggests that the melting of magnetic order may be independent of the metal-insulator transition. Quantum Monte Carlo simulations  suggested the existence of an intermediary spin liquid phase, with a single-particle gap but no long-range magnetic order, that separates the semimetal from the magnetic insulator. Similar conclusions have been drawn for the related π-flux model on the square lattice . The results of  have been challenged by more recent studies: Entropy calculations do not favour a degenerate ground state as expected for the Z2 spin liquid , and  demonstrates that the use of significantly larger system sizes leads to an almost complete disappearance of the purported spin-liquid from the phase diagram.
|Figure 3: Data collapse of the staggered magnetic moment m at λ = 0, measured for different system sizes L. Here we have used a critical value Uc/t = 3.78 , and critical exponents obtained from an ε-expansion and Gross-Neveu-Yukawa theory [15,16]. Results taken from .
Alternatively, we can start from the weakcoupling semimetalic state. Upon increasing U/t, we eventually expect a phase transition to an insulating state with antiferromagnetic order [15,20, 21]. Because gapless, linear fermionic excitations exist below the transition, the transition is expected to fall into the Gross-Neveu universality class [15,16]. This expectation is based on Gross-Neveu-Yukawa theory [5,15,16], which describes massless Dirac fermions with a Yukawa coupling to a vector boson describing the bosonic Goldstone mode of the antiferromagnetic state. The fermionic mass (the gap of the insulating state) is generated by the condensation of the bosonic mode (antiferromagnetic order). The theory has an upper critical dimension of d = 3 where the Gaussian approximation becomes exact. Within a first-order ϵ-expansion, one can estimate the critical exponents for the present twodimensional case. These theoretical predictions are confirmed by recent quantum Monte Carlo simulations . Using the predicted critical exponents, we observe an excellent data collapse [see Fig. 3] and an identical finite-size scaling of both the single-particle gap and the staggered magnetization . Importantly, Gross-Neveu theory predicts that the long-ranged tail of the Coulomb repulsion, relevant for possible future experiments, will not alter the nature of the transition. Numerical work to confirm this hypothesis is currently in progress.
Impact of U for λ > 0
The transition from the quantum spin Hall phase to the xy antiferromagnetic Mott insulator is an insulator to insulator transition. Because the onset of magnetic order breaks the time-reversal symmetry that protects the quantum spin Hall state, a direct transition is possible. Since the system is insulating on both sides of the transition, the latter involves only the spin degrees of freedom which order in the xy plane at the critical point. Such a transition, involving a U(1) order parameter in two spatial dimensions, is expected to fall into the same universality class as the three-dimensional XY model. Using the quantum Monte Carlo method, we have successfully verified this prediction in .
|Figure 4: Integrated spectral weight of the spin fluxon states created by inserting a pair of π fluxes into the honeycomb lattice. A response occurs in the correlated topological state [see Fig. 2], but not in the topologically trivial antiferromagnetic insulator. Here, λ = 0.2t. Results taken from .
An important question in the context of topological insulators is how a correlated topological state can be detected numerically . One attractive possibility, demonstrated in , is to exploit a unique property of quantum spin Hall insulators, namely the quantization of the spin Hall conductivity . As a consequence, the insertion of a magnetic flux of the size of half a flux quantum (a π flux) gives rise to two types of states localized around the flux: A so-called Kramers pair of spin fluxons carrying spin ±1/2, and a pair of charge fluxons carrying charge ±e [22, 23]. The existence of these states and their quantum numbers are directly related to the topological invariant [22, 23]. Electronic interactions remove charge fluxons from the low-energy excitation spectrum of the system, leaving spin fluxons which can be detected via a Curie-law-type contribution to the bulk magnetic susceptibility . Alternatively, as shown in Fig. 4, it is possible to measure the spectral weight of spin fluxon excitations in real space. The combination of π fluxes and quantum Monte Carlo simulations provides an efficient tool to identify correlated quantum spin Hall states.
Quantum Monte Carlo simulations of the Hubbard model with additional spin-orbit coupling allow us to study the interplay between band topology and correlation effects. Our results on finite lattices are free from systematic errors and therefore also serve as benchmarks for approximate theories. The continuous algorithmic progress enables us to address exciting open questions, such as the effect of a longrange Coulomb interaction on the phase diagram and the phase transitions.
Parts of this work were done in collaboration with M. Bercx, F. Goth, I. Herbut, T. Lang, Z.Y. Meng, A. Muramatsu and S. Wessel. Funding from the DFG under the grant number AS 120/9-1, AS120/10-1, and Ho 4489/2-1 (Forschergruppe FOR 1807) is acknowledged. We thank the Jülich Supercomputing Centre for generous allocation of CPU time.
 Kane, C.L., Mele, E.J. E.J.
Phys. Rev. Lett. 95, 146802, 2005
 König, M., Wiedmann, S., Brüne, C., Roth, A., Buhmann, H., Molenkamp, L.W.,Qi, X.-L., Zhang, S.-C.-L., Zhang, S.-C.
Science 318, 766, 200
 Meng, Z.Y., Lang, T.C., Wessel, S., Assaad, F.F., Muramatsu, A.
Nature 464, 847, 2010
 Sorella, S., Otsuka, Y., Yunoki, S.
Sci. Rep. 2, 992, 2012
 Assaad, F.F., Herbut, I.F.
Phys. Rev. X 3, 031010, 2013
 Hohenadler, M., Assaad, F.F.
J. Phys.: Condens. Matter 25, 143201, 2013
 Uehlinger, T., Jotzu, G., Messer, M., Greif, D., Hofstetter, W., Bissbort, U., Esslinger, T.
 Assaad, F.F., Evertz, H.G.
Computational Many Particle Physics, Lecture Notes in Physics, Vol. 739, edited by H. Fehske, R. Schneider, and A. Weiße (Springer Verlag, Berlin, 2008) pp. 277–356
 Assaad, F.F.
Phys. Rev. Lett. 83, 796 (1999)
 Hohenadler, M., Lang, T.C., Assaad, F.F.
Phys. Rev. Lett. 106, 100403 (2011)
 Hohenadler, M., Meng, Z. Y., Lang, T. C., Wessel, S., Muramatsu, A., Assaad, F. F.
Phys. Rev. B 85, 115132, 2012
 Berg, E., Metlitski, M.A., Sachdev, S.
Science 338, 1606, 2012
 Rachel, S., Le Hur, K.
Phys. Rev. B 82, 075106, 2010
 Assaad, F.F., Bercx, M., Hohenadler, M.
Phys. Rev. X 3, 011015, 2013
 Herbut, I.F.
Phys. Rev. Lett. 97, 146401, 2006
 Herbut, I.F., Juričić V., Vafek, O.
Phys. Rev. B 80, 075432, 2009
 Yang, H.-Y., Albuquerque, A.F., Capponi, S., Läuchli, A.M., Schmidt, K.P.
New Journal of Physics 14, 115027, 2012
 Chang, C.-C., Scalettar, R.T.
Phys. Rev. Lett. 109, 026404, 2012
 Clark, B.K.
 Sorella, S., Tosatti, E.
Europhys. Lett. 19, 699, 1992
 Paiva, T., Scalettar, R.T., Zheng, W., Singh, R.R.P., Oitmaa, J.
Phys. Rev. B 72, 085123, 2005
 Qi, X.-L., Zhang, S.-C.
Phys. Rev. Lett. 101, 086802, 2008
 Ran, Y., Vishwanath, A., Lee, D.-H.
Phys. Rev. Lett. 101, 086801, 2008
• Fakher F. Assaad