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 manybody 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 quantummechanical description of the electronic motion in the periodic potential of the ions in terms of Bloch states. From these states, we can construct the manybody 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 [1], band theory also permits more exotic insulating states referred to as topological insulators. These states, which arise from strong spinorbit coupling and timereversal 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 [2]. In two dimensions, the edge states are helical, that is the direction of motion of electrons is coupled to the direction of their spin [1]. 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 electronelectron 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 hightemperature 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 spinorbit coupling, see [6] for a review. The Mott transition of fermions on the honeycomb lattice has recently been realized experimentally using cold atoms in an optical lattice [7].
The Numerical Challenge
For quantum manybody 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 auxiliaryfield quantum Monte Carlo method [8]. 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 nontrivial 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 KaneMeleHubbard model describes electrons on the honeycomb lattice with nearestneighbour hopping t, spinorbit coupling λ, and an onsite repulsion U between electrons occupying the same lattice site. 
The KaneMeleHubbard 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 [3] 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 spinorbit term in the form of a complex secondneighbour hopping λ [1] leads to the KaneMeleHubbard model for correlated topological insulators [13]. The spinorbit term breaks sublattice symmetry and thereby opens a topological mass gap resulting in a quantum spin Hall insulator (a widely used name for twodimensional topological insulators). At half filling (one electron per lattice site), the above symmetries permit us to apply the auxiliaryfield 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 zerotemperature phase diagram of the KaneMeleHubbard 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. 
Results
In the limits of weak and strong coupling, the lowenergy physics of the KaneMeleHubbard 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 spinorbit 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 nearestneighbour exchange. Since the lattice is nonfrustrated, we expect longrange antiferromagnetic order at zero temperature, and a gapless Goldstone mode corresponding to spinwave excitations. This antiferromagnetic Mott insulator (AFM) state is labelled xyz AFM in Fig. 2. Spinorbit 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 [13]. Longrange 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 intermediatecoupling 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 strongcoupling 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 KaneMeleHubbard model, as it emerges from largescale 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 higherorder ringexchange terms proliferate with decreasing U/t [17]. This point of view suggests that the melting of magnetic order may be independent of the metalinsulator transition. Quantum Monte Carlo simulations [3] suggested the existence of an intermediary spin liquid phase, with a singleparticle gap but no longrange magnetic order, that separates the semimetal from the magnetic insulator. Similar conclusions have been drawn for the related πflux model on the square lattice [18]. The results of [3] have been challenged by more recent studies: Entropy calculations do not favour a degenerate ground state as expected for the Z2 spin liquid [19], and [4] demonstrates that the use of significantly larger system sizes leads to an almost complete disappearance of the purported spinliquid 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 [5], and critical exponents obtained from an εexpansion and GrossNeveuYukawa theory [15,16]. Results taken from [5]. 
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 GrossNeveu universality class [15,16]. This expectation is based on GrossNeveuYukawa 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 firstorder ϵexpansion, one can estimate the critical exponents for the present twodimensional case. These theoretical predictions are confirmed by recent quantum Monte Carlo simulations [5]. Using the predicted critical exponents, we observe an excellent data collapse [see Fig. 3] and an identical finitesize scaling of both the singleparticle gap and the staggered magnetization [5]. Importantly, GrossNeveu theory predicts that the longranged 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 timereversal 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 threedimensional XY model. Using the quantum Monte Carlo method, we have successfully verified this prediction in [11].

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 [14]. 
An important question in the context of topological insulators is how a correlated topological state can be detected numerically [6]. One attractive possibility, demonstrated in [14], is to exploit a unique property of quantum spin Hall insulators, namely the quantization of the spin Hall conductivity [1]. 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 socalled 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 lowenergy excitation spectrum of the system, leaving spin fluxons which can be detected via a Curielawtype contribution to the bulk magnetic susceptibility [14]. 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.
Outlook
Quantum Monte Carlo simulations of the Hubbard model with additional spinorbit 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.
Acknowledgments
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/91, AS120/101, and Ho 4489/21 (Forschergruppe FOR 1807) is acknowledged. We thank the Jülich Supercomputing Centre for generous allocation of CPU time.
References
[1] Kane, C.L., Mele, E.J. E.J.
Phys. Rev. Lett. 95, 146802, 2005
[2] 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[3] Meng, Z.Y., Lang, T.C., Wessel, S., Assaad, F.F., Muramatsu, A.
Nature 464, 847, 2010
[4] Sorella, S., Otsuka, Y., Yunoki, S.
Sci. Rep. 2, 992, 2012
[5] Assaad, F.F., Herbut, I.F.
Phys. Rev. X 3, 031010, 2013
[6] Hohenadler, M., Assaad, F.F.
J. Phys.: Condens. Matter 25, 143201, 2013
[7] Uehlinger, T., Jotzu, G., Messer, M., Greif, D., Hofstetter, W., Bissbort, U., Esslinger, T.
arXiv:1308.4401
[8] 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
[9] Assaad, F.F.
Phys. Rev. Lett. 83, 796 (1999)
[10] Hohenadler, M., Lang, T.C., Assaad, F.F.
Phys. Rev. Lett. 106, 100403 (2011)
[11] Hohenadler, M., Meng, Z. Y., Lang, T. C., Wessel, S., Muramatsu, A., Assaad, F. F.
Phys. Rev. B 85, 115132, 2012
[12] Berg, E., Metlitski, M.A., Sachdev, S.
Science 338, 1606, 2012
[13] Rachel, S., Le Hur, K.
Phys. Rev. B 82, 075106, 2010
[14] Assaad, F.F., Bercx, M., Hohenadler, M.
Phys. Rev. X 3, 011015, 2013
[15] Herbut, I.F.
Phys. Rev. Lett. 97, 146401, 2006
[16] Herbut, I.F., Juričić V., Vafek, O.
Phys. Rev. B 80, 075432, 2009
[17] Yang, H.Y., Albuquerque, A.F., Capponi, S., Läuchli, A.M., Schmidt, K.P.
New Journal of Physics 14, 115027, 2012
[18] Chang, C.C., Scalettar, R.T.
Phys. Rev. Lett. 109, 026404, 2012
[19] Clark, B.K.
arXiv:1305.0278
[20] Sorella, S., Tosatti, E.
Europhys. Lett. 19, 699, 1992
[21] Paiva, T., Scalettar, R.T., Zheng, W., Singh, R.R.P., Oitmaa, J.
Phys. Rev. B 72, 085123, 2005
[22] Qi, X.L., Zhang, S.C.
Phys. Rev. Lett. 101, 086802, 2008
[23] Ran, Y., Vishwanath, A., Lee, D.H.
Phys. Rev. Lett. 101, 086801, 2008
• Fakher F. Assaad
Martin Hohenadler
Institut für
Theoretische Physik
und Astrophysik
Universität Würzburg
top

