Numerically exploring the 1D-2D dimensional crossover on spin dynamics in the doped Hubbard model
Y. F. Kung, C. Bazin, K. Wohlfeld, Yao Wang, C.-C. Chen, C.J. Jia, S. Johnston, B. Moritz, F. Mila, T. P. Devereaux
NNumerically exploring the 1D-2D dimensional crossover on spin dynamics in thedoped Hubbard model
Y. F. Kung , , C. Bazin , , K. Wohlfeld , , Yao Wang , , C.-C. Chen ,C.J. Jia , S. Johnston , , B. Moritz , , F. Mila , and T. P. Devereaux , Department of Physics, Stanford University, Stanford, CA 94305, USA Stanford Institute for Materials and Energy Sciences,SLAC National Accelerator Laboratory and Stanford University, Menlo Park, California 94025, USA Institute of Physics, Ecole Polytechnique Federale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland Institute of Theoretical Physics, Faculty of Physics,University of Warsaw, Pasteura 5, PL-02093 Warsaw, Poland Department of Applied Physics, Stanford University, California 94305, USA Department of Physics, University of Alabama at Birmingham, Birmingham, Alabama 35294, USA Department of Physics and Astronomy, University of Tennessee, Knoxville, TN 37996, USA Joint Institute for Advanced Materials, The University of Tennessee, Knoxville, TN 37996, USA Department of Physics and Astrophysics, University of North Dakota, Grand Forks, ND 58202, USA and Geballe Laboratory for Advanced Materials, Stanford University, Stanford, CA 94305, USA
Using determinant quantum Monte Carlo (DQMC) simulations, we systematically study the dop-ing dependence of the crossover from one to two dimensions and its impact on the magnetic prop-erties of the Hubbard model. A square lattice of chains is used, in which the dimensionality canbe tuned by varying the interchain coupling t ⊥ . The dynamical spin structure factor and staticquantities, such as the static spin susceptibility and nearest-neighbor spin correlation function, arecharacterized in the one- and two-dimensional limits as a benchmark. When the dimensionalityis tuned between these limits, the magnetic properties, while evolving smoothly from one to twodimensions, drastically change regardless of the doping level. This suggests that the spin excita-tions in the two-dimensional Hubbard model, even in the heavily doped case, cannot be explainedusing the spinon picture known from one dimension. The DQMC calculations are complemented bycluster perturbation theory studies to form a more complete picture of how the crossover occurs asa function of doping and how doped holes impact magnetic order. I. INTRODUCTION
A common feature of many strongly correlated elec-tron systems, such as the cuprate high-temperature su-perconductors, is layers of planes where most of the in-teresting physics occurs. This has motivated extensivestudies of two-dimensional models. However, quasi-one-dimensional materials whose internal crystalline struc-ture is known to be made of weakly coupled chains, suchas SrCuO , KCuF , and the organic Bechgaard salts, alsoexist and provide an alternative perspective on propertiessuch as magnetism. Indeed, the dimensionality of thesystem under consideration plays a crucial role in its be-havior. At the microscopic level, dimensionality impactsthe role of interactions. In two dimensions, electrons havea much larger number of paths to avoid one another thanin one dimension, where they have to interact. This dif-ference drastically modifies the physics, as single-particleexcitations can be described in terms of Landau quasi-particles in two dimensions but not in one. For example,on a chain, only collective spin and charge excitationsare possible, leading to spin-charge separation that hasbeen observed experimentally and has important conse-quences for the magnetic properties. Hence elucidatinghow the system changes as a function of dimensional-ity can provide a deeper understanding of the propertiesthemselves.The Hubbard model provides a simple, unified framework that describes one-, two-, and quasi-one-dimensional correlated electron systems, incorporatingthe effects of electron hopping and Coulomb interac-tions. It can be solved analytically in one dimension, withone of two approaches. Because the low-energy physicsof the Hubbard chain belongs to the one-dimensionalLuttinger liquid (LL) universality class (except for thecharge sector at half-filling), the first approach uses theapproximate bosonisation scheme to calculate the spec-trum of the model and show some of its most promi-nent properties, such as the aforementioned spin-chargeseparation. Low-energy properties and asymptotic cor-relation functions can also be evaluated. The second ap-proach uses the analytically exact Bethe ansatz, whichprovides a way not only to compute the spectrum of thesystem, but also to evaluate the LL parameters. To-gether, these two methods demonstrate that spin-chargeseparation and its consequences for physical observablesare the signature of one-dimensional physics. Yet, al-though the Hubbard chain can be solved exactly, obtain-ing information about static and dynamical correlationsstill relies strongly on numerical methods, and the spindynamics of the arbitrarily doped system has not beenexplored in detail. While the one-dimensional Hubbard model can besolved analytically, the two-dimensional case lacks exactsolutions and is not yet fully understood due to the com-plexity of the physics it describes, necessitating the use a r X i v : . [ c ond - m a t . s t r- e l ] J u l of numerical techniques such as quantum Monte Carlo(QMC), exact diagonalization, cluster pertur-bation theory, and dynamical mean-field theory. Computational studies have established that the half-filled system is a ( π/a, π/a ) antiferromagnetic (AF)insulator and that the ground state is a Mott insula-tor in the strong coupling limit. By the Mermin-Wagnertheorem, long-range order (LRO) cannot exist at finitetemperatures, implying that a gap opens and correla-tion functions are exponentially damped at ( π/a, π/a ).However, strong AF correlations are still present and in-fluence physical observables even at finite temperatures.Indeed, at ( π/a, π/a ), the static correlation functions ex-hibit a dominant mode and the spin excitation spectrumhas strong intensity. Hence, for simplicity, we will usethe term “LRO” when referring to cases in which the AFcorrelations extend across the finite-size cluster. Probingthe spin dynamics at half-filling corresponds to the cre-ation of one magnon, so the physics is approximated wellby linear spin-wave theory. A good understanding of thehalf-filled strongly coupled case is possible because spinsremain localized.When holes are doped into the system, though, theinterplay between AF order and hole delocalization com-plicates the situation. Upon doping, AF LRO rapidlydisappears and the doped system exhibits a large vari-ety of phases that compete or cooperate with one an-other. As a consequence, understanding the spin dy-namics of the doped two-dimensional Hubbard model re-mains a nontrivial task. In fact, it was only due to recentresonant inelastic x-ray experiments on doped cuprates,which suggested (rather surprisingly) that collective spinexcitations may persist up to high doping level in someregions of the Brillouin zone, that this problem wasstudied in greater detail. Because the Hubbard model exhibits markedly differ-ent behavior in one and two dimensions, examining howthe system crosses over between them provides insightsinto properties such as magnetism. The crossover canbe modeled with a system of coupled chains, where aninterchain coupling, t ⊥ , tunes a transition from effec-tively decoupled chains with confined electrons ( t ⊥ = 0)to a deconfinement of electrons throughout the lattice( t ⊥ = t for a two-dimensional system). Studies of quasi-one-dimensional systems use a variety of analytical andnumerical methods to calculate single-particle and two-particle processes. The analytical approaches generallyrely on renormalization group procedures and methodssimilar to the field theories used in one dimension, whilethe numerical methods include DMFT, QMC, and vari-ational cluster approximation (VCA). The interesting question is at which point the transi-tion from one- to two-dimensional character occurs. Be-cause the one-dimensional system is a Mott insulatorwith a gapped charge sector at half-filling but turns intoa LL when doped, the two cases must be studied sepa-rately. In the half-filled system with an intermediate U value ( ∼ t ), tuning t ⊥ can trigger a phase transition. A DMFT study has shown that for sufficiently small U and t ⊥ , a Fermi surface forms as t ⊥ increases. In ad-dition, sufficiently strong next-nearest-neighbor hoppingcan prevent spin-density waves from opening a gap in twodimensions when U has an intermediate value. Thus,the frustrated two-dimensional system is gapless at half-filling, and a metal-insulator transition occurs as t ⊥ isincreased, as shown by QMC. VCA and cluster DMFThave also studied the impact of the dimensional crossoveron Mott quantum criticality. When the system is doped, it has gapless charge sec-tors in both one and two dimensions. Renormaliza-tion group and perturbative approaches show that in-terchain single-particle motion is controlled by the pa-rameter α = ( K ρ + 1 /K ρ − /
4, where K ρ is the LLparameter. When α ≤ t ⊥ > When α >
1, the particlesremain confined to the chains. Analytical results sug-gest that the Hubbard model with finite U should have α < / t ⊥ . How-ever, numerical studies demonstrate that the situation ismore complicated. A QMC study finds that electronscan be confined for intermediate values of α smaller than1. Another study shows that as t ⊥ increases, coher-ent interchain motion develops, and the spectral functionevolves from a LL form with decoupled chains and spin-charge separation towards a Fermi-liquid-like one withtwo-dimensional character and well-defined quasiparticlepeaks. However, for intermediate values of t ⊥ , LL fea-tures remain present at high energies. Similarly, a DMFTstudy demonstrates that with an intermediate U value( U = 4 t ), lowering the temperature can trigger a tran-sition from a LL to a Fermi liquid for t ⊥ > Thesenumerical results show that LL features can be observedeven for finite interchain coupling.Thus far, studies of magnetism in the dimensionalcrossover regime have focused on the half-filled system.As discussed already, the half-filled two-dimensional, butnot one-dimensional, system is an AF insulator that de-velops LRO at zero temperature. A question of greatinterest is the value of t ⊥ necessary to induce this LROin the chain. As yet there is no definitive answer, buta renormalization group study and multiple numeri-cal studies of the anistropic Heisenberg model, aswell as a QMC study of the intermediate- U Hubbardmodel, have suggested that any t ⊥ > How-ever, despite these studies, the interplay between dopingand magnetism has not yet been examined in the frame-work of the dimensional crossover.The goal of this study is to shed light on the mag-netism of the strongly correlated doped Hubbard modelas it transitions from one to two dimensions. In additionto providing a means of comparison between the dopedmagnetism and magnetic excitations in one and two di-mensions, the effect of the dimensionality on the dopedmagnetism of the Hubbard model is interesting in itself.Hence, the magnetic properties are examined on a latticeof coupled chains where the interchain coupling is varied,in order to build upon previous results in one and two di-mensions, and to elucidate magnetic properties such asthe spin dynamics that are not yet well understood. Assuggested in an earlier study, short-range spin correla-tions can also provide insight into the effect of doping onboth magnetic order and spin excitations. These quanti-ties are computed via DQMC and the maximumentropy method (MaxEnt) of analytic continuation. The paper is organized as follows. In Section II, theHubbard model is presented together with the numeri-cal methods used to carry out the simulations. SectionIII discusses the static and dynamic spin properties inone and two dimensions as a benchmark before SectionIV explores the doping dependence of the dimensionalcrossover. Finally, Section V summarizes the main out-comes of this study and discusses further perspectives.
II. MODEL AND METHODS
The single-band Hubbard Hamiltonian de-scribes strongly correlated electrons on a lattice: H = (cid:88) (cid:104) i,j (cid:105) σ t ij ( c † iσ c jσ + h.c. ) − µ (cid:88) iσ n iσ + U (cid:88) i (cid:18) n i ↑ − (cid:19)(cid:18) n i ↓ − (cid:19) , (1)where c † iσ ( c iσ ) creates (annihilates) a particle with spin σ on site i , and n iσ = c † iσ c iσ is the number operator.The nearest-neighbor hoppings along the same chain andbetween chains are controlled by t ij ≡ t and t ij ≡ t ⊥ , re-spectively. The longer range hoppings are all set to zeroexcept in the two-dimensional system, for which we mayalso consider the case with a finite next-nearest neighborhopping t ij ≡ t (cid:48) . U is the on-site Coulomb interactionthat penalizes double occupancy, and a = 1 is the unitof length. We work with U = 8 t , so the ground stateis a strongly correlated Mott insulator in the undopedsystem, and measure energies in units of t . The chem-ical potential µ is adjusted to give the desired doping.The model exhibits particle-hole symmetry, and the holedoping level can be defined as p = 1 − n , where n is theelectron density.In this study, properties of the Hubbard model arecalculated using DQMC, a numerically exact,auxiliary-field technique that computes observables fromimaginary-time Green’s functions as (cid:104) ˆ O (cid:105) = tr[ ˆ Oe − βH ]tr[ e − βH ] , (2)with the imaginary-time interval [0 , β ] divided into L slices of width ∆ τ . The Hamiltonian can be rewritten in terms of the non-interacting and interacting pieces,and the exponential decomposed using the Trotter ap-proximation e − L ∆ τH ≈ ( e − ∆ τK e − ∆ τV ) L , (3)where K contains quadratic terms and V is the quar-tic interaction term. Terms in the expansion of order O (∆ τ ) and higher are dropped. In this study, a suffi-ciently small time slice was used such that no significant∆ τ errors were found.A Hubbard-Stratonovich (HS) transformation e − ∆ τU ( n i ↑ − )( n i ↓ − ) = 12 e U ∆ τ/ (cid:88) s i,l = ± s i,l e − ∆ τλs i,l ( n i ↑ − n i ↓ ) , (4)is used to rewrite V in quadratic form, at the cost ofintroducing a new HS field s i,l = ± i and time slice l . The relation cosh (∆ τ λ ) = exp (∆ τ U/ λ . The partition function can now be calculatedas Z = (cid:88) s i,l = ± det M + det M − , (5)where M σ = I + B σL B σL − ...B σ (6)and B ± l = e ∓ ∆ τλv ( l ) e − ∆ τK , (7)and v ( l ) is a diagonal matrix with s i,l as the i th ele-ment. The Monte Carlo sampling is performed over theHS field configurations, each of which has a weight of P ( s ) = det M + det M − /Z . This is used to computethe Green’s function, which is in turn used to com-pute all other quantities via Wick’s theorem. Since allobservables are calculated in terms of imaginary time,they must be analytically continued to real frequenciesfor comparison to experiments. In this study we em-ploy MaxEnt, which uses Bayesian statistical inferenceto determine the most probable spectral density givenan imaginary-time correlator. DQMC has the advantages of being numerically ex-act and of accessing relatively large system sizes, but ingeneral it suffers from a fermion sign problem.
Be-cause the algorithm does not track the order of the oper-ators, a negative sign from the fermion anticommutationrelations remains undetermined and all observables mustbe divided by the average fermion sign (cid:104) f sgn (cid:105) as (cid:104) ˆ O (cid:105) = (cid:80) s m,l ˆ Of sgn P ( s ) (cid:80) s m,l f sgn P ( s ) = (cid:104) Of sgn (cid:105)(cid:104) f sgn (cid:105) ,f sgn = sign(det M + det M − ) . (8) FIG. 1: The doping dependence of the average fermion sign isshown for different values of interchain coupling at β = 3 /t . Statistical fluctuations become more significant as theaverage sign decreases; hence its value controls accessibleparameter regimes.Figure 1 systematically explores the average fermionsign for different doping levels and interchain hoppings.At half-filling, particle-hole symmetry protects the signsuch that it is always 1. Away from half-filling, however,the average sign is suppressed exponentially as the tem-perature decreases and the system size increases. De-spite doping, the average sign remains close to 1 in onedimension due to the small number of available hoppingpathways (and opportunities for ambiguities in sign). Asthe interchain coupling t ⊥ and dimensionality increase,the system smoothly evolves from one to two dimensions,where the average sign is strongly reduced upon doping.Doped two-dimensional systems require higher tempera-tures to have a non-vanishing average sign, constrainingour simulations to β = 3 /t in general. III. SPIN PROPERTIES IN ONE AND TWODIMENSIONSA. One-Dimensional Case
We begin with a review of the one-dimensional case,systematically studying its static and dynamic magneticproperties for comparison to the two-dimensional case.The static spin properties elucidate the effects of dop-ing on magnetic order and short-range correlations. Thepeak position of the static spin susceptibility χ ( q ) shiftswith doping as 2 k F = nπ . The peak intensity is highestat half-filling, q = π , suggesting strong AF correlations,and decreases upon doping, indicating that the 2 k F spin-density wave weakens as its mode changes. This weak-ening is confirmed by a decrease in the magnitude ofthe spatial spin-spin correlation function (cid:104) S i · S j (cid:105) upondoping. Finite temperature destroys the quasi-LRO, im-plying an exponential decay of correlations with a fi-nite correlation length ξ , which can be extracted from (cid:104) S i · S j (cid:105) . At half-filling, ξ is in good agreement withthe Bethe ansatz. The small correlation length suggeststhat short-range correlations play an important role inthe observed magnetic properties.Short-range correlations can be examined via thenearest-neighbor (NN) spin correlation (cid:104) S S (cid:105) , whichshows that low levels of hole doping reduce the local den-sity of spins without additionally destroying magnetic or-der, as is consistent with the importance of short-rangecorrelations upon doping. However, for higher dopinglevels, the NN correlation decreases more rapidly thanexpected from a local static picture. Naively, one wouldthink that once most electrons have been removed, holedelocalization should have a reduced impact on the mag-netic order. However, spin-charge separation complicatesthe situation, as doping a hole corresponds to addingboth a holon and a spinon. They can be thought ofas acting independently (due to spin-charge separation);thus the holon would delocalize without affecting (anddestroying) the magnetic order, explaining the behaviorof the NN spin correlations up to intermediate dopinglevels. However, spinons may act like domain walls andcause additional destruction of magnetic order. This pic-ture would be consistent with the faster decay of the NNspin correlations at large doping. The transition meansthat the way in which doped holes destroy magnetic orderevolves as electron density decreases.The dynamical spin structure factor, S ( q , ω ), provides FIG. 2: False-color plots of the dynamical spin structure fac-tor are shown at (a) n = 1 .
0, (b) n = 0 .
8, and (c) n = 0 . U = 8 t and β = 15 /t to access the expected low-temperature behavior. energy- and momentum-resolved information about mag-netic properties. Fig. 2(a) shows that the half-filled sys-tem exhibits a continuum, confirming that a spin flipdecays into independent spinons, as expected in the pres-ence of spin-charge separation. The spectral intensityis well described by the two-spinon continuum computedfrom the Bethe ansatz for the Heisenberg model (solidlines), with spectral weight concentrated at the lowerboundary due to suppression of itinerancy effects anddominant spin excitations at low energy due to virtualelectron hopping for larger U . The spectrum is broad-ened by finite temperature, which also destroys quasi-LRO and opens a small gap at π . Nevertheless, the highspectral intensity at 2 k F = π confirms the presence ofstrong AF correlations. In addition, despite finite-sizeand finite-temperature effects, the spin velocity arisingfrom the LL formalism agrees well with the Bethe ansatzprediction, suggesting that at half-filling, the spin dy-namics in the strong coupling limit is described well bythe Bethe ansatz solution of the Heisenberg chain at zerotemperature.In the doped system [Figs. 2(b) and (c)], spectral in-tensity both decreases and broadens, while the spectrumhardens overall. In addition, the damping of intensityis the largest at π . The continuum indicating the pres-ence of spin-charge separation still can be seen, suggest-ing that the doped chain also exhibits characteristics ofa LL, which is further confirmed by the soft mode at2 k F = nπ (where n is the electron density). This shiftin 2 k F away from π with doping also explains why thespectrum becomes increasingly gapped at k = π : It nolonger corresponds to the soft mode of the spin-densitywave of the doped system. As in the half-filled case, thespin velocities can be extracted and are found to be ingood agreement with theoretical values. The qualitative behavior of the doped spectra is verydifferent from that of the two-spinon continuum, whichis symmetric about 2 k F . As four-spinon processes areexpected to contribute more as doping increases, thisqualitative observation seems to confirm that processesother than two-spinon ones are involved in the spin dy-namics and are enhanced by doping. Moreover, as op-posed to the half-filled case, they do not seem to overlapsimply with the two-spinon continuum.In order to explain the overall hardening of the spin ex-citation spectrum upon doping, naively one could studyit in terms of the local static picture. Although there is noLRO in one dimension, the robustness and the relativelystrong intensity of the continuum at high doping levelsare likely related to strong short-range correlations. B. Two-Dimensional Case
Extensive DQMC studies have already been performedto characterize static and dynamic magnetic propertiesin two dimensions, which we review briefly. Both thestatic spin susceptibility and real-space spin-spin corre- lation function provide evidence for the presence of AForder at half-filling and its destruction with doping, asexpected. The NN spin correlation (cid:104) S · S (cid:105) explores the inter-action of doped carriers with the magnetic background.When t (cid:48) = 0, a discrepancy in (cid:104) S · S (cid:105) from whatwould be expected in the simple local static pictureis observed up to intermediate hole-doping levels. Infact, this result is not surprising, since in the spin po-laron picture , introducing holes into a magnetically or-dered background causes a strong reduction of magneticcorrelations that is larger than the simple reduction inmagnetic moments. However, once the number ofdoped holes exceeds 50%, the density of electrons is smallenough for the holes to delocalize without breaking AFbonds. Then hole delocalization no longer competes withmagnetic correlations and the NN spin correlation corre-sponds to the local static picture.A distinct situation occurs once the next-nearest-neighbor (NNN) hopping t (cid:48) = − . t is switched on (notethat such a value of the NNN hopping is typically chosenso that the Hubbard model constitutes a more realisticdescription of the cuprates). On the hole-doped side,the effect of t (cid:48) is to favor sublattice mixing, which leadsto an enhanced destruction of AF order. On the electron-doped side, finite t (cid:48) supports AF correlations, so thatthe local static picture can explain the hardening of thespin excitation spectrum. The latter, rather surprising,result strongly depends on the choice of t (cid:48) : It can only beobtained for strong enough t (cid:48) , i.e. a negative but muchsmaller value | t (cid:48) | (cid:28) . t would not be enough.Figure 3 shows S ( q , ω ) at different doping levels, withdots indicating the maxima of spectral intensity. Despitethermal broadening and renormalizations from quantumeffects, in general the spin excitation spectrum exhibitsthe features expected at half-filling [Fig. 3(a)]. The in-tensity maxima reproduce the linear spin-wave disper-sion (calculated for a nearest-neighbor Heisenberg modelwith J = 4 t /U = 0 . t for U = 8 t ) up to a multiplicativefactor. However, ( π, π ) has the strongest intensity, indi-cating that AF order builds up at half-filling despite theabsence of true LRO in the numerical simulation.As the doping level increases [Figs. 3(b) and (c)],the spectrum hardens and spectral intensity decreases.These effects are most pronounced at ( π, π ), which isaffected by the destruction of AF order. Along the di-rection (0 , → ( π, t (cid:48) case reportedin Ref. 24. Moreover, the spectrum hardens along the(0 , → ( π,
0) direction. This is a counterintuitive re-sult, since an overall softening has been suggested in theliterature.
However, a study comparing spin sus-ceptibilities calculated by DQMC and the random phaseapproximation, which was developed for weakly interact-ing systems, demonstrated that this discrepancy may becaused by significant correlations that persist to higherdoping levels. FIG. 3: False-color plots of the dynamical spin structure fac-tor S ( q , ω ) along high-symmetry directions are shown at (a) n = 1 .
0, (b) n = 0 .
8, and (c) n = 0 . ×
10 cluster with U = 8 t and β = 3 /t . The spectraat n = 0 . n = 0 . J = 4 t /U = 0 . t . IV. MAGNETIC PROPERTIES IN THEDIMENSIONAL CROSSOVER
Depending on the dimensionality, the magnetic re-sponse of the Hubbard system can be very different –e.g. the spin dynamical structure factor shows the onsetof spinon continua in one dimension, whereas surpris-ingly stable S = 1 collective spin excitations (magnons)are observed in two dimensions, even at high doping lev-els. Here we present results for the crossover from oneto two dimensions and explore the doping dependenceof magnetic properties as a function of dimensionality.Unless otherwise noted, calculations are performed on a10 ×
10 cluster with U = 8 t , β = 3 /t , and interchain cou-pling values t ⊥ /t = 0 . , . , . , . , . , . , . t ⊥ = 0. Thus, physical properties only have spatialdependence along the chain and are independent of thetransverse direction. In reciprocal space, the direction(0 , → ( π,
0) effectively corresponds to a single chain.Due to the transverse-direction independence, any prop-erties along (0 , k ) → ( π, k ) will be identical to the thosealong the chain, and properties along ( k, → ( k, π )will be identical to those at ( k, A. Static Properties
First we study how the static spin properties evolveupon tuning t ⊥ . As before, the static spin susceptibilityprovides an energy-integrated perspective, while the NNspin correlation function shows how doped holes interactwith the magnetic background.
1. Static Spin Susceptibility
Figure 4 shows the static spin susceptibility at half-filling and 40% hole doping for different values of t ⊥ . For t ⊥ < . t , the susceptibility only has momentum depen-dence along the chains and remains almost flat along thetransverse direction. This again indicates that the mag-netic order mostly retains one-dimensional character forsmall t ⊥ .Figure 4(a) shows that at half-filling, magnetic orderis transferred smoothly from ( π, π, π ) as inter-chain coupling increases. As soon as t ⊥ >
0, the staticsusceptibility is larger at ( π, π ) than at ( π, t ⊥ , the intensity at ( π, π ) is strongest, showing that AForder dominates and suggesting that increasing t ⊥ leadsto spinon confinement and AF LRO. Therefore, as hasbeen deduced from the transfer of the spectral intensityin the spin dynamics, the system smoothly develops AForder as the chains are increasingly coupled.In comparison with Fig. 4(a), the AF order at 40%doping [Fig. 4(b)] has been reduced strongly. In addition,the susceptibility at wave vectors that dominate in bothone and two dimensions at half-filling are weakened sig-nificantly by doping. Along the (0 , → ( π, → ( π, π )direction, the soft mode in the spin dynamics expectedin one dimension at q = (2 k F ,
0) is not observed becauseof the high temperature. Nevertheless, for t ⊥ ≤ . t , thestatic spin susceptibility has a peak along the chain corre-sponding to the 2 k F one-dimensional mode. This obser-vation clearly shows that for small enough t ⊥ , the chainsretain their one-dimensional character at 40% doping.We note that this 2 k F mode is directly related tothe soft mode expected in the spin dynamics. Thus, FIG. 4: Static spin susceptibility χ ( q ) for different values of t ⊥ at (a) half-filling and (b) 40% hole doping. The mag-netic order is transferred from one to two dimensions. Thecrossover only occurs for t ⊥ ≥ . t because of the high tem-perature. the static susceptibility provides us with a quantita-tive means of tracking the evolution from one- to two-dimensional character in the doped system. Indeed, as t ⊥ increases above 0 . t , the spectral intensity of the one-dimensional-like peak at (2 k F ,
0) shifts smoothly toward( π, π − q ), which is expected in two dimensions for highdoping.Therefore, studying the static spin susceptibility al-lows us to track more closely the evolution of the one-dimensional properties as t ⊥ is varied in the doped sys-tem. At half-filling, low temperatures can be reachedso that most physical observables can be used to trackthe dimensional crossover. However, in the doped case,the sign problem restricts simulations to higher temper-atures, making it difficult to observe one-dimensional be-havior. Thus the static susceptibility enables us to quan-tify the dimensional crossover. FIG. 5: NN spin correlation as a function of electron density n for different values of the interchain hopping t ⊥ , shown bothalong (top row) and perpendicular to (bottom row) the chain.Perpendicular to the chain, the spin correlation function for t ⊥ = 0 . t is shown instead of that for t ⊥ = 0 t , which is identi-cally zero. The values are very small close to one dimension,leading to larger error bars. The solid lines show the (1 − p )curve predicted by the local static picture.
2. Nearest-Neighbor Spin Correlations
As already discussed, dimensionality has a strong im-pact on the way in which hole delocalization interactswith magnetic order as a function of doping. In one di-mension, hole delocalization destroys magnetic order byreducing spin density at low doping levels but causes agreater destruction at high doping levels. In two dimen-sions, the opposite trend has been observed. To examinehow the behavior crosses over between these two lim-its, the longitudinal (along the chain) and transverse NNspin correlation functions are calculated as a function ofdoping for different values of t ⊥ , as shown in Fig. 5.Although the doping trend of the longitudinal correla-tions [Fig. 5 (top row)] interpolates smoothly from one totwo dimensions, the transverse correlations [Fig. 5 (bot-tom row)] show an unexpected trend. After a rapid de-crease to a plateau, upon increasing doping the transverseNN spin correlation function recovers, and its magnitudeeven exceeds the prediction from the local static picture.This plateau suggests a transitory regime where dopedholes induce processes that compensate the reduction ofmagnetic order caused by local spin density reduction. Itcould be related to the spin excitations along the trans-verse direction: As shown in Fig. 7, the dispersion ini-tially lifts up as the system is hole doped, but does notevolve significantly upon further doping (up to 40%).To check whether this doping dependence is a ther-mal effect, the same measurements could be performed atlower temperatures, decreasing t ⊥ to ameliorate the signproblem. We note that it would be surprising for hightemperatures to enhance correlations above the simplelocal spin density reduction curve. Thus, the doping de-pendence of the NN spin correlation function is intriguingand highlights the importance of short-range correlationsin understanding the properties of the doped system. B. Spin Dynamics: Half-Filling
Figure 6 shows the dynamical spin structure factor, S ( q , ω ), for different values of t ⊥ at half-filling. The solidlines correspond to the boundaries of the two-spinon con-tinuum for t ⊥ = 0, and to the linear spin-wave dispersionfor t ⊥ > ω ( k x , k y ) = (cid:113) ( J + J ⊥ ) − ( J cos ( k x ) + J ⊥ cos ( k y )) , (9)where J = 4 t /U and J ⊥ = 4 t ⊥ /U . The spectra arebroad due to the high simulation temperature, so dotsindicating the intensity peaks provide a guide to the eye.When t ⊥ = 0 [Fig. 6(a)], the system consists of 10-sitechains. At β = 3 /t , the two-spinon continuum is barelydistinguishable along (0 , → ( π, β = 15 /t (notshown), indicating that the 10-site chain is long enoughto capture one-dimensional spin dynamics at sufficientlylow temperatures.As t ⊥ increases [Figs. 6(b)-(f)], magnetic correlationstransition from the ones known for the chain to thoseknown for the two-dimensional lattice. Spectral inten-sity shifts from ( π,
0) toward ( π, π ), suggesting the tran-sition to AF LRO that is predicted at zero temperaturefor small transverse hopping. In addition, peaks in thespectral intensity increasingly follow the linear spin-wavedispersion (up to a multiplicative factor), suggesting thatcoherent spin-wave excitations have replaced indepen-dent spinons, which have become confined. The only ex-ception is along the direction (0 , → (0 , π ), where thespectral intensity does not follow the linear spin-wavedispersion up to large values of t ⊥ . Since low-energyspin excitations are present along this direction, strongthermal fluctuations probably impede the development oftwo-dimensional magnetism. This transition from one-to two-dimensional magnetism was observed previouslyin Sec. IV A in the static spin susceptibility.For small values of t ⊥ , the spectrum mostly retainsone-dimensional characteristics, such as the two-spinoncontinuum and weak momentum dependence along trans-verse directions. This behavior is due to the simula-tion temperature, as high temperatures partially washout the increase in dimensionality. Indeed, coherentspin waves form between t ⊥ = 0 . t [Fig. 6(c)] and t ⊥ = 0 . t [Fig. 6(d)] along certain directions. In ad-dition, the spectral intensity along the transverse direc-tion (0 , → (0 , π ) begins to disperse more strongly at t ⊥ = 0 . t .The effect of the dimensional crossover on the physicalobservables depends on the energy scale, as also reportedin experiments. For example, the cuprates exhibit a spin-wave dispersion at low energies, and a spinon con-tinuum at higher energies. This transition temperaturecan be estimated naively as T transition ∼ t ⊥ ; however, thetransition affects single- and two-particle interchain pro-cesses differently and is renormalized by interactions sothat the effective crossover temperature is lower. Despitethe Coulomb interaction in this calculation, the naivenon-interacting prediction of T transition ∼ t ⊥ = t/ β = 3 /t . The temperature dependenceof the transition is confirmed by performing a simulationat β = 5 /t for t ⊥ = 0 . t . As expected, spectral intensityat ( π, π ) is enhanced and the positions of the peak inten-sity track the linear spin-wave dispersion more closely at β = 5 /t than at β = 3 /t . At the energy scales probedhere, this crossover occurs smoothly. These observationsagree well with previous studies. C. Spin Dynamics: Doping and Dimensionality
In one dimension, the spin excitation spectrum hard-ens upon doping and develops a soft mode at 2 k F . Intwo dimensions, the spectrum exhibits an overall hard-ening but with a relatively persistent intensity along theaxes in reciprocal space. In order to study the effect ofdimensionality on the doping dependence of the spin exci-tation spectra, Fig. 7 shows the dynamical spin structurefactor for different doping levels and interchain couplingstrengths.In one dimension, due to high temperature, the softmode at 2 k F is not seen even though 20% and 40% holedoping correspond to 2 k F = nπ = π and π , respectively.However, when the temperature is decreased to β = 15 /t (not shown), the soft mode is observed, suggesting againthat one-dimensional physics is present on a 10-site dopedchain even when obscured by high temperatures. In fact,as discussed below, the 2 k F mode is seen in the staticspin susceptibility even at β = 3 /t .The temperature has the same effect on the crossoverin the doped system as in the half-filled case, with thetransition between one- and two-dimensional magnetismoccurring between t ⊥ = 0 . t and t ⊥ = 0 . t . Indeed,for t ⊥ < . t , the energy scales of the spectra are al-most unaffected and the momentum dependence alongtransverse directions of reciprocal space remains small.Generally, the dimensional crossover occurs smoothly inthe spin dynamics, independent of the doping level, withthe spectra transitioning gradually from the continuumin one dimension to showing more coherent spin excita-tions in two dimensions. This suggests that the natureof the spin excitations gradually changes from spinon-likebehavior at small t ⊥ to a magnon-like response at large t ⊥ . We will explore specific aspects of the crossover ingreater detail in the following sections. FIG. 6: False-color plot of the dynamical spin structure factor, S ( q , ω ), along the main symmetry directions at half-filling, fordifferent values of t ⊥ . In one dimension, the solid line corresponds to the two-spinon continuum; for t ⊥ (cid:54) = 0, it represents thelinear spin-wave dispersion. The dots indicate the maximum intensity at each momentum point (in units of π ).
1. Transverse Dispersion
Along the transverse direction, (0 , → (0 , π ), finite t ⊥ causes the dispersion to lift up upon doping and the spec-tral intensity to spread more. Doping the chains enhancesthe spin excitations and hence their two-dimensionalcharacter (they are dispersionless in one dimension butdisperse in two dimensions). The effect of doping on thedimensionality dependence is clearly seen in Fig. 7. Athalf-filling, the spectral intensity peaks show a weak dis-persion away from the linear spin-wave dispersion up tolarge values of t ⊥ , but at high doping levels, spectralhardening is more sensitive to dimensionality, even forfairly small t ⊥ .A simple local picture offers insights into the enhanceddispersion upon doping. Once holes are doped, the t − J three-site term, in addition to the Heisenberg term, con-tributes to the spin dispersion. In the transverse direc-tion, these two terms depend differently upon t ⊥ . Forthe Heisenberg exchange, the magnetic coupling involvesa virtual hopping between nearest neighbors only, so J ⊥ = 4 t ⊥ /U and is suppressed for small t ⊥ . On the otherhand, there are three possible channels for the three-site term, which involves a virtual hopping between thenearest and next-nearest neighbors. In one of the chan-nels, the coupling between spins depends linearly on t ⊥ : J ⊥ = 4 tt ⊥ /U . Thus, at half-filling, the spin exchangescales as t ⊥ and is reduced, while in the doped case itscales as t ⊥ . This linear dependence may explain theenhancement of the spin dispersion along (0 , → (0 , π ). The evolution of the peak intensity position of the dy-namical spin structure factor at (0 , π ) as a function of t ⊥ is shown in Fig. 8(a) for different doping levels. Near half-filling, the evolution is closer to t ⊥ . At high doping levels,the dependence becomes more linear even though it is ob-scured by variations from thermal broadening, which in-creases the variability of the peak position. Nevertheless,a clear transition occurs in the t ⊥ dependence betweenthe half-filled and the doped cases.At half-filling, for small t ⊥ , the energy scale of trans-verse spin excitations is too small compared to the tem-perature to observe a dispersion along (0 , → (0 , π ).However, once holes are doped, the contribution of thethree-site term induces higher-energy excitations that aremore easily probed at high temperatures, showing theimportance of the term in understanding the physics ofdoped systems. It also suggests that doping the systemmakes the crossover of the low-energy part of the spin ex-citation spectrum occur at a smaller interchain couplingvalue. Thus, an anistropic lattice allows us to disentangledifferent processes and separate their contributions.
2. Momentum-Dependent Crossover and Persistent SpinExcitations
The evolution of the spin dynamics at ( π, π ) and ( π, t ⊥ ≥ . t , their energy scales0 FIG. 7: False-color plot of the dynamical spin structure factor S ( q , ω ) along the main symmetry directions, for different dopinglevels and values of the interchain coupling, t ⊥ . At half-filling in one dimension, the solid line corresponds to the two-spinoncontinuum; for t ⊥ (cid:54) = 0, it represents the linear spin-wave dispersion. The dots indicate the maximum intensity for eachmomentum point, which is given in units of π . can be described with a linear spin-wave dispersion. At( π, π ), AF order sets the energy scale, which has no sig-nificant t ⊥ dependence and remains nearly gapless. At( π, ω ( π,
0) = 2 √ JJ ⊥ ∼ t ⊥ .Although linear spin-wave theory can predict the half-filled behavior, it does not describe the 40% hole-dopedsystem. As at half-filling, the energy scale of the ( π, π )spin excitations has only a small t ⊥ dependence. In onedimension, ( π, π ) is equivalent to ( π,
0) and the hard-ening of the spectrum at ( π, π ) thus corresponds to thehardening of the spectrum of the chain upon doping. Itis remarkable that despite the very different natures ofthe spin excitations in one and two dimensions, the en-ergy scale of this point remains the same as interchaincoupling is increased. On the other hand, the ( π,
0) spinexcitations soften as t ⊥ increases. Hence, as the spinonsstart to bind together, the energy cost of a spin flip alongthe chain direction decreases, and the energy scale of thespectrum smoothly interpolates from one to two dimen- sions. Therefore, there is almost no crossover in the en-ergy scale of the spin excitations at ( π, π ) while one existsat ( π,
3. Hardening of the Spectrum
In two dimensions, the spectrum has been observedto harden upon doping, along the (0 , → ( π,
0) direc-tion. In fact, as shown in Fig. 7, this behavior occurs forany given value of t ⊥ . Moreover, the energy scale of thehardening is always larger than that in two dimensions(ranging from ∼ . t in two dimensions to ∼ t in onedimension). Although the local static picture does notfully apply, it can still shed light on why the spectrumhardens more for t ⊥ < t . For a given value of t ⊥ , theenergy costs of a single local spin flip in the doped andundoped cases can be compared to determine whetherthey explain the greater hardening observed for small t ⊥ . In the undoped case, the Heisenberg model gives the1 FIG. 8: Energy of the spin excitations at different momenta asa function of t ⊥ for different doping levels. Panel (a) focuseson the detailed doping dependence at k = (0 , π ), while panel(b) compares the evolution of the energy of spin excitationsat k = ( π,
0) and k = ( π, π ) at two different doping levels. following energy cost for a spin flip:∆ E Undoped = 4 U (cid:0) t + t ⊥ (cid:1) . (10)When one hole is doped on the NN site of the flippedspin, whether the doped hole is on the same chain as theflipped spin or on its neighboring chain must be takeninto account.While the Heisenberg model only includes terms of or-der t and t ⊥ , the three-site term includes three differentchannels: ∆ E Doped = 18 4 U (cid:0) t ⊥ + 4 tt ⊥ + 7 t (cid:1) . (11)Thus, the local picture predicts that for sufficiently small t ⊥ , ∆ E Doped < ∆ E Undoped , so the spectrum shouldsoften upon doping and harden only for large enough t ⊥ . Moreover, as t ⊥ gets closer to the two-dimensionallimit, the spectrum should harden more. This is clearly very different from what is observed in Fig. 7, where thespin excitations harden less and less upon doping as t ⊥ increases.The local picture may fail because of the high tem-perature. As high-energy spin excitations in the low- t ⊥ regime retain their collective one-dimensional properties,a picture of local magnon creation does not apply. Inorder to examine the prediction of the transition from asoftening to a hardening of the spin excitations, a futurestudy could examine the small- t ⊥ regime at lower tem-peratures. If this transition is observed, it will validatethe local picture. D. Comparison with Cluster Perturbation Theory
The DQMC studies of the dimensional crossover in thespin excitation spectrum can be complemented with clus-ter perturbation theory (CPT) calculations. This nu-merical technique combines exact diagonalization (ED)and perturbation theory, dividing the infinite plane intosmaller identical clusters that are solved exactly usingED. Hopping between the clusters is treated to lead-ing order in perturbation theory. CPT is exact in thelimits of strong and weak coupling as the number ofBrillouin zone sites L → ∞ . Unlike DQMC, it isgenerally performed at zero-temperature, thus avoidingfinite-temperature effects. It also complements DQMCwith its fine momentum resolution. However, becausethe ED solver works in the canonical ensemble, dop-ings are limited to discrete levels, as opposed to thecontinuous doping evolution accessible to DQMC, whichworks in the grand canonical ensemble. In this section,the CPT simulations are performed with U = 8 t and t ⊥ /t = 0 . , . , . , . , . , . , . , . , . , . , .
0. Boththe 12-site C4 and 2 × t ⊥ /t = 0 , . , . t ⊥ >
0, the numerical calculation agreeswell with the linear spin-wave dispersion (dashed lines)throughout the Brillouin zone, again similar to the be-havior of the DQMC calculation. The only significantdiscrepancy between the CPT and DQMC spin excita-tion spectra occurs in the quasi-one-dimensional systemat (0 , π ) (the Y-point), where the dynamical spin struc-ture factor from CPT follows the spin-wave dispersionwhile that from DQMC has spectral intensity at lower en-ergies up to large values of t ⊥ . As DQMC is performed ata significantly higher temperature, the difference is mostlikely due to a thermal effect.When the system is doped away from half-filling, theCPT and DQMC calculations continue to agree well. Re-gardless of doping level, the dimensional crossover occursin a smooth transition as t ⊥ is varied. When the in-2 FIG. 9: False color plot of the dynamical spin structure factor S ( q , ω ) along the main symmetry directions, at half-filling (leftcolumn), 16 .
7% (center column), and 33 .
3% (right column) hole doping, for different values of the transverse hopping integral t ⊥ . The first row corresponds to t ⊥ = 0, the second to t ⊥ = 0 . t , and the third to t ⊥ = t . The color scale is the same forall plots. At half-filling, the solid line corresponds to the two-spinon continuum and the dashed line to the linear spin-wavedispersion. terchain hopping is increased, the transverse dispersion(Γ → Y ) hardens, and spectral weight shifts from ( π, π, π ) as two-dimensional character is enhanced.In addition, increasing hole doping causes the spectra toharden along the longitudinal direction (Γ → X ). Theclose agreement of results from these two techniques sug-gests that they access the same physics despite the dif-ference in simulation temperatures.Section IV has examined the dimensional crossoverin half-filled and doped systems using DQMC, comple-mented by CPT calculations. The evolution of the dopingdependence with increasing interchain coupling has beenexplored systematically. The crossover appears to occursmoothly for all doping levels and momentum points. Bycomparing the momentum dependence of the crossoverin the spin dynamics at half-filling and at 40% dop-ing, the presence of persistent magnon-like excitationshas been confirmed. Moreover, the importance of thethree-site term in understanding spin dynamics has beenhighlighted. The mechanism behind spectral hardeningupon hole doping remains unclear, however; future stud-ies could be performed at lower temperatures to studythe applicability of the local picture. The doping trendof the NN spin correlation function as the system evolvesfrom one to two dimensions also sheds light on how dopedholes interact with the magnetic order. Finally, a com-parison of spin dynamics calculated by DQMC and CPTforms a more complete picture of the doping-dependentdimensional crossover. V. CONCLUSIONS
This study has systematically examined the magneticproperties of the Hubbard model in one and two dimen-sions and explored the dimensional crossover in the half-filled and doped systems.In one dimension, the spin dynamics of a strongly-correlated Hubbard chain has been explored for differ-ent doping levels. To understand the role of short-rangecorrelations in the spin dynamics, the doping evolutionof the NN spin correlation function has also been exam-ined. In contrast to the two-dimensional system, dopedholes appear to interact weakly with the magnetism atlow doping levels, but have an enhanced impact at highdoping levels. A simple picture of spin-charge separationcannot explain the trend for all doping levels.In two dimensions, the doping evolution of the spindynamics has been calculated and compared to an ear-lier study in which NNN hopping t (cid:48) was included. Atthe AF wavevector, the spin excitations almost disap-pear upon doping, as ( π, π ) AF order is destroyed. Alongother directions, the spectrum hardens with a persistentintensity. The NN spin correlation function reveals thatdelocalization of doped holes destroys magnetic order in amore subtle way than predicted by a simple local pictureof AF exchange. At lower doping levels, magnetic orderis suppressed below what is expected, while at higherdoping levels, the local picture appears to apply.The crossover of magnetism from one to two dimen-sions provides a means of elucidating the processes in-volved in spin dynamics. When the dimensionalityis tuned between these limits, the magnetic propertiesdrastically change regardless of the doping level. Cru-cially, this suggests that the spin excitations in the two-3dimensional doped Hubbard model cannot be explainedusing the spinon picture known from one dimension.More precisely, we note that:Firstly, doping modifies the t ⊥ dependence of thecrossover at low, but not high, energies. Indeed, dopingenhances the spin dispersion perpendicular to the chains,which can be understood with the three-site term of the t − J model and demonstrates its importance when study-ing doped systems. Comparing the evolution of spin dy-namics at half-filling to that at 40% hole doping demon-strates that persistent coherent spin excitations developat intermediate t ⊥ and smoothly interpolate to the two-dimensional case. Moreover, the sensitivity of spin ex-citations to interchain coupling is momentum dependent.Secondly, dimensionality also affects the way in whichpersistent spin excitations harden upon doping. Indeed,the hardening is enhanced with decreasing t ⊥ . The localpicture used in two dimensions to explain this hardeningis adapted to the anisotropic case. However, for smallinterchain coupling, it predicts a softening of the spinexcitations upon doping, and for large interchain cou-pling, it predicts that the hardening should increase with t ⊥ . Evidently, the local picture does not fully accountfor the hardening mechanism; it may be confounded bythermal effects. Future studies could simulate the small t ⊥ part of the crossover at lower temperatures in order tosearch for this transition from a softening to a hardening.Finally, the dimensional crossover of the NN spin cor-relation sheds light on the role of short-range correla-tions. Along the chain, it smoothly evolves from theone-dimensional result toward the two-dimensional oneas t ⊥ increases, highlighting the very different nature of the interplay between doping and magnetic order in oneand two dimensions. Perpendicular to the chain, the NNcorrelation function exhibits a plateau over a large rangeof doping levels, suggesting that processes may exist tocompensate the local reduction of spin density by thedoped holes.In this work, a coupling-driven dimensional crossoverapproach has been used to calculate the evolution of thespin excitation spectrum for all momentum points. Analternative strategy would be to study a geometric di-mensional crossover tuned by increasing the number oflegs in a ladder. This future work would provide anadditional perspective on the interplay between dimen-sionality, doping, and magnetism in the Hubbard model.This research was supported by the U.S. Departmentof Energy (DOE), Office of Basic Energy Sciences, Divi-sion of Materials Sciences and Engineering, under Con-tract No. DE-AC02-76SF00515. Y.F.K. was supportedby the Department of Defense (DOD) through the Na-tional Defense Science and Engineering Graduate Fel-lowship (NDSEG) Program and by the National Sci-ence Foundation (NSF) Graduate Research Fellowshipunder Grant No. 1147470. K.W. acknowledges sup-port by Narodowe Centrum Nauki (NCN, National Sci-ence Centre) under Projects No. 2012/04/A/ST3/00331and No. 2016/22/E/ST3/00560.S.J. acknowledges par-tial support from The Joint Directed Research and De-velopment (JDRD) program with Oak Ridge NationalLaboratory. The computational work was partially per-formed at the National Energy Research Scientific Com-puting Center (NERSC), supported by the U.S. DOEunder Contract No. DE-AC02-05CH11231. S. Maekawa, T. Tohyama, S. E. Barnes, S. Ishihara,W. Koshibae, and G. Khaliullin,
Physics of transitionmetal oxides (Springer, 2004). Y.-J. Kim, J. P. Hill, H. Benthien, F. H. L. Essler, E. Jeck-elmann, H. S. Choi, T. W. Noh, N. Motoyama, K. M. Ko-jima, S. Uchida, et al., Phys. Rev. Lett. , 137402 (2004). T. Valla, P. D. Johnson, Z. Yusof, B. Wells, Q. Li, S. M.Loureiro, M. M. R. J. Cava, Y. Mori, M. Yoshimura, andT. Sasaki, Nature , 627 (2002). A. Zheludev, M. Kenzelmann, S. Raymond, E. Ressouche,T. Masuda, K. Kakurai, S. Maslov, I. Tsukada, K. Uchi-nokura, and A. Wildes, Phys. Rev. Lett. , 4799 (2000). B. J. Kim, H. Koh, E. Rotenberg, S.-J. Oh, H. Eisaki,N. Motoyama, S. Uchida, T. Tohyama, S. Maekawa, Z.-X.Shen, et al., Nat. Phys. , 397 (2006). T. Giamarchi,
Quantum Physics in One Dimension (Ox-ford University Press, 2003). J. Voit, Rep. Prog. Phys. , 977 (1995). J. Solyom,
Fundamentals of the Physics of Solids, Nor-mal, Broken-Symmetry, and Correlated Systems , vol. 3(Springer, 2003). B. Sciolla, A. Tokuno, S. Uchino, P. Barmettler, T. Gia-marchi, and C. Kollath, Phys. Rev. A , 063629 (2013). Z. N. C. Ha and F. D. M. Haldane, Phys. Rev. Lett. , 2887 (1994), URL https://link.aps.org/doi/10.1103/PhysRevLett.73.2887 . R. Preuss, A. Muramatsu, W. von der Linden, P. Dieterich,F. Assaad, and W. Hanke, Phys. Rev. Lett. , 732 (1994). F. Forte, M. Cuoco, C. Noce, and J. v. d. Brink, Phys.Rev. B , 245133 (2011). R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Phys.Rev. D , 2278 (1981). J. E. Hirsch, Phys. Rev. B , 4403 (1985). F. Assaad, Helvetica Physica Acta , 580 (1990). S. R. White, D. J. Scalapino, R. L. Sugar, E. Y. Loh, J. E.Gubernatis, and R. T. Scalettar, Phys. Rev. B , 506(1989). C. N. Varney, C.-R. Lee, Z. J. Bai, S. Chiesa, M. Jarrell,and R. T. Scalettar, Phys. Rev. B , 075116 (2009). E. Dagotto, Rev. Mod. Phys. , 763 (1994). R. B. Lehoucq, D. C. Sorensen, and C. Yang,
ARPACKUsers’ Guide: Solution of Large-Scale Eigenvalue Prob-lems with Implicitly Restarted Arnoldi Methods (SIAM,Philadelphia, 1998). D. S´en´echal, D. Perez, and M. Pioro-Ladri`ere, Phys. Rev.Lett. , 522 (2000). D. S´en´echal, D. Perez, and D. Plouffe, Phys. Rev. B ,075129 (2002). A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg,Rev. Mod. Phys. , 13 (1996). D. I. Khomskii,
Basic aspects of the Quantum theory ofsolids, Order and Elementary Excitations (Cambridge Uni-versity Press, 2010). C. J. Jia, E. A. Nowadnick, K. Wohlfeld, Y. F. Kung,C. C. Chen, S. Johnston, T. Tohyama, B. Moritz, andT. P. Devereaux, Nat. Comm. , 3314 (2014). E. Dagotto, A. Moreo, F. Ortolani, D. Poilblanc, and J. Ri-era, Phys. Rev. B , 10741 (1992). M. Le Tacon, G. Ghiringhelli, J. Chaloupka, M. M. Sala,V. Hinkov, M. Haverkort, M. Minola, M. Bakr, K. Zhou,S. Blanco-Canosa, et al., Nature Physics , 725 (2011). M. Le Tacon, M. Minola, D. Peets, M. M. Sala, S. Blanco-Canosa, V. Hinkov, R. Liang, D. Bonn, W. Hardy, C. Lin,et al., Phys. Rev. B , 020501 (2013). M. Dean, A. James, R. Springell, X. Liu, C. Monney,K. Zhou, R. Konik, J. Wen, Z. Xu, G. Gu, et al., Phys.Rev. Lett. , 147001 (2013). M. Dean, G. Dellea, R. Springell, F. Yakhou-Harris,K. Kummer, N. Brookes, X. Liu, Y. Sun, J. Strle,T. Schmitt, et al., Nature Materials , 1019 (2013). W. Lee, J. Lee, E. Nowadnick, S. Gerber, W. Tabis,S. Huang, V. Strocov, E. Motoyama, G. Yu, B. Moritz,et al., Nature Physics , 883 (2014). Y. F. Kung, E. A. Nowadnick, C. J. Jia, S. Johnston,B. Moritz, R. T. Scalettar, and T. P. Devereaux, Phys.Rev. B , 195108 (2015). M. Raczkowski and F. F. Assaad, Phys. Rev. B , 085120(2013). S. Biermann, A. Georges, A. Lichtenstein, and T. Gia-marchi, Phys. Rev. Lett. , 276405 (2001). M. Raczkowski, F. F. Assaad, and L. Pollet, Phys. Rev. B , 045137 (2015). B. Lenz, S. R. Manmana, T. Pruschke, F. F. Assaad, andM. Raczkowski, Phys. Rev. Lett. , 086403 (2016). M. Tsuchiizu, Y. Suzumura, and C. Bourbonnais, Phys.Rev. Lett. , 126404 (2007). M. Raczkowski and F. F. Assaad, Phys. Rev. Lett. ,126404 (2012). D. Boies, C. Bourbonnais, and A.-M. S. Tremblay, Phys.Rev. Lett. , 968 (1995). S. Capponi, D. Poilblanc, and F. Mila, Phys. Rev. B ,17547 (1996). M. Fabrizio, A. Parola, and E. Tosatti, Phys. Rev. B ,3159 (1992). H. J. Schulz, Phys. Rev. Lett. , 2831 (1990). D. Poilblanc, H. Endres, F. Mila, M. G. Zacher, S. Cap-poni, and W. Hanke, Phys. Rev. B , 10261 (1996). P. Ribeiro, P. D. Sacramento, and K. Penc, Phys. Rev. B , 045112 (2011). I. Affeck, M. P. Gelfand, and R. R. P. Singh, Journal ofPhysics A: Mathematical and General , 7313 (1994). Y. J. Kim and R. J. Birgeneau, Phys. Rev. B , 6378(2000). A. W. Sandvik, Phys. Rev. Lett. , 3069 (1999). R. dos Santos, Brazilian Journal of Physics , 277 (2003). M. Jarrell and J. Gubernatis, Phys. Rep. , 133195(1996). P. W. Anderson, Phys. Rev. , 2 (1959). J. Hubbard, in
Proceedings of the Royal Society of LondonA: Mathematical, Physical and Engineering Sciences (TheRoyal Society, 1963), vol. 276, pp. 238–257. P. A. Lee, N. Nagaosa, and X.-G. Wen, Rev. Mod. Phys. , 17 (2006). V. I. Iglovikov, E. Khatami, and R. T. Scalettar, Phys.Rev. B , 045110 (2015). M. Iazzi, A. A. Soluyanov, and M. Troyer, Phys. Rev. B , 115102 (2016). J. E. Hirsch and D. J. Scalapino, Phys. Rev. B , 7169(1983). H. Tsunetsugu, J. Phys. Soc. Jap. , 1460 (1991). H. Benthien and E. Jeckelmann, Phys. Rev. B , 205128(2007). M. J. Bhaseen, F. H. L. Essler, and A. Grage, Phys. Rev.B , 020405 (2005). G. Martinez and P. Horsch, Phys. Rev. B ,317 (1991), URL https://link.aps.org/doi/10.1103/PhysRevB.44.317 . J. Bona, P. Prelovek, and I. Sega, EPL (EurophysicsLetters) , 87 (1989), URL http://stacks.iop.org/0295-5075/10/i=1/a=015 . S. Winterfeldt and D. Ihle, Phys. Rev. B ,9402 (1998), URL https://link.aps.org/doi/10.1103/PhysRevB.58.9402 . A. A. Vladimirov, D. Ihle, and N. M. Plakida, Phys. Rev.B , 104425 (2009), URL https://link.aps.org/doi/10.1103/PhysRevB.80.104425 . S. Kar and E. Manousakis, Phys. Rev. B ,205107 (2011), URL https://link.aps.org/doi/10.1103/PhysRevB.84.205107 . J. Bala, A. M. Oles, and J. Zaanen, Phys. Rev. B , 4597(1995). W. Chen and O. P. Sushkov, Phys. Rev. B ,184501 (2013), URL http://link.aps.org/doi/10.1103/PhysRevB.88.184501 . M. Daghofer, K. Wohlfeld, A. M. Oles, E. Arrigoni, andP. Horsch, Phys. Rev. Lett. , 066403 (2008). L. Liu, H. Yao, E. Berg, S. R. White, R. Steven, and S. A.Kivelson, Phys. Rev. Lett.108