Magnetic order in the repulsive Fermi-Hubbard model in three-dimensions and the crossover to two-dimensions

2.6 MB·English

by  Jie Xu

Magnetic order in the repulsive Fermi-Hubbard model in three-dimensions and the
crossover to two-dimensions
Jie Xu, Simone Chiesa, Eric J. Walter, Shiwei Zhang
Department of Physics, College of William and Mary, Williamsburg, VA 23187, USA
Systemsoffermionsdescribedbythethree-dimensional(3D)repulsiveHubbardmodelonacubic
lattice have recently attracted considerable attention due to their possible experimental realiza-
tion via cold atoms in an optical lattice. Because analytical and numerical results are limited
away from half-filling, we study the ground state of the doped system from weak to intermediate
interaction strengths within the generalized Hartree-Fock approximation. The exact solution to
the self-consistent-field equations in the thermodynamic limit is obtained and the ground state is
shown to exhibit antiferromagnetic order and incommensurate spin-density waves (SDW). At low
3 interaction strengths, the SDW has unidirectional character with a leading wave-vector along the
1 (cid:104)100(cid:105)-direction, and the system is metallic. As the interaction increases, the system undergoes a
0 simultaneous structural and metal-to-insulator transition to a unidirectional SDW state along the
2 (cid:104)111(cid:105)-direction,withadifferentwavelength. Wesystematicallydeterminethereal-andmomentum-
space properties of these states. The crossover from 3D to two-dimensions (2D) is then studied
n
by varying the inter-plane hopping amplitude, which can be experimentally realized by tuning the
a
J distancebetweenastackofsquare-latticelayers. Detailedcomparisonsaremadebetweentheexact
numerical results and predictions from the pairing model, a variational ansatz based on the pair-
5
ing of spins in the vicinity of the Fermi surface. Most of the numerical results can be understood
2
quantitatively from the ansatz, which provides a simple picture for the nature of the SDW states.

]
l PACSnumbers: 75.30.Fv,71.15.Ap,71.45.Lr,71.10.Fd,75.10.Lp,75.50.Ee
e
-
r
t
s I. INTRODUCTION larly important and timely to understand such phases in
. the Hubbard model. Somewhat surprisingly, apart from
t
a half-filling (one particle per site), which displays anti-
m Over the past several years, optical lattices have be-
ferromagnetic (AFM) order, the nature of the magnetic
come an increasingly powerful tool for emulating many
d- systems in condensed matter physics1–4. An optical lat- properties in the 3D Hubbard model has not been char-
acterized, even at the mean-field level. In this work, we
n tice can provide exceptionally clean access to a vari-
study the magnetic properties in the ground states of
o ety of model many-body Hamiltonians in which param-
c eters can be systematically tuned and controlled. Thus, the 3D Hubbard model and in the crossover regime, us-
[ ing generalized Hartree-Fock (HF) theory. It is shown
they make possible quantitative experimental study of
that the system has a tendency to form unidirectional
1 thepropertiesofinteractingelectronmodels, whichhave
v proven extremely challenging for analytic and numerical spin-density wave (SDW) states with AFM order and a
6 approaches alone. The combination of these approaches modulating wave along either the (cid:104)100(cid:105)- (at low U/t)
8 or the 111 -direction (at higher U/t). We examine the
1 presents unprecedented opportunities for improving our evolutio(cid:104)n of(cid:105)the SDW wavelength in the full mean-field
understanding of interacting electron systems, by test-
6
solution as U, density and t vary and characterize the
1. ing theoretical concepts and increasing the accuracy and ground state by its propert⊥ies in real and momentum
predictivepowerofnumericalapproachesviacomparison
0 space.
with experiment.
3
1 The one-band Hubbard model is one of the most fun- Despite the simple nature of the mean-field approach,
: damental models in condensed matter physics. It has thedeterminationofthecorrectequilibriumpropertiesin
v been widely studied in two dimensions (2D)5–20 as the these models is not straightforward20,22. The main chal-
i
X simplest model for the Cu-O plane in cuprate supercon- lenge lies in finding an unbiased strategy to determine
r ductors. Forthethree-dimensional(3D)Hubbardmodel, the leading wave-vector(s) characterizing the spatial de-
a however,considerablylessisknownfromboththeoretical pendence of the order parameter. Calculations are per-
and experimental sides. Optical lattices play, in this re- formed in a real space simulation cell and most choices
spect, aparticularlyfundamentalroleastheyallowfora of the cell will return solutions that are biased by finite-
clean experimental realization of the 3D model and offer size effects. This is often further complicated by shell
theinterestingpossibilityoftuningthehoppingparame- effects and sensitivity of the solution to the topology of
teralongonedirection,t ,therebyallowingasystematic the Fermi surface, which often lead to local minimum
studyoftheevolutionof⊥propertiesasthesystemcrosses solutions.
over from 3D to 2D. To overcome the difficulties, it is necessary to move to
Thankstoadvancesintheabilitytodirectlycoolatoms larger and larger cells and gain insights from the evolu-
in optical lattices21, experiments are nearing the realiza- tion of the corresponding solutions. This line of attack
tion of phases with magnetic order. It is thus particu- hasbecomeincreasinglypossiblebecauseofthedramatic
2
increase in computing power and continuous algorithmic the Hubbard Hamiltonian reads
progress. In the present work, we combine such an ap-
proach with more targeted searches to obtain the global (cid:88) (cid:16) (cid:17)
minimumsolutionoftheself-consistent-field(SCF)equa- H = −t c†rzσcr(cid:48)zσ+c†r(cid:48)zσcrzσ
tionsinthethermodynamiclimit. Furthermore,weshow rr(cid:48) ,z,σ
(cid:104) (cid:105)
(cid:88) (cid:16) (cid:17)
htioownatlhaennsautmzebraicsaeldroensutlhtsecpaanirbineguonfdseprisntsooindtbhyeavivcainriitay- −t⊥ c†rzσcrz(cid:48)σ+c†rz(cid:48)σcrzσ
r, zz(cid:48) ,σ
of the Fermi surface. Detailed comparisons are made (cid:104) (cid:105)
(cid:88)
between the direct numerical solutions and the pairing +U nrz nrz , (1)
↑ ↓
ansatzpredictions. Theexcellentagreementhelpstopro- r,z
vide a simple, predictive picture for the properties of the
SDW states.
wheretheoperatorc (c )creates(annihilates)apar-
†rzσ rzσ
The mean-field approach is often the starting point
ticlewithspinσ(σ = , )atsite(r,z)andn isthecor-
rzσ
in the study of strongly interacting systems such as the ↑ ↓
respondingnumberoperator. Thehoppingamplitudetis
Hubbard model. Although the approximations involved
between nearest neighbor sites within a plane (denoted
can lead to significant errors, mean-field theory often
by rr in the summation), t is the inter-plane hop-
(cid:48)
provides insights into qualitative aspects of the behav- ping(cid:104)am(cid:105)plitudebetweenneares⊥tneighborsitesbelonging
ior of many-body systems. Moreover, comparisons with
to different planes (denoted by zz in the summation),
(cid:48)
quantum Monte Carlo results19 have shown20 that, in (cid:104) (cid:105)
and U > 0 is the on-site interacting strength. Through-
2D, the mean-field solution captures the basic physics of
out this work, energy is quoted in units of t and we set
SDW states at intermediate interaction strengths, and
t=1. TheHamiltonianinEq.(1)describesthe3Dcubic
providesagoodqualitative(orevenquantitativeinsome
Hubbard model when t =1, the crossover between the
aspects) description of the magnetic correlations in the squareandcubiclattice⊥swhen0<t <1andastackof
true ground state. Because it is reasonable to expect a decoupled 2D Hubbard planes when⊥t =0. Only unpo-
similar level of accuracy for the models presently consid- larizedsystemsareconsidered, i.e., th⊥eaveragedensities
ered,weexpectthatourfindingswillprovideguidanceto
of the two spin species are kept equal: n = n . The
many-body approaches and experimental studies alike. nature of the ground state is thus character↑ized b↓y three
WehavelimitedourstudytoU (cid:46)6t,wherethemean- parameters: the inter-plane hopping amplitude t , the
fieldapproachcanbeexpectedtoofferusefulinsight. Be- on-site repulsion U and the doping (hole density)⊥
low we will discuss the mean-field predictions and their
implications (and the caveats) on the true many-body
states drawing from the comparison in 2D19,20. Clearly, h 1 (n +n ). (2)
≡ − ↑ ↓
the form of generalized mean-field theory considered in
thisworkwillnotcaptureexoticinstabilities,suchasun- Theparticle-holetransformation,c ( 1)x+y+zc ,
conventionalpairingorder. Indeed,asU increases,itwill maps the h < 0 sector into the h >†rzσ0→one−, regardlesrszσof
become increasingly inadequate for magnetic properties the value of t or U, and we therefore confine our study
as well. to h>0. ⊥
The remainder of the paper is organized as follows. In
Athalf-filling,h=0,thenon-interactingFermisurface
Sec.II,weintroducetheHamiltonian,andbrieflyoutline
is given by 2(cosk +cosk +t cosk ) = 0. Despite
some of its basic properties to facilitate the ensuing dis- thelackofs−ymmetryxbetweenythe⊥z-andzther-directions
cussion. In Sec. III, we summarize the strategies used to
for any t = 1, perfect nesting via Q ( π, π, π)
solvethemean-fieldequations. Resultsforthe3Dmodel remains t⊥hr(cid:54)oughout the whole surface,≡an±d ca±uses±an
are presented in Sec. IV; numerical results for the 100 -
(cid:104) (cid:105) AFM instability for any t = 0 and arbitrary small U
and the (cid:104)111(cid:105)-SDW are followed by a discussion where values. The evolution of ⊥the(cid:54) non-interacting half-filled
the pairing ansatz is first introduced and then used to
Fermi surface as t varies is shown in Fig. 1. In the first
helpunderstandthenumericalfindings. Thedimensional column, representi⊥ng the 2D limit, the Fermi surface has
crossover results are then presented in Sec. V, followed
no dependence on k and any wave-vector of the form
again by a discussion based on insights from the pairing z
( π, π,q)isperfectlynestedonit. Thearbitrarinessof
ansatz. We conclude in Sec. VI. ± ±
q is reflected in the complete lack of correlation between
the r-planes. The large nesting degeneracy is abruptly
interrupted as soon as t = 0, and Q remains the only
II. BACKGROUND nesting vector as the sy⊥ste(cid:54)m evolves from the 2D limit
toward 3D. The middle and bottom rows illustrate pro-
Given our goal to study both the 3D case and the jections of the Fermi surface along the 100 - and 111 -
(cid:104) (cid:105) (cid:104) (cid:105)
crossover from 3D to 2D, it is most convenient to define directions; as we shall see in Sec.’s IV and V, the pro-
the3DHubbardHamiltonianasastackofsquare-lattice jected area of the Fermi surface plays a central role in
planes. We will use r (x,y) to denote in-plane coor- determining the character of the SDW in the proximity
≡
dinates and z to label the planes. With this convention of h=0.
3
are N N matrices with elements
×
[H (k)] = t (k)+δ (UD µ),
σ βγ βγ βγ βσ
− − (5)
[S±(k)]βγ =UδβγSβ±,
(cid:80)
where tβγ(k)= Lexp(ik·L)tβ,γ+L, and Dβσ, Sβ± and
µaredeterminedbytherequirementthatthefreeenergy
F = H TS is a minimum for the targeted average
0 0
(cid:104) (cid:105) −
density n=n +n . This amounts to the following SCF
(gap) equation↑s ↓
(cid:90)
1
Dβ,−σ =(2π)3 dk(cid:104)c†βσ(k)cβσ(k)(cid:105)0
(cid:90)
1
Sβ± =− (2π)3 dk(cid:104)c†β,±σ(k)cβ,∓σ(k)(cid:105)0 (6)
(cid:90)
1 (cid:88)
n=N(2π)3 dk(cid:104)c†βσ(k)cβσ(k)(cid:105).
β,σ
To locate the ground state we proceed with two com-
plementary approaches. In the first approach we select
the L ’s so that they span a large supercell containing
i
(5000) sites. A twisted boundary condition23,24 is ap-
O
plied, namely, using a single randomly selected k-point
FIG. 1: (Color online) Non-interacting half-filled Fermi sur- in place of the integrals in Eq. (6). The iterative process
face from different view angles: 3D (top), along [010] (mid- is started with various initial states, including random
dle),along[11¯1](bottom). Fromlefttoright,thecolumnsare
ones, and multiple annealing cycles are performed. In
fort =0,0.5and1,respectively. Notethatonlythesurface
⊥ each cycle a random perturbation (whose strength can
inoneoctantisshowninthebottomrow. Perfectnestingvia
be controlled) is applied to a converged solution and the
Q=(±π,±π,±π) holds for any t at half-filling.
⊥ self-consistent process is repeated. Separate calculations
for different k-points are done to check for consistency.
Once an understanding of the character of the ground
III. METHOD
state is gained, we use a second approach to target the
specific family of states compatible with the results of
The following mean-field formalism in real space is the random search. For instance, suppose the random
used in this work. A simulation cell of N sites is defined searchfindsaunidirectionalSDWatsmallU valueswith
by three vectors, L1, L2 and L3, whose components are wave-vector along the 100 -direction. We then choose a
(cid:104) (cid:105)
integers. Bloch states are then introduced as cluster of L = (L,0,0), L = (0,1,0) and L = (0,0,1)
1 2 3
with G=π(( 1)L+2l,1,1), where l is the number of os-
(cid:88) (cid:2) (cid:3) −
c (k) c exp ik L , (3) cillationsoftheorderparameter, chosentobeaninteger
β β+L
∝ · orhalfaninteger. Foragivensetofthethreeparameters
L
(t , U, h), Lisfinelyscanned(withLontheorderof50
whereLarevectorsoftheformL=n L +n L +n L , an⊥d step size of 1) until the energy minimum is found.
1 1 2 2 3 3
k is a reciprocal lattice vector that is free to vary within A large number of k-points is used (on the order of 100
the first Brillouin zone (BZ) defined by the L ’s, and β in the two short directions and a few in the other) so
i
labels sites inside the simulation cell. Using these states, that the character and properties of the targeted states
the mean-field Hamiltonian can be decoupled into a sum can be accurately determined. This approach allows us
(cid:80)
of k-dependent pieces, H = H (k), with each piece to study different forms of SDW and long wavelength
0 k 0
of the form modes without increasing the computational cost.
In our study, we mix the two numerical approaches
H0(k)=[c†↑c†↓](cid:20)HS↑(+k) H↓(kS−−G)(cid:21)[c↑c↓]T, (4) aexsanmepedlee,dcoamndpaursiseonthoefmeninergcioemspalmemonegntsaervyerwaalyfas.miFlioesr
of SDW is made with the second approach. To confirm
where c (c ) represents a row of operators c (k) thecorrectnessofthegroundstate,thesolutionsarethen
β
(c (k ↑G))w↓ithindexβ runningthroughtheN site↑sof checkedagainstdifferentinitialstatesandannealingpro-
β
the↓cel−l. A non-zero value of G causes the spin densities cedures using the first approach on a supercell commen-
at β and β +L to be related via a rotation by G L surate with the optimal wavelength.
i i
·
around the z-axis. Charge and spin densities along z- Various observables are computed to characterize the
direction obey periodic boundary conditions. H and S converged solutions. The local charge density ρ and the
±
4

16 16 16
14
12 12 12
z 8 z 8 z 8 z 7
4 4 4
16 12y8 4 4 8x 12 16 16 12y8 4 4 8x 12 16 16 12y8 4 4 8x 12 16 16 12y8 4 4 8x 12 16
0.896 0.898 0.9 0.902 −0.1 0 0.1 −0.2 0 0.2 −0.2 0 0.2
16 16
FIG. 3: (Color online) Spatial dependence of the order
12 12 parameter m for h = 13/128, same as in Fig. 2, but with
U =2.9,ona16×16×16supercell(left)anda16×16×14
z 8 z 8 supercell (right). Uniform AFM order in the xy-plane disap-
pears for L = 16, but linear SDW along z-direction is seen
z
4 4 again on the right with L =14.
z
0.894 0.897 0.9 0.903 −0.12 −0.06 0 0.06 0.12
ρ m with c (cid:80) exp( ik R)c . We use n to compare
kσ ∝ R − · Rσ k
theconvergedmean-fieldsolutionwiththepairingansatz
FIG.2: (Coloronline)Chargedensityρ(left)andorderpa- predictionandAktocharacterizetheFermisurfaceofthe
rameterm(right)ofthesolutionforthe3DHubbardmodel. ordered phase.
Shownisa16×16×16supercell,withh=13/128atU =2.5.
Alinearwaveisseenalongthez-direction,withuniformAFM
order in the xy-plane. The bottom panel shows a line cut IV. 3D RESULTS
along z-direction, with dots the actual data and the line a
sinusoidal fit.
A. SDW correlation in the (cid:104)100(cid:105)-direction
local order parameter, identified as the local staggered At half-filling, the existence of perfect nesting allows
magnetization m, are defined as an AFM solution for any U >0. Away from half-filling,
perfect nesting ceases to exist and a finite critical value
ρ(R) n + n , (7)
rz rz of the interaction is needed to cause the onset of order.
≡ (cid:104) ↑(cid:105) (cid:104) ↓(cid:105)
m(R) ( 1)x+y+z( n n ), (8) The critical value U depends on h. Using the first ap-
rz rz c
≡ − (cid:104) ↑(cid:105)−(cid:104) ↓(cid:105)
proach described in Sec. III, we have determined that,
and used to characterize the state in real space (here justaboveU , thegroundstateofthesystemisanSDW
c
R (r,z)). Since all the minimum energy solutions withmodulationalongthe 100 -direction. Figure2illus-
we≡find are unidirectional spin/charge density waves tratesthespatialdependen(cid:104)ceof(cid:105)ρandmina16 16 16
(SDW/CDW),itisnaturaltocharacterizethembytheir supercell at h = 13/128 0.10 and U = 2.5. T×he S×DW
modulation wavelength along a relevant Cartesian axis, is characterized by a sin(cid:39)gle wave-vector and λ = 8.
100
λSDW/CDW, defined by the leading component of the TheamplitudeoftheSDWis 0.1,roughlythi(cid:104)rty(cid:105)times
Fourier transform of ρ and m, respectively. The min- larger than that of the CDW(cid:39). The simple form of the
ima in the CDW are found to coincide with nodes in order found for m(R) is indicative of the proximity of U
the SDW, thus 2λCDW = λSDW. Below we will some- to Uc. All of the observations above are consistent with
times discuss our results in terms of a single wavelength the pairing model, as discussed below in Sec. IVC.
λ λCDW, which can also be identified as the distance We next examine the evolution of λ as the inter-
≡ 100
between two consecutive nodes of the order parameter. action strength changes. Keeping h an(cid:104)d t(cid:105)he simulation
Whenwerefertothedirectionofthemodulation,wewill cellunchangedandincreasingU from2.5to2.9,theran-
use 100 todenotesymmetry-equivalent[100]-directions, dom search returns the state displayed on the left panel
(cid:104) (cid:105)
and similarly for [110] and [111]. of Fig. 3, which suggests that a more complicated type
To characterize the system in momentum space, we of order is seemingly settling in. We apply our second
use the momentum distribution nk and the momentum- approach,usingadensek-pointgridanda1 1 Lsim-
resolved single-particle spectral function Ak(ω), defined ulationcell,tosearchfortheoptimalwavelen×gth.×Weuse
as large L (containing about 8 nodes in the cell) and vary
its value until an energy minimum is found. Figure 4
nkσ =(cid:104)c†kσckσ(cid:105), (9) shows the result of such minimization in terms of λ 100 ;
Akσ(ω)= π1Im(cid:104)ckσ(ω−H +E0−iη)−1c†kσ(cid:105) (10) theminimumoccurswhenλ(cid:104)100(cid:105) =7,indicatingtha(cid:104)tth(cid:105)e
5
−1.403
1 1 L h=1/32
16× 1×6 16
−1.4031 16×16×14 0.85 h=1/16
× × h=2/25
h=3/32
h=13/128
−1.4032
0.8 h=3/25
e h=1/8
t i
si −1.4033 00 h=3/20
E/ α1h0.75 h=1/5
−1.4034
0.7
−1.4035
−1.4036 0.65
6.4 6.6 6.8 7 7.2 7.4 7.6 7.8 8 8.2 2 2.5 3 3.5 4 4.5 5
λ 100 U
h i
FIG.4: (Coloronline)Energyof(cid:104)100(cid:105)-SDW(blue)vs. λ FIG. 5: (Color online) Characteristic wavelength as a func-
(cid:104)100(cid:105)
for a system of h=13/128 and U =2.9. Horizontal lines are tionofU atvariousdoping. α givesthemodulationwave-
(cid:104)100(cid:105)
theenergiesofthecalculationsshowninFig.3. Theminimum length, λ , in units of 1/h. As U is increased, the value
(cid:104)100(cid:105)
of(cid:104)100(cid:105)-SDWisreachedwhenλ =7. Thestateintheleft of α converges to approximately 2/3 at small h (slightly
(cid:104)100(cid:105) (cid:104)100(cid:105)
panel of Fig. 3 leads to an energy higher than the minimum larger at larger h).
but lower than the energy with λ =8.
(cid:104)100(cid:105)
in Fig. 5, with α defined as
100
(cid:104) (cid:105)
16 16 16supercellisnotcommensuratewiththewave- α =hλ . (11)
len×gth×of the minimum energy SDW state and that the (cid:104)100(cid:105) (cid:104)100(cid:105)
patternintheleftpanelofFig.3isaresultoffrustration Whenthedopingissmall,thewavelengthofthemodula-
from an incommensurate supercell size. We next return tion is proportional to 1/h, with α almost indepen-
100
to our first approach, and perform a new mean-field cal- dent of U and roughly equal to 2/3(cid:104). Fo(cid:105)r larger h, α
100
culation, with random initial guess and annealing, on a converges to a slightly larger value. There is a gen(cid:104)era(cid:105)l
16 16 14supercell,asizewhichiscommensuratewith trend of an increase of α as U approaches U from
100 c
the×wav×elengthoftheminimumenergysolution. Andin- above. We will be able to(cid:104)rat(cid:105)ionalize these trends within
deedwefindthepredictedstatewithλ =7correctly the pairing model in Sec. IVC.
100
reproduced (right panel in Fig. 3). (cid:104) (cid:105) The evolution of the properties of the 3D SDW state
with interaction U is similar to what is observed in 2D.
In Fig. 4, we report the energies of the two large su-
Figure 6 shows 1D cuts of m and ρ in the z-direction, at
percell calculations of Fig. 3, to verify that the energy
h=0.05 and U values such that α has saturated to
odbutcaeidnebdyinthtehe1161×16L×c1lu4stseurpesercaerlclhiswciothrreλctly re=pro7-. ∼2/3. The figure illustrates the ex(cid:104)i1s0t0e(cid:105)nce of Uc and the
The energy of th×e 16× 16 16 calculation, on(cid:104)t1h00e(cid:105)other increase of the SDW and CDW amplitudes with U. It
× × also shows the crossover from a regime where the order
hand, falls between those of λ =7 and 8. This gives
100 parameter has a smooth sinusoidal modulation and the
a clear illustration of the chara(cid:104)cte(cid:105)ristics of the two types
holes are delocalized, to one where it is characterized by
of approaches. The supercell being incommensurate pre-
domainwallsorstripes, withholeslocalizedinthenodal
ventsthesolutionfromcollapsingontothelowestenergy
regions.
SDW state of λ = 7. The self-consistent solution in
100 There exists an important difference in the physics of
a large supercel(cid:104)l th(cid:105)en finds a different pattern that cor-
the mean-field ground state of the 3D system and its 2D
responds to the true ground state compatible with the
counterpart. While, inthelatter, thesystemremainsin-
imposed constraint. The energy of this state, computed
sulating when lightly doped, the 3D model immediately
by fixing the density and converging the value using a
turns metallic. The difference is a consequence of the
dense k-point mesh, is higher than the global minimum,
differentbehaviorsofthemodulatingwavelength. In2D,
but lower than that of the SDW state with an imposed
α is unity, independent of h and U, while in 3D it varies
wavelength λ =8.
100 with parameters and has a non-integer value. To illus-
(cid:104) (cid:105)
We proceed to determine the exact dependence of the trate this in a simple case, consider a value of doping h
wavelengths on h and U by explicit solutions of the SCF such that λ ( α/h) is an integer. For such a system to
≡
equationsin1 1 Lclusters. Verificationsoftheresults beaninsulator,thenumberofparticlesina1 1 λcell,
× × × ×
are done on large supercells whose sizes are commensu- (1 h)λ=λ α,willhavetobeaninteger. However,be-
− −
rate with the wavelengths. Our results are summarized cause α 2/3 in the limit of small doping, the condition
∼
6

1

0.96
22
0.92
ρ0.88 16
0.84
12
0.25
U =2.0
0.2 U =2.5 z11 z 8
0.15 UU ==33..05 4
0.1
0.05 16
12 16
m 0 16 8 12
−−0.00.51 1 2y8 4 4 8x 12 16 y 4 4 8x
−0.15
−0.2 −0.5 0 0.5 −0.5 0 0.5
−0.25

1 14 27 FIG. 7: (Color online) Order parameter m at h = 1/8 and
z
U =5.0 in a 16×16×22 supercell (left) and a 16×16×16
supercell (right). Though the supercell of 16×16×22 is
FIG.6: (Coloronline)Chargedensityρ(top)andorderpa- commensurate with the optimal (cid:104)100(cid:105)-wavelength, the ran-
rameterm(bottom)vs.U. Thesystemhasdopingh=0.05. dom search produces a lower energy solution. The minimum
Eachcurveisa1Dcutalongz-direction, thedirectionofthe energy solution is an SDW along the (cid:104)111(cid:105)-direction, which
modulation. Beyond Uc, the (cid:104)100(cid:105)-SDW/CDW amplitudes is correctly reproduced by a random search in a 16×16×16
increase with U and the solution evolves from a sinusoidal supercell as shown on the right.
wave to domain walls. The CDW amplitude is much weaker
than that of the SDW.
U = 4.0 U = 4.5
−1.231 −1.158
cannot be satisfied, and the system is necessarily metal- e −1.232 −1.159
t
lic. A related way to see this is to consider the case of si −1.233 −1.16
/
domainwallstates,forexampleU =3.5inFig.6. Inside E −1.234
−1.161
eachdomainwall(nodalregion)arelocalizedholeswhose
−1.235
integrated (along the direction of the modulation) den-
5 6 7 8 9 5 6 7 8
sity is α. Thus the domain wall as a whole will act as a
U = 5.0 U = 5.5
quasi-2Dliquidofholeswithnon-integerdensity. Wewill −1.034
discuss the corresponding momentum space signature in −1.094
the next section. te −1.036
i
s
/ −1.096 −1.038 100
E h110i
B. SDW modulation along the (cid:104)111(cid:105)-direction −1.04 h111i
−1.098 h i
5 6 7 8 5 6 7 8
AsshownaboveandfurtherdiscussedinSec.IVC,an λ λ
orientation of the SDW other than 100 is not the so-
lution in the proximity of U . Howev(cid:104)er, w(cid:105) hen the inter- FIG. 8: (Color online) Energies per site of SDW in (cid:104)100(cid:105)-,
c
(cid:104)110(cid:105)- and (cid:104)111(cid:105)-directions vs. λ for h=1/8. At and above
action grows larger, other Fermi liquid instabilities be-
U =4.5,thelowestenergystateismodulatedalongthe(cid:104)111(cid:105)-
come possible. This fact is clearly displayed when cal-
instead of (cid:104)100(cid:105)-direction.
culations on supercells commensurate with the optimal
100 -wavelength for a given U do not yield a state with
(cid:104) (cid:105)
100 -SDWorder. Figure7showstheoccurrenceofsuch
(cid:104)a cas(cid:105)e in a calculation with h = 1/8 and U = 5.0, for lowest energy state is given by a 111 -SDW, instead of
(cid:104) (cid:105)
which λ = 5.5. The 16 16 22 supercell should the 100 -order atlowerU. For U =5, theminimum en-
have pre(cid:104)c1i0s0e(cid:105)ly accommodated×4 no×dal planes of the or- ergy(cid:104)solu(cid:105)tioniscorrectlyreproducedbyarandomsearch
der parameter, but rather than doing so, the random in a 16 16 16 supercell as shown in the right panel of
× ×
search produces the lower energy solution shown in the Fig. 7. In our searches, 110 -direction SDW’s are never
(cid:104) (cid:105)
left panel of the figure. found to be the global ground state.
To search for the solution at higher U, we investigate By repeating the same procedure, we construct the
unidirectionalSDW’swithmodulationlyingalongeither equation of states for U = 3 and U = 4 contained in
the 110 - or 111 -direction. An example is given in Fig. 9. At U = 3, there is no density regime where the
(cid:104) (cid:105) (cid:104) (cid:105)
Fig. 8. The energies from constrained searches using the 111 -SDW is the global ground state. In contrast, for
(cid:104) (cid:105)
second approach are shown as a function of λ for a scan U = 4, a discontinuous transition from 100 to 111
(cid:104) (cid:105) (cid:104) (cid:105)
of U values. It is seen that, at and above U = 4.5, the occurs around n = 0.9 with a small coexistence region.
7
5 1.0 themodelhelpstoexplainthenumericalresultsandpro-
0.02 vides a simple conceptual framework that captures the
4 [111]
0.01 [100] essential physics of the SDW states in 3D and in the
3 0 [110] 0.5 crossover regime discussed in the next section. We first
-3 t ]10 2 0.99 0.995 1 soufm1m00ar-iSzeDWtheaftormmoadleissmt U, a,nfodlltohweendabpypltyheit 1to11th-SeDcaWse.
y [ 0.0 A(cid:104)t lo(cid:105)w U and small h, the pairing model i(cid:104)s defi(cid:105)ned by
g 1
ner spin-orbitals of the form
e
0
-0.5
-1 φ†kσ =ukc†kσ+σvkc†k+qkσ. (12)
UUU === 333 UU == 44
The construction requires excitations to outside the
-2 -1.0
Fermi sea. It is reminiscent of the ansatz used to con-
0.88 0.92 0.96 0.88 0.92 0.96 1.00
struct the BCS pairing states for attractive interactions,
density density
except that here the tendency for small excitations is
FIG. 9: (Color online) Ground-state energy per site from perhaps more “natural”, because of the repulsive inter-
constrainedsearchof(cid:104)100(cid:105)-,(cid:104)110(cid:105)-and(cid:104)111(cid:105)-SDWatU =3 action. A collection of spin pairs in orbitals given by
(left) and U = 4 (right). A linear common shift has been Eq. (12) leads to a uniform charge density, ρ(R) = n,
applied to the energies to highlight the convexity and the and a spin density of the form
different trends. At U = 3, (cid:104)111(cid:105)-SDW is not the global
groundstate. (Theinsetshowsazoomoftheenergydifference 4 (cid:88)
s(R)= a cos(q R) (13)
between (cid:104)111(cid:105)- and (cid:104)100(cid:105)-SDW states as n→1.) At U =4, N k k·
the ground state is (cid:104)111(cid:105)-SDW for n(cid:38)0.92. k
∈R
with a = u v . The region over which k is summed
k k k
R
will be closely related to the non-interacting Fermi sea,
In both cases the low-doping ground state is character- andpreservingthevolumeof4π3n,butwillingeneralbe
ized by a linear energy-density dispersion. This bears slightlymodifiedfromthevariationaloptimization,aswe
two important consequences. First, contrary to what is further discuss below. To ensure orthogonality amongst
observedusingvariationalstateswithauniformspiralor- the spin-orbitals, q must be such that k+q and
k k
derparameter25,thereisnosignofphaseseparationinto k+qk = k(cid:48)+qk(cid:48). The corresponding potentia(cid:54)∈l eRnergy
a half-filled, AFM region and a hole-rich region. Second, per site(cid:54)is then given by
the effective interaction between domain walls is short-
rangedandtheirpreciselocationinthehole-dilutedlimit U (cid:88)
V =Un2 s2(R). (14)
is therefore irrelevant as long as they stay sufficiently far − 4N
R
apart.
We find α =1 at any density for which the 111 - The potential energy lowering relative to the paramag-
111
SDW is the(cid:104)gro(cid:105)und state. The 111 -SDW state(cid:104)s ar(cid:105)e netic (PM) solution is thus:
(cid:104) (cid:105)
fullygapped,owingtotheintegervalueofα ,incon-
111
trast to the metallic behavior of the 100 s(cid:104)tate(cid:105)s. Upon U (cid:88)
increase of U at a constant h, the st(cid:104)ruct(cid:105)ural transition ∆V = −N2 akak(cid:48)[δ(qk+qk(cid:48))+δ(qk−qk(cid:48))],
k,k(cid:48)
is therefore always accompanied by a metal-to-insulator ∈R
(15)
transition. We have verified, for selected cases, that ran-
dom searches on larger supercells with sizes commensu-
where the Kronecker δ is intended as periodic on the
rate to the optimal wave-vector always return unidirec- reciprocal lattice, i.e., modulo 2π in any direction.
tional SDW’s with the same predicted wavelength and
Equation(15)makesitclearthatthemaximumreduc-
orientation. This provides a strong indication that the
tioninV isachievedbyhavingasmanypairsaspossible
character of intermediate U instabilities remains that of
withq ’swhichareparalleloranti-paralleltoeachother.
k
a unidirectional SDW. Thus, as we increase U at con-
NotingthatthevectorQisperfectlynestedwhenh=0,
stant density, the system is always expected to undergo
let us consider the following explicit construction for :
a discontinuous transition from a 100 - to a 111 -SDW R
displacethehalf-filledFermisurfaceineachoctantofthe
(cid:104) (cid:105) (cid:104) (cid:105)
ground state.
first BZ by ∆q/2, choosing the direction that shrinks
±
the Fermi sea and with a length ∆q such that the en-
closed volume is reduced to 4π3n. This construction is
C. A variational pairing ansatz illustrated in the top panel in Fig. 10 for the case of
the 100 -SDW (discussed next in Sec. IVC1). The sur-
(cid:104) (cid:105)
The pairing model is a variational ansatz that has face of can now anchor spin pairs with one common
R
proved extremely helpful in rationalizing the properties pairing vector, by making u less than 1 in a small layer
k
of SDW’s in the mean-field treatment of the electron immediatelyinsidethesurfaceof ,andcorrespondingly
(cid:112) R
gas22,26 and the 2D Hubbard model20. Similarly here, v = 1 u 2 >0.
k k
−| |
8
√2π 0 √2π estimate of α (Eq. (11))
−π π
1 (cid:90)
∆q 1
2 h100i 0.9 α∆q = 4π2 e e (cid:98)e∆q·dS. (17)
(cid:98)∆q (cid:98)
·
∆q S
100
h i
Q
0.8 where (cid:98)e is a relevant Cartesian unit vector.
z0 0 Different directions of ∆q lead to different reconstruc-
k
tions of the non-interacting doped Fermi surface. The
q
h100i 0.7 kinetic energy cost of the pairing ansatz can therefore
be thought of as the combination of two contributions:
the reconstruction energy due to using rather than
0.6 R
the true non-interacting Fermi sea, and the kinetic en-
π π ergychangeduetomovingparticlesfromk(inside )to
−√2π 0 √2−π R
− k[110] 0.5 k+qk (outside), i.e., the non-zero vk’s. It is easy to see
that,similarto2D20,atsufficientlylargeU thepotential
π 0 π
π− π energy lowering will overtake the kinetic energy increase
0.4 for the constructions discussed here. The correct state
is determined by maximizing the gain in the potential
energy from pairing (larger areas near the Fermi surface
0.3
participating with parallel ∆q) while minimizing the ki-
netic energy cost.
y0 0
k 0.2 Theansatzgivesaclearpicturefortheonsetofthein-
stability. First, by its form, the model captures how the
energetically costly CDW can be suppressed compared
0.1 to the SDW. Second, it indicates that amongst different
possiblechoicesofq ,theoneinvolvingonlyparalleland
k
anti-parallel vectors are optimal. Third, the direction of
π π
− π 0 π− ∆qmustbesuchastoleadtotheminimalpossiblere-
− ±
kx construction of the doped Fermi surface. These findings
are in qualitative agreement with the numerical results
FIG. 10: (Color online) Illustration of the pairing model for obtained by the solution of the SCF equations: 1) the
the(cid:104)100(cid:105)-SDWstatein3D.Theschematicdiagramisdrawn SDW is much stronger than the accompanying CDW,
on the actual momentum distribution from the exact numer- 2) the Fermi surface reconstructs in a way to enhance
ical solution for h = 1/8 and U = 4.0. The top panel shows
pairing and the SDW order tends to be unidirectional as
thepairingconstructiononthecontourplotofthe(1¯10)cut.
a result, 3) more drastic reconstructions are only possi-
Thewhitedashedlinesrepresentthehalf-filledsurface,across
ble with larger U. Much quantitative information can
which is the nesting vector Q. The reconstructed surface,
be obtained with straightforward calculations using this
shown as magenta solid lines, is obtained by displacing the
half-filled one along the z-direction by a distance ∆q /2, model, as we discuss next for the 100 - and 111 -SDW
(cid:104)100(cid:105) (cid:104) (cid:105) (cid:104) (cid:105)
as given by Eq. (16). The nesting vector across the shifted states, respectively.
surface,q ,isshownbythelongsolidlinewitharrow. The
(cid:104)100(cid:105)
bottom panel shows n(k) of the k = π plane. This is in a
z
regionwheretheFermisurfacesurvivestheonsetoforderand 1. Analysis of the (cid:104)100(cid:105)-SDW
where it differs more severely from the pairing construction.
The actual Fermi surface, seen distinctly inside the recon-
Among all directions, only a ∆q along the 100 -
structed surface, differs little from the non-interacting Fermi (cid:104) (cid:105)
direction causes the Fermi surface in each octant to
surface(blacksolidline). n dropssharply,andnopairingis
k
shrink equally and this, it can be shown, leads to the
present in this region.
minimum kinetic energy increase at small h. An SDW
with 100 modulation is thus the lowest energy solution
(cid:104) (cid:105)
For small h, we can determine ∆q directly from the at low h and just above Uc, consistent with the results
construction: fromexplicitsolutionsoftheSCFequationsinSec.IVA.
(cid:90) We numerically calculate the projected area along the
Ω
∆q e dS=h BZ, (16) 100 -direction(shownintherightmiddlepanelinFig.1)
(cid:98)∆q· 2 (cid:104)and(cid:105)obtain α 0.63, in very good agreement with
100
S the exact resu(cid:104)lts(cid:105)fr(cid:39)om direct solutions shown in Fig. 5,
where S is the half-filled surface, e is the direction whereα 0.66. Thattheestimatedvalueisslightly
(cid:98)∆q 100
of ∆q, and Ω = (2π)3 is the volume of the first BZ. smalleri(cid:104)sco(cid:105)n(cid:39)sistentwiththepresenceofsurvivingFermi
BZ
Equation(16)impliesalinearrelationshipbetweenhand surface inside the reconstructed doped Fermi surface as
∆q and, using λ = π/(∆q e), provides the following seen in Fig. 10.
(cid:98)
·
9
√2π 0 √2π
−π π
30
π π
2 2
25
z 0 0
k
20
π π
−2 −2
15
FIG. 12: (Color online) (a) Half-filled Fermi surface for the
π π 3D Hubbard model in one of the octants (left) and (b) aver-
−√2π 0 √2−π
− k[110] aged∆q(cid:104)100(cid:105) overdifferentareasvs.doping(right). Theblue
π kz = 0 kz = π/2 kz = π 10 stholeidblluineediont(ibn)(ias)c.aTlchueladteadshaetd/(−doπt/t2ed,−liπn/e2s,iπn/(2b)), sahroewanvears-
aged over varying areas as indicated by the circles with the
samecolor/stylein(a). Theblacklinein(b)showstheaver-
5
y0 agedvalueovertheentiresurface,whichleadstotheestimate
k
of α ∼0.63 discussed in the text.
(cid:104)100(cid:105)
π
− π 0 π π 0 π π 0 π
− k − k − k part of the Fermi surface survives the onset of order; the
x x x
parts that do not survive are in areas around the “hot
FIG. 11: (Color online) Contour plots of the single-particle spots” k ( π/2, π/2, π/2), where the pairing con-
spectralfunctionevaluatedattheFermienergy;onthe(1¯10) struction(cid:39)and±the d±oped F±ermi surface are most similar,
cut (top) and on the kz =0,π/2,π planes (bottom from left andthechangein ∂(cid:15)k isataminimumimplyingacloser
toright)forasystemwithU =2.7andh=1/8. Alargepart proximity to a perf∂ekczt common paring vector.
oftheFermisurfacesurvives,exceptforareasaroundthehot
These understandings allow quantitative explanations
spots where it is most energetically favorable for pairing.
of all the features of the data on α shown in Fig. 5.
100
For example, Fig. 12 shows that t(cid:104)he(cid:105)distance along z-
directionbetweenthedopedandhalf-filledFermisurfaces
A direct comparison between the pairing model and
is at a minimum around the hot spots. Given that such
the exact SCF solution can also be made in momentum
distance equals ∆q/2 in the pairing construction, one
space. We will identify the Fermi surface in the numeri-
finds that the local ∆q value is smaller at the hot spots
calsolutionfrommean-fieldtheoryasthelocusofpoints
thanwhencomputedastheaveragedistanceovertheen-
wheren =0.5. Dependingonthesystemandthevalue
kσ tiresurface. Hence,whenonlythehotspotsareinvolved
ofU,themomentumdistributionoftheexactmean-field
in pairing, a larger α results, as seen at lower values
ground state can maintain a true Fermi surface, char- ofU justaboveU . O(cid:104)b10v0i(cid:105)ously,thepairingreconstruction
acterized by a discontinuity in n , or have it smeared c
k becomesincreasingly accurateas h approaches 0. In this
out by large pairing amplitudes (i.e., u close to 1/√2
k limit, one therefore finds the increasingly smaller U and
neartheboundaryof ). Thetwoscenarioscanoccurin c
R the faster convergence (in U) to the saturated value of
the same system at different k values. The identification
α 2/3 shown in Fig. 5.
using nkσ =0.5 is consistent with both. (cid:104)100(cid:105) ∼
Figure 10 shows, in particular, that the portion of the
Fermi surface where pairing takes place is in very good
2. Application to the (cid:104)111(cid:105)-SDW
agreement with the construction based on the pairing
model,whichindicatesthattheansatzcapturesthedom-
inant ingredient of the physics of the SDW state. The Bytakingadisplacement∆qalongthe 111 -direction,
(cid:104) (cid:105)
figure also provides a direct explanation for the survival we can apply the pairing construction straightforwardly
of the Fermi surface around k = π as it is there that to the diagonal-modulated SDW. Analogously to the
z
thepairingconstructionshowslargediscrepancywiththe 100 case, Fig. 13 shows a remarkable agreement be-
(cid:104) (cid:105)
true Fermi surface. This, in turn, implies that pairing in tween the shifted half-filled surface and the calculated
that region would be associated with too large a kinetic Fermi surface of the SDW. The two cuts in the figure
energycosttobefavorable. Theabsenceofanygapfrom clearly show broken cubic symmetry, where the Fermi
pairing at the Fermi surface in the k = π plane is con- surface in one pair of the octants is further away from
z
sistent with the picture discussed earlier of a quasi-2D the half-filled one so as to share the common modula-
liquid within each domain wall. These effects are ampli- tionwave-vector∆q withtheotherthreepairs. The
111
fied at smaller U’s as shown in Fig. 11, where a larger 111 -SDW state off(cid:104)ers,(cid:105) in this respect, a particularly
(cid:104) (cid:105)
10
√2π 0 √2π this is not energetically too costly. When doing so the
−π π
system benefits from both pairing and band-insulating
∆q
h111i mechanisms to gap the entire Fermi surface and further
q′111 0.9 lower the ground-state energy.
h i
0.8
z0 Q 0 V. DIMENSIONAL CROSSOVER RESULTS
k
q
h111i 1 0.7 A. Results from full numerical HF solutions
∆q
Q 2 h111i
∆qh111i 0.6 Themean-fieldgroundstateofthedoped2DHubbard
model shares many similarities with its 3D counterpart.
π π
−√2π 0 √2−π JustaboveUc,the2DsystemdevelopsasinusoidalSDW
− k 0.5 with a modulating wave along the 10 -direction and a
[110] (cid:104) (cid:105)
√2π 0 √2π much weaker accompanying CDW. As U is increased,
−π π
the SDW increases its amplitude before the SDW state
0.4 eventuallychangesintoacollectionofweaklyinteracting
domainwalls. AboveacertainU,thereisadiscontinuous
transition to a phase where the modulation is along the
0.3
11 -direction. ThecrossoverfromSDWtodomainwalls
(cid:104) (cid:105)
occurs before the 10 to 11 transition at small h, but
kz0 0 0.2 after at larger h20(cid:104). A(cid:105)pec(cid:104)ulia(cid:105)rity of the 2D case, due to
the special topology of the 2D half-filled surface, is that
α=1andthesystemisaninsulatorregardlessofdoping,
0.1 U or direction of the modulation wave-vector apart from
a region close to U .
c
π π Bycontrollingthedistancebetweensquarelatticelay-
−√2π 0 √2−π ers,opticallatticeexperimentsallowthestudyoftheevo-
− k[1¯10] lutionofthesystemasitcrossesoverfrom2Dto3D.This
situationistheoreticallydescribedbyanincreaseoft in
FIG.13: (Coloronline)Schematicillustrationofthepairing Hamiltonian(1)andthequestion,withinmean-field⊥the-
modelforthe(cid:104)111(cid:105)-SDWstate,shownonthecontourplotsof
ory, concerns the ensuing evolution of the ground-state
the (1¯10) (top) and the (110) (bottom) cuts of n for a sys-
k↑ properties. The pairing model and the arguments de-
tem of h=1/8 and U =5.0. The white dashed lines are the
scribedinSec.IV.Cremainvalidinthecrossoverregime.
half-filled surface, across which is the nesting vector Q. The
WethusrestrictourinvestigationtounidirectionalSDW
magenta solid lines show the pairing construction, obtained
by shifting the half-filled surface along [¯1¯11]-direction by a ground states, although we did use the first approach
distance of ∆q /2. In the upper panel, q and q(cid:48) to carry out some searches, finding no additional struc-
(cid:104)111(cid:105) (cid:104)111(cid:105) (cid:104)111(cid:105)
give two equivalent representations (differing by a reciprocal tures. As in 3D, we verify that the SDW solution with
lattice vector) of the pairing vector across the reconstructed minimum energy, identified using the second approach,
Fermisurface. Notetheasymmetrybetweenthetwodiagonal can be obtained by the first approach in a large super-
directions in the upper panel. cell that is commensurate with the optimal wavelengths,
evenwhenstartingfromrandominitialguesses. SDWin
directions different from 100 or 111 are not found to
(cid:104) (cid:105) (cid:104) (cid:105)
clearexamplewhereFermisurfacereconstructioncanbe be the global ground state for any value of t .
observed. It also shows how an accurate experimental Results are summarized in the t -U mean⊥-field phase
characterization of the momentum distribution in opti- diagram of Fig. 14. An overall incr⊥ease in the critical U
callatticescanbeusedtocharacterizethebandstructure valuesisseenasdopingincreases, asaresultofagreater
andpairingattheFermisurfacewhich, inturn, provides deformation from the perfectly nested half-filled Fermi
momentum space evidence onthe real space character of surface and the need for more excitations to achieve re-
the SDW. construction of the Fermi surface for pairing. As before,
Using Eq. (17) we have estimated the wavelength of numerical calculations focus on small doping (h (cid:54) 0.2)
the 111 -SDW,andfoundα =0.93. TheexactSCF and low to intermediate interactions (U (cid:54) 5.5), where
111
calc(cid:104)ulati(cid:105)ons in Sec. IVB sho(cid:104)wed(cid:105), instead, that α is mean-field theory can be expected to be more accurate.
111
precisely pinned at 1 in a fairly large regime of U(cid:104). T(cid:105)his Upon increase of U, and similarly to 3D, the system un-
somewhat large discrepancy is a consequence of the nat- dergoes a first transition to a 100 -SDW state followed
(cid:104) (cid:105)
uraltendencyofthesystemto“lock”theintegratedden- by a second, discontinuous transition to a 111 -SDW
(cid:104) (cid:105)
sityofholespernodalregionto1wheneverthetopology state. The absence of cubic symmetry away from t =1
of the non-interacting doped Fermi surface is such that causes the modulation wave-vector for 100 -SDW⊥to lie
(cid:104) (cid:105)