Ground state phase diagram of the half-filled bilayer Hubbard model
Michael Golor, Timo Reckling, Laura Classen, Michael M. Scherer, Stefan Wessel
GGround state phase diagram of the half-filled bilayer Hubbard model
Michael Golor, ∗ Timo Reckling, Laura Classen, Michael M. Scherer, † and Stefan Wessel Institut f¨ur Theoretische Festk¨orperphysik, JARA-FIT and JARA-HPC, RWTH Aachen University, 52056 Aachen, Germany Institut f¨ur Theoretische Physik, Universit¨at Heidelberg, 69120 Heidelberg, Germany (Dated: October 2, 2018)Employing a combination of functional renormalization group calculations and projective deter-minantal quantum Monte Carlo simulations, we examine the Hubbard model on the square latticebilayer at half filling. From this combined analysis, we obtain a comprehensive account on theground state phase diagram with respect to the extent of the system’s metallic and (antiferromag-netically ordered) Mott-insulating as well as band-insulating regions. By means of an unbiasedfunctional renormalization group approach, we exhibit the antiferromagnetic Mott-insulating stateas the relevant instability of the free metallic state, induced by any weak finite onsite repulsion.Upon performing a careful analysis of the quantum Monte Carlo data, we resolve the difficulty ofidentifying this antiferromagnetic ground state for finite interlayer hopping in the weak-couplingregime, where nonmonotonous finite-size corrections are shown to relate to the two-sheeted Fermisurface structure of the metallic phase. On the other hand, quantum Monte Carlo simulations arewell suited to identify the transition between the Mott-insulating phase and the band insulator inthe intermediate-to-strong coupling regime. Here, we compare our numerical findings to indicationsfor the transition region obtained from the functional renormalization group procedure.
I. INTRODUCTION
The Hubbard model provides a fundamental descrip-tion of correlation effects in many-body condensed mat-ter systems. On the square lattice bilayer, it exhibitsseveral different routes for the transition out of normalmetallic behavior by an independent tunability of (i) thefree band-structure via a single tight-binding parame-ter, i.e., the interlayer hopping amplitude, (ii) the on-site repulsion, as well as (iii) the filling. For example,at half filling, the interlayer hopping controls a contin-uous transition from a Mott-insulating state to a bandinsulator . Furthermore, at finite doping the modelshows the appearance of superconducting pairing insta-bilities and by tuning the interlayer hopping, a transitionbetween distinct types of pairing symmetries can be ac-complished, cf. e.g. Ref. 12. These pairing symmetries,of d -wave and extended s -wave character, are frequentlydiscussed in the context of the cuprate and pnictide com-pounds, respectively . Finally, the bilayer Hubbardmodel was also discussed as a simple model with multi-ple Fermi surfaces and electron and hole pockets, cf. e.g.Ref. 18, as well as in the context of copper oxide bilay-ers featuring high- T c superconductivity . The bilayerHubbard model might as well become a rich and instruc-tive example system with a well-defined and at the sametime complex phase diagram, given the progress in thecontrolled manipulation of ultracold atoms in optical lat-tices .The bilayer Hubbard model at half filling has beenstudied independently using various methods recently,ranging from dynamical mean field theory (DMFT) andcluster extensions and variational Monte Carlo , tofinite-temperature determinantal quantum Monte Carlo(DQMC) studies . At present, however, it appearsthat some of the results on its phase diagram remaininconclusive – especially since the persistence of an ex- tended paramagnetic metallic phase at small onsite inter-actions has been put forward based on cluster DMFT aswell as finite temperature DQMC calculations – a result,which contrasts the perfect nesting property of the non-interacting system at half filling throughout the wholemetallic regime. The above outlined perspectives for thephysics of the bilayer Hubbard model hence call for aclarification concerning its ground state phase diagramat half filling.Here, we provide a comprehensive account on theground state phase diagram by combining complemen-tary methods, which – taken together – allow us tocover the full range of the local interaction strength. Forthis purpose, we first review the standard Hartree-Fockmean-field theory (HFMFT) approximation to properlysettle the mean-field character of the antiferromagneticstate as a function of the onsite interaction, as arisingfrom a Stoner instability of the metallic noninteractingstate. In a second step, which leads systematically be-yond mean-field theory as well as random phase approxi-mations (RPA), we implement a functional renormaliza-tion group (fRG) method in the discrete patching scheme,which provides a well-controlled approach at weak inter-actions and allows us to obtain reliable results for inter-mediate interaction strengths . An important bene-fit provided by this method is that it takes into accounteffects from possibly competing correlations and there-fore allows us to detect the appearance of instabilitiesin an unbiased way, i.e. without a priori assumptionsconcerning the nature of the emerging order. Finally,we employ an unbiased and numerically exact method,zero-temperature (projective) DQMC, which is particu-larly powerful in the regime of intermediate to strongcoupling and which allows us to identify actual groundstate correlations on finite systems. It turns out that es-pecially the mutual strengths of these approaches help toprovide a coherent picture of the low-energy physics of a r X i v : . [ c ond - m a t . s t r- e l ] N ov the bilayer Hubbard model.The rest of this paper is organized as follows: In Sec. II,we introduce the Hubbard model on the square lattice bi-layer and also summarize our central findings from thesubsequent sections by presenting the resulting groundstate phase diagram of the half-filled system. Sec. IIIprovides a detailed account on the weak-coupling regimefrom the perspective of three different many-body meth-ods: HFMFT, fRG, and DQMC. Our projective DQMCcalculations exhibit rather severe, nonmonotonous finite-size effects in the weak-coupling regime, that most likelyexplain conflicting previous conclusions on the persis-tence of the metallic state drawn from finite-temperatureDQMC simulations (which also require control of thesimulation temperature appropriately in order to sam-ple ground state correlations). In Sec. IV, we examinethe transition from the antiferromagnetic Mott insulatorto the strong-interlayer-hopping band insulator and con-trast its identification within the DQMC and the fRGapproach. Final conclusions will be presented in Sec. V. II. MODEL AND PHASE DIAGRAM
In the following, we examine the Hubbard model on thesquare lattice bilayer with intra-layer nearest-neighborhopping t , inter-layer hopping t ⊥ and a local Coulombrepulsion U , as described by the Hamiltonian H = − t (cid:88) (cid:104) ij (cid:105) sλ (cid:16) c † iλs c jλs + h . c . (cid:17) (1) − t ⊥ (cid:88) is (cid:16) c † i s c i s + h . c . (cid:17) + U (cid:88) iλ n iλ ↑ n iλ ↓ , where the c ( † ) iλs denote annihilation (creation) operatorsfor electrons of spin s ∈ {↑ , ↓} on site i in layer λ ∈ { , } ,and the n iλs denote local occupation number operators.At half-filling, our calculations reveal the phase dia-gram shown in Fig. 1. The metallic phase only persistsin the noninteracting case ( U = 0) for interlayer hoppings t ⊥ /t <
4. For larger values of t ⊥ /t , the system becomesinsulating due to the opening of a band gap. In the weakcoupling regime ( U (cid:46) t ), the system is an antiferro-magnetic insulator state, triggered by a Stoner instabil-ity, which is due to the perfect nesting property of thetwo-sheeted Fermi surface at wave vector Q = ( π/a, π/a )with the lattice constant set to a = 1 in the following, seeFig. 2. The phase transition to the band insulator stilloccurs at t ⊥ ≈ t . For larger couplings U/t , this phaseboundary is continuously shifted to t ⊥ = 1 . t in the U → ∞ limit: in this regime, the low-energy spin physicscan be mapped onto that of a spin-1 / .The following sections provide details of our analyticaland numerical findings on which the identification of thephase diagram is based. III. WEAK-COUPLING REGIMEA. Tight-Binding and Hartree-Fock Approximation
Since the physics in the weak-coupling regime isstrongly influenced by the structure of the Fermi surfaceof the noninteracting system, we start our analysis bybriefly reviewing this case. For U = 0, Eq. (1) reduces toa tight-binding model, which can be solved by exact di-agonalization and exhibits the single-particle dispersion (cid:15) ± ( k , t ⊥ ) = − t (cos( k x ) + cos( k y )) ± t ⊥ . (2)This corresponds to two copies of the square latticesingle-layer dispersion ( t ⊥ = 0), shifted with respect toeach other linearly in t ⊥ . Since for the single-layer latticethe energies lie within the range − t ≤ (cid:15) ≤ t , a bandgap opens for t ⊥ > t and the system becomes band-insulating at half filling. The Fermi surface for variousvalues of t ⊥ /t is shown in Fig. 2. Both bands remainperfectly nested for all values t ⊥ /t , with a nesting vector Q = ( π, π ), i.e. (cid:15) +0 ( k + Q , t ⊥ ) = − (cid:15) − ( k , t ⊥ ) . (3)As a consequence of Eq. (2), the density of states(DOS) of the bilayer lattice, ρ ( (cid:15), t ⊥ ), is the sum of twodisplaced single-layer DOS contributions, see Fig. 3. Thiscauses the logarithmic van Hove singularity of the single-layer square lattice to be shifted away from the Fermilevel to energies ± t ⊥ . As a consequence, one expects asuppression of interaction effects in the bilayer system. t ⊥ / t U/t tt ⊥ fRGQMC0123456 ∞ AntiferromagneticInsulatorBand InsulatorMetal
FIG. 1. (color online) Phase diagram of the Hubbard modelon the square lattice bilayer (inset) at half filling. Red dotsdisplay the phase boundary obtained from projective DQMC(the red solid line is a guide to the eye). The blue dashedline indicates the reduction of the sharply pronounced an-tiferromagnetic momentum structure upon increasing t ⊥ atintermediate values of U/t in the fRG approach, hinting to-wards a phase boundary of the antiferromagnetic instability(cf. Sec. IV B for details). − π − π π π − π − π π π k y k x Q t ⊥ = 0 − π − π π πk x Q t ⊥ = t − π − π π πk x Q t ⊥ = 2 t − π − π π πk x Q t ⊥ = 3 t FIG. 2. (color online) Fermi surfaces of the square lattice bilayer tight-binding model ( U = 0) at half filling for different valuesof t ⊥ /t = 0 , , , Q = ( π, π ).The grid points denote the momenta included in a finite system of linear length L = 6 with periodic boundary conditions. Forthe system sizes accessible in the DQMC simulations, these represent the only values of t ⊥ /t with discrete momenta locatedon the Fermi surface (larger dots). However, for t ⊥ < t , the DOS stays finite at the Fermisurface. This, combined with the nesting property, leadsto a divergence of the zero-temperature static spin-spinsusceptibility χ + − ( q ) = − N (cid:88) k f (cid:0) (cid:15) +0 ( k + q , t ⊥ ) (cid:1) − f (cid:0) (cid:15) − ( k , t ⊥ ) (cid:1) (cid:15) +0 ( k + q , t ⊥ ) − (cid:15) − ( k , t ⊥ ) , (4)where f ( (cid:15) ) denotes the Fermi-Dirac distribution. At thenesting vector, χ + − ( Q ) → ∞ (for temperature T → n iλ ↑ n iλ ↓ → n iλ ↑ (cid:104) n iλ ↓ (cid:105) + (cid:104) n iλ ↑ (cid:105) n iλ ↓ − (cid:104) n iλ ↑ (cid:105)(cid:104) n iλ ↓ (cid:105) (5)and expressing the occupation expectation values interms of a homogeneous staggered magnetization m = |(cid:104) n iλ ↑ (cid:105) − (cid:104) n iλ ↓ (cid:105)| /
2. The resulting effective Hamiltoniancan again be solved analytically (see App. A for details)and the dispersion is now given by four bands labeled byan index α ∈ { , ..., } , (cid:15) α HF ( k , t ⊥ ) = ± (cid:113)(cid:0) (cid:15) ± ( k , t ⊥ ) (cid:1) + ∆ , ∆ ∝ t e − C t ⊥ t/U , (6)with a positive coefficient C t ⊥ that increases with t ⊥ /t .Here, ∆ = U m corresponds to the energy gap at theFermi surface and decays exponentially with decreasingcoupling strength
U/t . Nonetheless, the gap and thestaggered magnetization stay finite for any finite interac-tion, indicating the onset of antiferromagnetism for theweakly interacting system. Previous finite-temperatureDQMC simulations concluded in favor of a finite extentof the paramagnetic, metallic region for finite t ⊥ , arguingthat the tendency of the interlayer coupling towards theformation of spin-singlet states, which indeed prevails in . . . . . . .
30 -6 -4 -2 0 2 4 6 ρ (cid:15)/t t ⊥ = 0 t ⊥ = 2 t FIG. 3. (color online) The U = 0 bilayer density of states(DOS) ρ ( (cid:15), t ⊥ ) for t ⊥ = 2 t (solid red line) in comparison withthe single-layer DOS (dashed black line). The dotted verticalline marks the position of the Fermi level. the large- U region, might compete with the HFMFT an-tiferromagnetic instability also at low values of U/t . Aswill be shown below, this conclusion reflects an inherentdifficulty of performing a systematic finite-size analysisof DQMC data in the weak-coupling region, due to non-monotonous finite-size effects, even when ground stateexpectation values are accessed. In the following section,we first employ an fRG approach that shows in an un-biased way how the antiferromagnetic instability indeedemerges in the weak-interaction regime from the nonin-teracting metallic state.
B. Functional Renormalization Group
Here, we employ the fRG method for the one-particle-irreducible (1PI) vertices of a fermionic many-body sys-tem to obtain their renormalization group (RG) evolutionfor the bilayer system, see e.g. Refs. 23 and 25 for detailson the general fRG approach. To that end, the bare prop-agator of the action corresponding to the model definedin Eq. (1) is modified by introducing an infrared regula-tor with energy scale Λ. Varying Λ generates the RG flowin terms of a functional differential equation. The flow isthen integrated out starting with the bare action at theinitial scale Λ corresponding to the bandwidth of thesystem. We integrate towards the infrared limit Λ → V Λ ( k , k ; k , k ) where k i con-sists of a wavevector k i , a Matsubara frequency ω i , aspin projection s i , and band indices b i or layer indices λ i .External frequencies are set to zero to study the groundstate and furthermore, the momentum dependence is dis-cretized. We also neglect selfenergy corrections to reducethe numerical effort. The vertex V Λ is discretized by di-viding the Brillouin zone into N p patches with a constantwavevector dependence within each patch, as shown inFig. 4. Representative momenta for these patches arechosen to reside at the Fermi level. In the following, weemploy a wavevector resolution with a patch number of N p = 32 and N p = 48. We have also checked our re-sults to be stable towards higher resolutions with up to N p = 96 patches. In summary, we obtain a vertex func-tion V Λ with N p · N b components, where N b = 2 is thenumber of energy bands, and a set of N p · N b coupleddifferential equations that has to be integrated.The approximations involved in this approach cor-respond to an infinite-order summation of one-loopparticle-particle and particle-hole diagrams, cf. Fig. 4,which allows the study of various competing correlationeffects. Eventually, this leads to instabilities in the RGflow at a critical scale Λ c , where components of the vertexgrow large. Such a divergence of the interaction vertexcan be seen as an artifact of our truncation due to theapproximations described above. On the other hand adivergence in the interaction vertex is a dependable in-dicator for a tendency towards the formation of an or-dered state and the well-pronounced emerging momen-tum structure close to the critical scale allows us to iden-tify an effective Hamiltonian for the low-energy regime.In practice, the flow is halted at a finite value for thelargest vertex component of the order of several timesthe bandwidth. This defines an energy scale which givesa reasonable approximation to the critical scale as thevertex diverges quickly close to the instability . Nearthe scale Λ c , the dominant correlations become clearlydetectable in the vertex function and we can use thisscale as an estimate for e.g. the energy scale below whichthe single-particle spectrum gets modified, e.g. by a gap. (cid:1)(cid:2)(cid:3) (cid:1)(cid:4)(cid:3)(cid:1)(cid:5)(cid:3) (cid:1)(cid:2)(cid:3)(cid:4) (cid:5) (cid:1)(cid:2)(cid:3)(cid:4) (cid:6) (cid:1)(cid:7)(cid:2)(cid:3)(cid:4) (cid:8) (cid:1)(cid:7)(cid:2)(cid:3)(cid:4) (cid:9) (cid:45) (cid:45) (cid:45) (cid:45) k x (cid:144) Π k y (cid:144) Π FIG. 4. (color online) Left panel on top: Interaction vertexand spin convention. Below: Loop contributions from theparticle-particle channel (a), the crossed particle-hole chan-nel (b) and the direct particle-hole channel (c). Right panel:Patching scheme of the Brillouin zone for a total of N p = 32patches. Depending on the energy band the coupling functionis evaluated on a wave vector at the Fermi level indicated bythe dots (red: upper band, blue: lower band). The procedure described here is well-controlled for smallinteractions and can be expected to be reliable also inthe regime of intermediate interaction strengths .Importantly, the fRG in this approximation scheme takesinto account effects from competing correlations and al-lows to identify the leading instabilities in an unbiasedway, i.e., no a priori assumptions concerning the natureof the appearing order are required. The effect of the in-clusion of self-energy feedback and frequency dependentinteractions has been studied in the single-layer Hubbardmodel in Ref. 32 where it has been shown that the lead-ing channel is not strongly affected by this extension ofthe truncation.In the fRG data at half filling and for arbitrary on-site repulsion, we observe as the prevailing divergencean antiferromagnetic spin density wave (AF-SDW) withmomentum transfer Q = ( π, π ), see Fig. 5 for a snapshotof the four-Fermi vertex close to the critical scale. Thisbehavior clearly reflects the perfect nesting of the Fermisurface. The leading part of the interaction close to thecritical scale can be expressed in terms of an effectiveinteraction Hamiltonian H SDW = − J (cid:88) λ,i,j e i Q · ( r i − r j ) (cid:0) S iλ · S jλ − S iλ · S j ¯ λ (cid:1) , (7)where J >
0, and ¯ λ denotes the layer opposite to λ .The spin operator S iλ is given by the relation S iλ =1 / (cid:80) s,s (cid:48) c † iλs σ ss (cid:48) c iλs (cid:48) in terms of the Pauli matrices. Asa result of the sharpness in momentum space the inter-action becomes long-ranged in position space. With thiseffective Hamiltonian, we can thus perform a controlledmean-field decoupling. For the resulting effective Hamil-tonian, the spins are aligned antiferromagnetically withineach layer and also between the layers.In Fig. 6, we compare the critical scales for severalcases, namely the situation with vanishing perpendicu-lar hopping amplitude t ⊥ = 0 which is equivalent to FIG. 5. (color online) Effective interaction vertex near thecritical scale for the AF-SDW in units of t , exhibiting asharply pronounced momentum structure. The axes are num-bered according to the number of the patch, cf. Fig. 4,where wavevectors k are depicted vertically and k are de-picted horizontally. We fix k to be on patch 1. Left Panel:Effective vertex or a combination of layer indices λ i where λ = λ = λ = λ . Middle Panel: λ = λ , λ = λ . Rightpanel: λ = λ , λ = λ . t ¶ ê t = t ¶ ê t = t ¶ ê t = t ¶ ê t = t ê U L c ê t U ê t = U ê t = U ê t = t ¶ ê t FIG. 6. (color online) Left panel: fRG critical scale Λ c asa function of t/U for the single-layer case (solid line) andthe bilayer case with t ⊥ = 1 t, t, t (from top to bottom)using N p = 48 patches, exhibiting exponential decrease ofΛ c ∼ e − const . · t/U . This functional dependence suggests aninstability for any value of U > c decreases by severalorders of magnitude when t ⊥ /t is increased from 0 to 4. Weshow this behavior for three choices of U/t = 3 , U/t = 2 and
U/t = 1 by the solid, dashed and dotted lines, respectively. the one-band Hubbard model as well as for perpendic-ular hoppings of t ⊥ /t = 1 , , U/t is an ex-ponential ∼ e − const . · √ t/U , cf. Fig. 6. In contrast, thecritical scales for the bilayer case are considerably lower,which can be understood given the reduced density ofstates at the Fermi level. Furthermore, the functionaldependence on U/t is different and rather follows a be-havior ∼ e − const . · t/U which is in accordance with theHFMFT expectation (cf. Sec. III A). Using the fRG ap- proach, treating all the fluctuation channels on equalfooting, we can also check whether competing instabil-ities are present for this choice of parameters. Here,we do not observe the appearance of any other strongcorrelation effects apart from the AF-SDW instability.This also constitutes a worthwhile cross-check for theDQMC approach, excluding competing correlations aspossible contributions to the observed finite-size scaling(cf. Sec. III C).In the right panel of Fig. 6 we compare the t ⊥ -dependence of the critical scale Λ c for different values of U/t ∈ { , , } . We observe a continuous decrease of Λ c as a function of increasing interlayer hopping, in accordwith the HFMFT approximation, cf. App. A. C. Quantum Monte Carlo
We next employ projective ( T = 0) DQMC simula-tions to study the ground state correlations on finite sys-tems. The DQMC approach enables us to obtain groundstate expectation values of an arbitrary observable O byprojecting a trial wave function | Ψ T (cid:105) to the interactingground state, (cid:104) Ψ | O | Ψ (cid:105)(cid:104) Ψ | Ψ (cid:105) = lim Θ →∞ (cid:104) Ψ T | e − Θ H O e − Θ H | Ψ T (cid:105)(cid:104) Ψ T | e − H | Ψ T (cid:105) . (8)Here, we take the ground state of the free system ( U = 0)as the trial wave function in all simulations. Possibledegeneracies are lifted by a slight dimerization of theintralayer hopping amplitude t of the tight-binding sys-tem ( δt/t = 10 − ). The projection parameter Θ hasto be chosen sufficiently large, such that convergence tothe ground state | Ψ (cid:105) is guaranteed. Depending mostlyon the linear system size L and the coupling strength U/t , values of Θ between Θ = 30 /t and 120 /t are re-quired to ensure convergence. We employ a symmet-ric Suzuki-Trotter decomposition with an imaginary-timediscretization of ∆ τ = 0 . /t , such that discretization er-rors are well below the size of the statistical errors. TheHubbard interaction is decoupled by a SU (2) symmet-ric Hubbard-Stratonovich transformation, so that spinrotational symmetry is preserved during the entire sim-ulation. For a detailed account of the projective DQMCalgorithm, see Ref. 34. The available computational re-sources allow us to simulate systems with linear lengthsof up to L = 24, i.e. lattices containing up to N = 2 · L =1152 sites. In all cases, periodic boundary conditions inboth dimensions are applied, ensuring translational sym-metry.In the weak-coupling regime the physics will be domi-nated by the structure of the Fermi surface and its nest-ing properties in particular. The finite systems with peri-odic boundary conditions, treated within DQMC, trans-late in momentum space to a finite set of k -points sam-pling the Brillouin zone. In order to include the relevantlow-energy processes, we need to ensure that this set con-tains points at the Fermi surface. Due to the specific − − − − (0 ,
0) ( L ,
0) ( L , L ) (0 , | (cid:104) S r λ · S r (cid:48) λ (cid:105) | r − r (cid:48) U = tU = 2 tU = 3 tU = 4 tL = 24 , t ⊥ = 2 t FIG. 7. (color online) Static spin-spin correlations along atriangular path (inset) for different coupling strengths
U/t .The calculations were performed for a bilayer system withinterlayer hopping t ⊥ = 2 t and linear length L = 24. Fermi surface structure of the bilayer square lattice, thisconsiderably restricts our choice of t ⊥ and L , see Fig. 2.In particular, only for interlayer hopping t ⊥ /t ∈ { , , } ,there are several system lengths L ≤
24 available such asto allow for at least a limited finite-size analysis.The DQMC method allows for a direct access to staticspin-spin correlations (cid:104) S r λ · S r (cid:48) λ (cid:48) (cid:105) between two sites atpositions r , r (cid:48) in layers λ, λ (cid:48) . For the bilayer lattice with t ⊥ = 2 t and L = 24, the spin-spin correlations along a tri-angular path within one layer of the lattice are shown inFig. 7. With increasing distance, the correlations decayquickly into the long-distance regime. The magnitudeof the asymptotic value falls off drastically upon reduc-ing U and for the lower values of U , the curves exhibitcharacteristic anomalies, which already hint towards thedifficulty of extrapolating the finite-size data to the ther-modynamic limit that will be considered in more detail,below.A spin density wave (SDW) state with ordering vector q can be identified by monitoring the structure factor S ± ( q ), which is obtained by a Fourier transform of themeasured spin-spin correlations, S ± ( q ) = 1 N (cid:88) r , r (cid:48) ,λ e i q · ( r − r (cid:48) ) [ (cid:104) S r λ · S r (cid:48) λ (cid:105) ± (cid:104) S r λ · S r (cid:48) ¯ λ (cid:105) ] , (9)where ¯ λ denotes the layer opposite to λ . Due of theperfect nesting and the bipartite nature of the bilayerlattice, we expect an antiferromagnetic (AF) ordering,i.e. a SDW with ordering vector Q = ( π, π ) and the corre-sponding antisymmetric structure factor S AF = S − ( Q ).The order parameter for the AF state, the staggeredmagnetization m , is then usually estimated by ¯ m = (cid:112) S AF /N . However, this estimator is disadvantageous in situations where the long distance correlations are of verylow magnitude, as is the case in Fig. 7 for small valuesof U/t . The AF structure factor will then be dominatedby local correlations, which are not suited to character-ize a long-range ordered antiferromagnetic state on finitelattices. We therefore restrict the summation in Eq. (9)to sites r , r (cid:48) of distance | r − r (cid:48) | > L/
4, where in allconsidered cases the correlations have decayed to theirasymptotic level. The prefactor 1 /N is modified accord-ingly. We denote this structure factor by S AF L/ and thecorresponding estimator for the staggered magnetizationby ¯ m L/ .The results for ¯ m L/ are shown in Fig. 8 as a functionof the inverse system size 1 /L for interlayer hoppings t ⊥ = t and 2 t . At this point a finite-size scaling analy-sis needs to be performed, e.g. in terms of a polynomialextrapolation to the thermodynamic limit at 1 /L = 0to access a thermodynamic limit estimate for the stag-gered magnetization. This procedure works reasonablywell at couplings U (cid:38) t , for which a linear extrapola-tion yields robust finite estimates for the staggered mag-netization, thus indicating an emerging AF state. How-ever, for weaker interactions, the finite-system estimatesfor the staggered magnetization exhibits a strong non-monotonous behavior upon varying the system size, thatdefies a similarly reliable extrapolation.This peculiar finite-size behavior can be traced back tothe coarse sampling of the Fermi surface for small finitesystems. In fact, a similar nonmonotonous dependenceon the system size can be observed already in the nonin-teracting system’s static spin-spin susceptibility, Eq. (4):While at T = 0 the susceptibility diverges at the nestingvector Q , at low finite temperatures it exhibits a verysimilar finite-size dependence as the one observed for theAF structure factor at small finite values of U . This isevident from Fig. 9, where both quantities are comparedfor the case of t ⊥ = 2 t . From the behavior of χ + − ( Q )it appears that if only system sizes with L being mul-tiples of 12 are considered, a quadratic behavior can bewell fitted to the finite-size susceptibility data, yieldinga finite estimate in the thermodynamic limit. In orderto perform a similar meaningful finite size scaling for thefinite- U structure factor data, results for system lengths L = 36 , , . . . would be required, which are currently outof reach. On the other hand, the low- U data for ¯ m L/ in Fig. 8 may also be argued to exhibit a leveling-off be-havior towards a small, finite value at the lower rangeof 1 /L values. While thus no definite statement can bemade based on our data about the size of the staggeredmagnetization at the lower values of U in the thermo-dynamic limit, the data in Fig. 8 is consistent with thepersistence of a small, but finite value of the staggeredmagnetization in the thermodynamic limit for both U = t and U = 2 t , which would be in accord with the weak cou-pling scenario from the fRG analysis. Previous DQMCsimulations suggested the existence of a paramagneticmetallic phase at small values of U/t ; however, the cor-responding calculations were performed at low but finite . . . . . . . .
00 0 .
02 0 .
04 0 .
06 0 .
08 0 . ¯ m L / /Lt ⊥ = t U = tU = 2 tU = 3 tU = 4 tU = 5 t . . . . . . .
00 0 .
02 0 .
04 0 .
06 0 .
08 0 .
10 0 .
12 0 . ¯ m L / /Lt ⊥ = 2 t U = tU = 2 tU = 3 tU = 4 tU = 5 t FIG. 8. (color online) Finite size scaling of the staggeredmagnetization for t ⊥ = t (top) and t ⊥ = 2 t (bottom). Forcouplings U/t ≥
4, linear fits were performed (solid lines). temperatures and were restricted to linear system sizes L ≤
10 only. Based on our ground state analysis with L up to 24, we cannot confirm the persistence of the param-agnetic phase. Our DQMC data are also consistent withthe (more conventional) scenario supported by the fRGcalculations, restricting the range of the paramagneticmetallic phase down to the U = 0 line.Within this scenario, the AF Mott insulating phase isaccompanied by a finite single-particle gap. Based on theDQMC calculations, we can extract this excitation gapfrom the large-imaginary time limit of the MatsubaraGreen’s function, G k ( τ ) = −(cid:104) c k ( τ ) c † k (0) (cid:105) τ →∞ ∝ e − ∆( k ) τ . (10)The single particle gap ∆ sp is then defined as the small-est occurring value, ∆ sp = min { ∆( k ) } . This quantity al-lows for a direct comparison to the gap obtained withinHFMFT as well as the critical scale obtained from the .
00 0 .
02 0 .
04 0 .
06 0 .
08 0 . S A F L / × χ + − ( Q ) t /L t ⊥ = 2 t S AF L/ /N, U = tS AF L/ /N, U = 2 tχ + − ( Q ) FIG. 9. (color online) Comparison of the antiferromagneticstructure factor of the weakly interacting system and thenoninteracting static susceptibility χ + − ( Q ) at finite temper-ature T = 10 − t (blue dots) for different linear lengths L . Aquadratic function may be fitted to the susceptibility data forlengths L that are multiples of 12 (blue solid line). fRG calculations. In Fig. 10, the corresponding resultsfrom the three different methods are compiled for thesingle layer system (left panel) and the bilayer systemat t ⊥ = 2 t (right panel). Generally, QMC and fRG re-sults indicate the exponential behavior obtained withinthe HFMFT calculations. The deviations that are ob-served at weak coupling in the finite-size QMC data arealso obtained within HFMFT calculations on such finitesystems. This indicates again the presence of sizablefinite-size effect in the low- U regime. For the interact-ing single-layer system, the immediate onset of antiferro-magnetism at any finite U is by now well established .Even though a reliable extrapolation of the QMC data islimited to couplings U ≥ t for the Hubbard bilayer, thecomparison in Fig. 10 supports the validity of the fRGscenario in the weak coupling regime and thus makes afurther case for the emergence of the AF Mott insulatorphase for any weak, finite interaction strength from thefree metallic region. IV. ANTIFERROMAGNET TO BANDINSULATOR TRANSITION
As already discussed in Sec. III A, the free system( U = 0) undergoes a transition from the metal to a bandinsulator at t ⊥ = 4 t . On the other end of the interactionrange, i.e. in the strong coupling limit U → ∞ , the low-energy spin dynamics in the bilayer Hubbard model canbe effectively described by a spin-1 / H = J (cid:88) (cid:104) ij (cid:105) λ S iλ · S jλ + J ⊥ (cid:88) i S i · S i , (11)with antiferromagnetic couplings J = 4 t /U and J ⊥ =4 t ⊥ /U . For this spin model, a quantum phase transi- − . . . . . . ∆ s p / t , Λ c / t t/Ut ⊥ = 0 0 . . . . t/Ut ⊥ = 2 t MFT, L = 24MFT, TDLQMC, L = 24QMC, TDLfRG FIG. 10. (color online) Comparison of the single-particle gapsfrom HFMFT (solid and dashed lines), QMC (red dots andblack squares) and the critical scale from fRG (blue triangles).Left: Single-layer case. Right: Bilayer system with t ⊥ = 2 t . tion from an antiferromagnet low- J ⊥ phase to a mag-netically disordered dimerized phase with spin singletspredominantly formed on the interlayer bonds occurs at J ⊥ /J = 2 . . With respect to the Hubbardmodel parameters, this corresponds to a critical inter-layer hopping of t ⊥ /t = 1 . U = 0 band insulator phase, while theAF Mott insulator is separated from the large- t ⊥ bandinsulator regime by a magnetic quantum phase transitionrelated to the breaking of the spin rotation symmetry inthe AF phase. In this section, we discuss how the extentof the AF phase can be extracted from DQMC and fRGfor finite values of the interaction strength U/t . A. Quantum Monte Carlo
We employ finite size scaling of the DQMC structurefactor data to locate the transition between the Mott-insulating phase and the band-insulating regime. In thefollowing, we locate the transition points upon varying U at fixed values of t ⊥ /t . Given the limits to be t ⊥ = 4 t for U = 0 and t ⊥ = 1 . t for U → ∞ , we extract the criticalcoupling U c within the intermediate interlayer hoppingrange at t ⊥ /t = 2 . , .
5, and 3 .
0. In order to determine U c , we make use of finite size scaling with the followingstandard scaling ansatz for the structure factor in theproximity of the quantum phase transition: S AF /N = L − βν F S (cid:0) uL ν (cid:1) , u = U − U c U c , (12)where F S is the scaling function of the structure fac-tor, and β , ν denote the critical exponents for the or- der parameter and the correlation length, respectively.At the critical coupling u = 0, the scaling function willbe evaluated at the same point for different L , so thatthe rescaled structure factors L βν S AF /N should cross at U = U c . Thus, by monitoring this observable over arange of interactions U/t , curves for different L shouldintersect at the critical coupling U c /t . Since the onsetof AF order in the ground state of this two-dimensionalquantum system breaks the SU (2) spin rotational sym-metry, we anticipate that the quantum phase transitionbelongs to the universality class of the three-dimensionalclassical Heisenberg model . We thus employ the crit-ical exponents of this universality class β = 0 . ν = 0 . . From the corresponding analysisshown in Fig. 11, we extract for the critical line U c ( t ⊥ /t )the three points U c (3 .
0) = 5 . t , U c (2 .
5) = 6 . t and U c (2 .
0) = 10 . t . The estimated phase boundarypoints have already been indicated in Fig. 1. The datacollapses shown in Fig. 12 demonstrate, that the finite-size data indeed follows the anticipated scaling ansatzin Eq. (12) within the vicinity of the quantum criticalpoints. Fig. 12 also exhibits an increase of further finite-size corrections for lower values of the critical interactionstrength, in accord with our findings in Sec. III C.At this point, it is interesting to observe, that the en-hanced susceptibility towards the free system’s Fermi sur-face structure at low values of U in effect leads to a changein the finite-size scaling behavior of the structure factor S AF . This is seen from e.g. the right panel of Fig. 13at t ⊥ = 2 t , which could be mistaken to indicate a mag-netic phase transition near U/t ≈
3. These apparentcrossings instead signal a crossover in the characteristicsof the AF state: while at low values of U , the AF orderarises from a Fermi surface instability, it relates for largervalues of U to the localized nature of the Mott insulatorwith the emerging Heisenberg spin physics. In fact, asimilarly misleading indication for such a transition isobserved in the decoupled layer ( t ⊥ = 0) data, shown inthe left panel of the same figure, near U/t ≈ .
5. Fromour DQMC simulations, we thus find that on the bilayer,this crossover behavior is even more pronounced than forthe single-layer model, due to the more complex Fermisurface structure in the bilayer case and the correspond-ing suppression of the van Hove singularity in the DOSaway from the half-filled system.
B. Functional Renormalization Group
It is interesting to assess, if and how the transitionout of the AF state upon increasing U can be identi-fied within the fRG approach. The fRG method is acontrolled approximation for weak interactions and hasbeen shown to provide reliable results also in the inter-mediate interaction regime. Here, we thus study the be-havior of the bilayer system over a range of U/t valuesfrom small interactions up to the order of the bandwidth.For small interlayer hopping t ⊥ /t and at any investigated . . . . . L β ν S A F L / / N U/tt ⊥ = 3 . t . . . . . . U/tt ⊥ = 2 . t . . . . U/tt ⊥ = 2 . tL = 12 L = 18 L = 24 L = 14 L = 18 L = 24 L = 12 L = 16 L = 24 FIG. 11. (color online) Finite size scaling analysis of therescaled antiferromagnetic structure factor around the AF-BItransition, assuming β = 0 .
369 and ν = 0 . U c /t of the transition. -4 -2 0 2 4 L β ν S A F L / / N uL ν t ⊥ = 3 . tU c = 5 . t -4 -2 0 2 4 uL ν t ⊥ = 2 . tU c = 6 . t -4 -2 0 2 4 uL ν t ⊥ = 2 . tU c = 10 . tL = 12 L = 18 L = 24 L = 14 L = 18 L = 24 L = 12 L = 16 L = 24 FIG. 12. (color online) Data collapse analysis of the rescaledantiferromagnetic structure factor around the AF-BI transi-tion, assuming β = 0 .
369 and ν = 0 . value of the onsite interaction U , we observe the appear-ance of the characteristic feature of a clear AF instabilityin the interaction vertex, which is sharply pronouncedin momentum space, as shown in Fig. 5. As explainedin Sec. III B, this specific vertex structure leads to theidentification of the AF insulator state. Interestingly, forlarger values of t ⊥ /t , we observe, as shown in Fig. 14,that these sharp structures in the vertex V Λ smear outas we increase U (cid:38) t . This smearing happens gradu-ally, and for fixed onsite interaction beyond a value of t ⊥ smaller than 4 t the characteristic AF structure van-ishes completely from the vertex function, cf. Fig. 14.The sizeable values of the onsite interactions investigatedhere, are usually expected to lie beyond the regime of ap-plicability of the presented fRG approach, operating onthis level of truncation. However, in light of the above re-sults from the DQMC simulations, we may interpret theabove behavior also as a signature for the breakdown ofthe AF order and use this to obtain an fRG estimate forthe transition region towards the band insulating state.Therefore, we resort to an effective parameterizationof the two-particle interaction vertex V Λ by using an ex- . . . . L β ν S A F L / / N U/tt ⊥ = 0 2 . . . . . U/tt ⊥ = 2 tL = 12 L = 16 L = 24 L = 12 L = 16 L = 24 FIG. 13. (color online) Apparent crossings of the rescaledstructure factor data for t ⊥ = 0 (left panel) and t ⊥ = 2 t (rightpanel) in the lower- U region, indicative of the crossover in thenature of the AF state. V (cid:76) (cid:144) t FIG. 14. (color online) left panel: Effective low-energy in-teraction vertex in units of t for a combination of layer in-dices λ i where λ = λ = λ = λ for three choices of t ⊥ ∈ { t, . t, . t } (from top to bottom) at fixed onsite in-teraction U = 6 t . We present the vertex at an RG scale Λ ∗ where its largest component has grown large, V Λ , max ∼ t .For these combinations of parameters we find Λ ∗ ∼ O (0 . t ).The conventions are identical to the ones in Fig. 5. Rightpanel: Interaction vertex profile for fixed k at patch 1 forthe t ⊥ from the left panel. This corresponds to a horizontalcut of the interaction vertices shown in the left panel. Theflat line in the bottom (green) shows the value of the initialinteraction. change propagator ∝ / ( C Λ q + m Λ ) following the idea ofa gradient expansion around ordering momenta . Theinverse of the maximum interaction component deter-mines m Λ = 1 /V Λ , max and the inverse width C Λ is relatedto the spin stiffness. From the flow of C Λ , we then readoff the AFMI to BI transition, see Fig. 15: In the AF0 (cid:76) C (cid:76) t (cid:166) (cid:144) t (cid:61) t (cid:166) (cid:144) t (cid:61) t (cid:166) (cid:144) t (cid:61) FIG. 15. (color online) Flow of the inverse width C Λ for threechoices of t ⊥ ∈ { t, . t, . t } (solid, dashed and dotted line,respectively) at fixed onsite interaction U = 6 t . regime C Λ grows large corresponding to a sharp momen-tum structure or a long-ranged correlation in positionspace. In contrast, the flow of C Λ stays finite in theband insulating phase, where the correlations becomeshort-ranged leading to dimers accross the two layers.The corresponding transition line is presented in Fig. 1.The above parameterization is reasonable in the regimewhere antiferromagnetic correlations become importantand otherwise is less reliable. In this way, transition lineshould only be taken as a rough estimate. Nevertheless,from our fRG data, we may thus ascertain the disap-pearance of the AF instability for sizable values of U and t ⊥ as expected for the Mott-insulator to band-insulatortransition and we may determine the approximate posi-tion of the transition region from the above procedure.Given the expected less quantitative accuracy of the fRGapproach in this coupling range, we consider the transi-tion as extracted above from the fRG to agree reasonablywell to the DQMC transition points, cf. Fig. 1.The above analysis reveals that the fRG approach al-lows to detect a softening of the AF order such as near thetransition to the band insulator phase. Hence, we takethe absence of any similar weakening in the AF structureof the interaction vertex in the low- U region as anotherargument against the relevance of interlayer singlet for-mation as a mechanism to stabilize a low- U paramagneticmetallic state. V. CONCLUSIONS
Based on a combined analysis using fRG and DQMCcalculations, we established the nature of the quantumphase diagram of the half-filled Hubbard model on thesquare lattice bilayer. In particular, we identified thedominant AF instability in the weak-coupling region asa Fermi surface instability of the metallic free state.Within both the fRG and the DQMC simulations, we didnot detect any competing instabilities at weak-coupling.Based on a careful analysis of the ground state DQMCdata for the staggered structure factor, we identified the main difficulty in performing a standard finite-size ex-trapolation in the weak-coupling region, due to an en-hanced susceptibility of the magnetic state on the finitesystem size, which relates to the complex Fermi surfacestructure in the free metallic region. Our DQMC datais in accord with the HFMFT and our fRG scenario inwhich the metallic region is restricted to the U = 0 line,reflecting the persistent perfect nesting in the bilayer lat-tice. Furthermore, we employed standard finite size scal-ing to trace out within DQMC the transition line betweenthe Mott insulator and the band-insulating regime, andalso obtained signatures for this transition from a care-ful analysis of the fRG interaction vertex flow. For thefuture, it will be interesting to address the stability ofthe AF phase with respect to extended interactions bothwithin the planes as well as upon the addition of inter-layer interactions. ACKNOWLEDGMENTS
We acknowledge discussions with Fakher F. Assaad,Stephan Haas, Carsten Honerkamp, Stefan A. Maier, JanM. Pawlowski, and Daniel D. Scherer. L.C. acknowledgessupport by the HGSFP. M.M.S. is supported by the grantERC-AdG-290623 and DFG grant FOR 723. M.G. andS.W. acknowledge support from the DFG within grantFOR 1807 and grant WE 3949/3-1. Furthermore, weacknowledge the allocation of CPU time through JARA-HPC at JSC J¨ulich and at RWTH Aachen.
Appendix A: Bilayer in Hartree-Fock approximation
We show in this appendix that the half-filled bilayersystem in HFMFT approximation is antiferromagneti-cally ordered for any finite value of
U/t for t ⊥ < t .In HFMFT, the interaction term is decoupled as n iλ ↑ n iλ ↓ → n iλ ↑ (cid:104) n iλ ↓ (cid:105) + (cid:104) n iλ ↑ (cid:105) n iλ ↓ − (cid:104) n iλ ↑ (cid:105)(cid:104) n iλ ↓ (cid:105) (A1)We now assume a staggered magnetization m on the sub-lattices σ = A, B . We divide the lattice in N C = N unitcells with four sites, labeled by ( λ, σ ), and make the fol-lowing ansatz for the mean field parameters of the AFstate: (cid:104) n A ↑ (cid:105) = + m (cid:104) n A ↓ (cid:105) = − m (A2) (cid:104) n B ↑ (cid:105) = − m (cid:104) n B ↓ (cid:105) = + m (A3) (cid:104) n A ↑ (cid:105) = + m (cid:104) n A ↓ (cid:105) = − m (cid:104) n B ↑ (cid:105) = − m (cid:104) n B ↓ (cid:105) = + m This, after Fourier transformation via c iλσs = (cid:112) /N C (cid:80) k e − i kr i c k λσs yields the following HFMFT1Hamiltonian H HF = (cid:88) k s c † k s M ( k ) c k s + N U m + N U , (A4) c k s = c k As c k Bs c k As c k Bs , M ( k ) = U m (cid:15) ( k ) 0 − t ⊥ (cid:15) ( k ) − U m − t ⊥ − t ⊥ U m (cid:15) ( k ) − t ⊥ (cid:15) ( k ) − U m . Here, (cid:15) ( k ) = (cid:15) +0 ( k ,
0) denotes the square lattice single-layer dispersion. The matrix M ( k ) is independent of thespin index s , and upon diagonalization, we arrive at fourbands α ∈ { , ..., } , (cid:15) α HF ( k , t ⊥ ) = ± (cid:113) ( U m ) + (cid:0) (cid:15) ± ( k , t ⊥ ) (cid:1) . (A5)Since U m corresponds to (half) the enery gap betweenthe upper and lower bands, we define the gap parameter∆ =
U m . The HF energy density E (∆) = (cid:104) H HF (cid:105) /N isthen expressed in terms of the gap as E (∆) = 1 N (cid:88) k αs f (cid:0) (cid:15) α HF ( k , t ⊥ ) (cid:1) (cid:15) α HF ( k , t ⊥ ) + ∆ U (A6)= − (cid:90) ∞ d(cid:15) ρ ( (cid:15), t ⊥ ) (cid:112) ∆ + (cid:15) + ∆ U , where f ( (cid:15) ) denotes the Fermi-Dirac distribution. To ob-tain the second line of Eq. (A6), we set T = 0 and ex- ploited the particle-hole symmetry of the system. In or-der to find the critical interaction U c /t , we expand theHF energy in ∆: E (∆) = E (0)+ ∆ U (1 − U χ )+ O (∆ ) , χ = (cid:90) ∞ d(cid:15) ρ ( (cid:15), t ⊥ ) (cid:15) (A7)The AF state will be favorable in energy, if the secondorder term becomes negative, i.e., at U c = χ − . Thus, inthe case of a diverging χ , any finite U leads to an AFinstability, i.e., U c = 0. And indeed, for a finite DOS at (cid:15) = 0, χ exhibits a logarithmic divergence, due to the1 /(cid:15) singularity of the integrand at (cid:15) = 0.Finally, we take a closer look at the functional behaviorof the gap ∆( U ). By demanding ∂E∂ ∆ = 0, we arrive atthe gap equation:1 = U (cid:90) ∞ d(cid:15) ρ ( (cid:15), t ⊥ ) √ ∆ + (cid:15) . (A8)In the single-layer system ( t ⊥ = 0), the gap equation canbe approximately solved and exposes a modified expo-nential behavior ∆ ∼ t exp( − π (cid:112) t/U ), while for finitevalues of t ⊥ , a numerical solution reveals a conventionalexponential scaling ∆ ∼ t exp( − C t ⊥ t/U ), cf. Fig. 10.The positive coefficient C t ⊥ increases continuously from C . t = 3 .
88 to C . t = 12 .
4, which upon approaching U → ∗ [email protected] † [email protected] A. Fuhrmann, D. Heilmann, and H. Monien, Phys. Rev. B , 245118 (2006). S. S. Kancharla and S. Okamoto, Phys. Rev. B , 193103(2007). K. Bouadim, G. G. Batrouni, F. H´ebert, and R. T. Scalet-tar, Phys. Rev. B , 144527 (2008). H. Hafermann, M. I. Katsnelson, and A. I. Lichtenstein,EPL (3), 37006 (2009). B .D. Napitu and J. Berakdar, Eur. Phys. J. B , 50(2012). S. Capponi, C. Wu, and S.-C. Zhang, Phys. Rev. B ,220505 (2004). A. Euverte, S. Chiesa, R. T. Scalettar, and G. G. Batrouni,Phys. Rev. B , 125141 (2013). L. Rademaker, S. Johnston, J. Zaanen, and J. van denBrink, Phys. Rev. B , 235115 (2013). L. Rademaker, J. van den Brink, J. Zaanen, andH. Hilgenkamp, Phys. Rev. B , 235127 (2013). R. R¨uger, L. F. Tocchio, R. Valent´ı, C. Gros, New J. Phys. , 033010 (2014). H. Lee, Y.-Z. Zhang, H. O. Jeschke, and R. Valent´ı, Phys.Rev. B , 035139 (2014). H. Zhai, F. Wang, and D.-H. Lee, Phys. Rev. B , 064517(2009). N. Bulut, D. J. Scalapino, and R. T. Scalettar, Phys. Rev. B , 5577 (1992). A. I. Liechtenstein, I. I. Mazin, and O. K. Andersen, Phys.Rev. Lett. , 2303 (1995). K. Kuroki, T. Kimura, and R. Arita, Phys. Rev. B ,184508 (2002). T. A. Maier and D. J. Scalapino, Phys. Rev. B ,180513(R) (2011). W. Cho, R. Thomale, S. Raghu, and S. A. Kivelson, Phys.Rev. B , 064505 (2013). S. Shimizu, H. Mukuda, Y. Kitaoka, A. Iyo, Y. Tanaka,Y. Kodama, K. Tokiwa, and T. Watanabe, Phys. Rev.Lett. , 257002 (2007). A. Damascelli, Z. Hussain, and Z.-X. Shen, Rev. Mod.Phys. , 473541 (2003). D. Fournier et al., Nat. Phys. , 905911 (2010). T. Esslinger, Annu. Rev. Condens. Matter Phys. , 129(2010). R. T. Scalettar, J .W. Cannon, D. J. Scalapino, andR. L. Sugar, Phys. Rev. B , 13419 (1994). W. Metzner et al., Rev. Mod. Phys. , 299 (2012). C. Honerkamp, Phys. Rev. Lett. , 146404 (2008). C. Platt, W. Hanke, R. Thomale, Advances in Physics ,453 (2013). K. Hida, J. Phys. Soc. Jpn. , 1013 (1992). A. J. Millis, and H. Monien, Phys. Rev. Lett. , 2810(1993). A. W. Sandvik, and D. J. Scalapino, Phys. Rev. Lett. , L. Wang, K. S. D. Beach, and A. W. Sandvik, Phys. Rev.B , 014431 (2006). C. Honerkamp and M. Salmhofer, Phys. Rev. B , 184516(2001). A. Eberlein and W. Metzner, Phys. Rev. B , 035126(2014). S. Uebelacker and C. Honerkamp, Phys. Rev. B , 235140(2012). M. Hohenadler, Z. Y. Meng, T. C. Lang, S. Wessel, A. Mu-ramatsu, and F.F. Assaad, Phys. Rev. B , 115132(2012). F. F. Assaad, and H. G. Evertz, Lecture Notes in Physics , 277 (2008). C. J. Halboth, and W. Metzner, Phys. Rev. B , 7364 (2000). C. N. Varney, C.-R. Lee, Z. J. Bai, S. Chiesa, M. Jarrell,and R. T. Scalettar, Phys. Rev. B , 075116 (2009). T. Sch¨afer, F. Geles, D. Rost, G. Rohringer, E. Ar-rigoni, K. Held, N. Bl¨umer, M. Aichhorn, and A. Toschi,arXiv:1405.7250 (2014). S. Chakravarty, B. I. Halperin, and D. R. Nelson, Phys.Rev. Lett. , 1057 (1988). F. D. M. Haldane, Phys. Rev. Lett. , 1029 (1988). M. Campostrini, M. Hasenbusch, A. Pelissetto, P. Rossi,and E. Vicari, Phys. Rev. B , 144520 (2002). S. A. Maier, A. Eberlein, C. Honerkamp, Phys. Rev. B ,035140 (2014). J. E. Hirsch, Phys. Rev. B31