Hybrid-space density matrix renormalization group study of the doped two-dimensional Hubbard model

0.82 MB·

by  G. Ehlers

Hybrid-space density matrix renormalization group study
of the doped two-dimensional Hubbard model
G. Ehlers,1 S. R. White,2 and R. M. Noack1
1 Fachbereich Physik, Philipps-Universita¨t Marburg, 35032 Marburg, Germany
2 Department of Physics and Astronomy, University of California, Irvine, California 92697, USA
(Dated: March 22, 2017)
The performance of the density matrix renormalization group (DMRG) is strongly influenced
by the choice of the local basis of the underlying physical lattice. We demonstrate that, for the
two-dimensional Hubbard model, the hybrid–real-momentum-space formulation of the DMRG is
computationally more efficient than the standard real-space formulation. In particular, we show
7 that the computational cost for fixed bond dimension of the hybrid-space DMRG is approximately
1 independentofthewidthofthelattice,incontrasttothereal-spaceDMRG,forwhichitispropor-
0 tional to the width squared. We apply the hybrid-space algorithm to calculate the ground state of
2 the doped two-dimensional Hubbard model on cylinders of width four and six sites; at n=0.875
filling, the ground state exhibits a striped charge-density distribution with a wavelength of eight
r
a sites for both U/t=4.0 and U/t=8.0. We find that the strength of the charge ordering depends
M onU/tandontheboundaryconditions. Furthermore,weinvestigatethemagneticorderingaswell
as the decay of the static spin, charge, and pair-field correlation functions.

1
2 PACSnumbers: 71.10.Fd,71.27.+a

]
l I. INTRODUCTION asmultiscaleentanglementrenormalization(MERA)[15,
e
- 16] and projected entangled pair states (PEPS) [17, 18]
r avoid this restriction in principle by adapting the topol-
t Although the one-band Hubbard model [1] in two di-
s ogy of the tensor network to the entanglement structure
. mensionshaslongbeentoutedasaleadingcandidatefor
t of the system. However, these methods suffer from the
a explainingthebasicphenomenonofhigh-temperaturesu-
m perconductivity in copper oxide planes [2], whether the factthatthescalingoftheircomputationalcostswiththe
dimension of the Hilbert space treated is much higher
- unmodified model provides sufficient features to do this
d is still controversial [3]. In order to clear up this issue, than that of the standard DMRG or other MPS-based
n algorithms.
enormouseffortisbeingmadetoobtainitsphasediagram
o
numerically [4]. As the model is doped away from half Another ansatz to broaden the applicability of the
c
[ filling,richbehavioremerges[5],withphasesthatinclude DMRGfortwo-dimensionalsystemsis,insteadofchang-
antiferromagnetic ordering near half filling, a metallic ingthetopologyoftheunderlyingnetwork,tochangethe
2
phase for weak on-site interaction, and a superconduct- local basis of the physical model [19]. A change of basis
v
ing phase for moderate interaction. Many details of the can influence three characteristics that drastically influ-
0
9 phase diagram, such as the presence of charge and spin ence the performance of the DMRG: the entanglement
6 density stripes [6, 7], which have been shown experimen- betweensubsystems,therangeandstructureoftheinter-
3 tally to play a role in high-temperature superconducting actions, and the number of quantities that are conserved
0 materials [8, 9], are yet to be unequivocally determined. within the specific representation.
.
1 One of the most important numerical methods used For instance, for translational invariant models, the
0 to address these questions is the density matrix renor- momentum-space DMRG takes advantage of the con-
7
malization group (DMRG) [10–12]. In particular, the served momentum quantum number and achieves a sig-
1
DMRG can provide unbiased results that are a very use- nificant speedup for a fixed size of the truncated Hilbert
:
v fulbenchmarkandcomplementforothermethods. While space [20, 21]. However, the scaling of the size of the
Xi the DMRG has been very successful in treating one- Hilbert space needed to maintain a given fixed accuracy
dimensional systems and ladders, its application to two- is problematic. In particular, the block entropy is zero
r
a dimensional systems, such as wider cylinders, is much in the noninteracting limit in the momentum-space rep-
more challenging [13] because the growth of the entan- resentation, as the noninteracting Fermi sea is a prod-
glement between the DMRG blocks, which is generally uct state. Unfortunately, short-range interactions in real
proportionaltothesystemwidthforashort-rangeHamil- space become long-range in momentum space so that
tonian,leadstoexponentialgrowthofthecomputational all parts of the system become strongly entangled and
cost. Furthermore, the DMRG treats two-dimensional the block entropy increases rapidly as the interaction is
systems by mapping the lattice to an intrinsically one- turned on; as a result, the entropy scales with the vol-
dimensional matrix product state (MPS) [14], resulting ume of the system for all nonzero interaction strengths.
in longer-range effective interaction even for short-range Thus, the computational cost scales exponentially in the
models. volume of the system, making treating systems of any
More recently developed tensor network methods such significant size prohibitively expensive [22].
2
Recently, it has been shown that, by choosing a mixed tiperiodicboundaryconditions,k =2π(j+ 1)/L . The
y 2 y
basis, one can partially take advantage of the perfor- resulting Hamiltonian in hybrid space,
mance benefits of the momentum-space DMRG but also
(cid:88)
retain the beneficial entanglement scaling of the real- H = t c† c (3a)
− xkyσ x(cid:48)kyσ
space representation [23]. Our goal hence is to fur- (cid:104)x,x(cid:48)(cid:105)kyσ
ther explore this approach. In particular, we apply (cid:88)
+ ε(k )n (3b)
the hybrid–real-momentum-space DMRG to the two- y xkyσ
dimensional Hubbard model on a lattice with cylindrical xkyσ
topology. + U (cid:88) c† c† c c ,
The remainder of the paper is organized as follows: in Ly xky+qy↓ xpy−qy↑ xpy↑ xky↓
Sec. II, we express the Hubbard model in the hybrid- xkypyqy
(3c)
space representation, discuss the structure of the ma-
trixproductoperator(MPO),andoutlinehowreal-space
consistsofthreeterms,alongitudinalhoppingterm(3a),
two-point correlation functions can be calculated in hy-
a transverse hopping term (3b) with dispersion rela-
brid space. In Sec. III, we discuss the computational
tion ε(k )= 2tcosk , and a long-range interaction
cost of the hybrid-space DMRG as applied to the two- y − y
term (3c). Note that the on-site Hubbard interaction
dimensional Hubbard model. In particular, we analyze
becomes long-range in the momentum direction, but re-
thescalingoftheCPU-timeandmemorycostsasafunc-
mains short-range in the real-space direction.
tion of the cylinder width, and verify the results us-
As will be described in Sec. IV, we find states with
ing realistic numerical calculations. Section IV then de-
stripedchargeandspindensitypatternsforsystemsthat
scribes our study of the ground state of the doped Hub-
are moderately doped away from half filling. In order to
bardmodelatfillingn=0.875andinteractionstrengths
stabilize and target states with a particular wavelength
U/t=4.0andU/t=8.0oncylindricallatticesofwidth4
of the charge-density stripes, λ, we sometimes apply an
and6. Inparticular,weaddressthequestionsofwhether
additional pinning field
the ground state exhibits stripe structures and whether
the pairing correlations are enhanced. Finally, Sec. V (cid:88)
H =P sin(φ+κx)n (4)
contains the conclusion. n xkyσ
xkyσ
tothedopedsystem. Thefieldcouplesdirectlytothelo-
II. MODEL AND METHOD
calchargedensityinhybridspaceandcanbetailoredus-
ing its amplitude P, wavenumber κ=2π/λ, and phase
A. Hubbard model in hybrid space
φ. Depending on the stability of the target state, we
either turn of the pinning field after the initial sweeps
We investigate the two-dimensional Hubbard model of the DMRG, or we keep the field amplitude finite but
with nearest-neighbor hopping and on-site Coulomb re- small throughout the calculation and subtract the con-
pulsion defined by the Hamiltonian tributionofthefieldtotheground-stateenergyoncethe
calculation is completed. In order to ascertain that the
(cid:88) (cid:88)
H = t c† c +U n n , (1) presence of a small pinning field has no significant effect
− rσ rσ r↑ r↓
(cid:104)r,r(cid:48)(cid:105)σ r on the physical results, we have compared calculations
with and without a pinning field for cases in which the
where r,r(cid:48) denotes nearest neighbors on a square lat- pinning field can be ramped to zero during the calcu-
tice wi(cid:104)th la(cid:105)ttice sites r=(x,y). Here, crσ (c†rσ) are lation without affecting the final convergence and have
creation (annihilation) operators for electrons with spin found the difference in observables and in energy (with
σ ∈{↑, ↓},andnrσ =c†rσcrσ istheparticle-numberop- the energy of the pinning field subtracted) to be negligi-
erator. We take the lattice geometry to be cylindrical, ble.
with cylinder length L , width (i.e., circumference) L ,
x y
a lattice spacing of unity, and periodic or antiperiodic
boundary conditions in the transverse direction. Such a B. Hybrid-space matrix product operator
geometry is favorable for hybrid–real-momentum-space
representation. We Fourier transform in the transverse
In its modern formulation, the DMRG algorithm
(y-) direction, i.e., write
is best described within the framework of MPSs and
MPOs [24, 25]. The MPS and MPO store the coeffi-
1 (cid:88)
c†rσ = (cid:112)L e−ikyy c†xkyσ (2) cientsofthestate|Ψ(cid:105)andHamiltonianH asproductsof
y ky matrices associated with each site of the DMRG chain.
TheMPSandMPObondsarethecontractionsoftherow
andcorrespondinglyforc . Forperiodicboundarycon- andcolumnindices,respectively,ofneighboringmatrices.
rσ
ditions, the transverse momentum points are given by ThedimensionsoftheMPSbonds,m,aredirectlyrelated
k =2πj/L with integer 0 j <L , whereas, for an- to the number of states kept in the traditional DMRG,
y y y
≤
3
andthedimensionsoftheMPObondsrefertotheopera- der (dashed line B in Fig. 1), the corresponding MPO
torsstoredintheleftandrightblockoftheDMRG.The bond has to carry all the long-range interactions in the
MPO also encodes the rules of how H Ψ must be cal- term(3c);foreachringthissumrunsoverthreeindepen-
| (cid:105)
culated within the DMRG-specific block-site-site-block dent momenta. If we were to construct the MPO with
decomposition of the system and how the DMRG blocks individual channels for each term, the resulting bond di-
are updated during the sweeping process. Therefore, it mension would be proportional to L3.
y
is crucial to find the optimal MPO representation for a Fortunately, for a given bipartition, the sum (3c) can
given model and lattice, i.e., to find the MPO with min- be factorized to reduce the effective number of terms as
imal bond dimension. follows: for fixed x, the part of (3c) for which the anni-
Roughly speaking, the dimension of each bond of the hilation and creation operators are in separated parts of
MPO depends on the number of terms of the Hamilto- the system can be rewritten as
nian acting simultaneously on both sites of the bond. If
the Hamiltonian is factorizable in an appropriate way, U (cid:88)A† A (5)
the bond dimension can be minimized by accumulating N xky xky
all interactions between opposing sides of a “composed” ky
operator, so that the same interaction can be expressed with composed operators
using fewer terms [20]. In this section, we describe how
(cid:88)
thisisdoneinprinciple;thedetailsoftheconstructionof A† = c† c† ,
the MPO matrices for the hybrid-space Hubbard Hamil- xky xpy↓ xky−py↑
py
tonian are given in Appendix A. (cid:88)
A = c c . (6)
For the two-dimensional Hubbard model in real space, xky xpy↑ xky−py↓
the bond dimension of the optimal MPO is proportional py
to the system width L . This can be understood as fol-
y Analogous steps can be carried out for all possible dis-
lows: Figure 1 shows the most common mapping of a
tributions of creation and annihilation operators of both
spin species onto the two subsystems, resulting in a for-
mulation in which the total number of terms connecting
the two parts of the system is (L ). Thus, without ap-
y
O
proximation, we obtain an MPO in which the dimension
ofeachbondisproportionaltothesystemwidth, justas
in real space.
This approach was introduced for the momentum-
A B space DMRG in Ref. [20]. The technique can be applied
to other models as long as the Hamiltonian is factoriz-
able. The hybrid-space Hamiltonian and the optimized
FIG. 1. (Color online) Mapping of the one-dimensional
MPO for the fermionic Hofstadter model are described
DMRG chain (black solid line) onto a 10×4 square lattice.
The bold gray lines depict nearest neighbors on the square in detail in Ref. [23]. Similar optimizations of MPOs
lattice,andthedashedlinesAandBdepicttwopossiblecuts have also been carried out for other systems with long-
throughtheDMRGchainandthecorrespondingbipartitions rangeinteractions,suchasthequantumchemicalHamil-
of the system. tonian [26].
The exact bond dimension for the Hubbard model in
two-dimensional square lattice onto the one-dimensional hybrid space varies with the position of the cut in the
MPS / MPO chain; the one-dimensional path is simply cylinder, whereas in real space it is constant within the
folded over the width of the lattice into two dimensions. bulk of the system. Assuming a constant MPS bond di-
Therefore, it is clear that, whenever the chain is cut, the
number of nearest-neighbor bonds, and thus the number
TABLE I. Average MPO bond dimension for the two-
of hopping term in the two-dimensional Hubbard model
dimensional Hubbard model in the real-space and in the
in real space, is proportional to the width of the system.
hybrid-space representation for different cylinder widths.
Since (almost) all of these terms act on different sites of
the lattice, the MPO must include an individual channel Ly =4 Ly =6 Ly =8 Ly =10
foreachterm; thusitsbonddimensionisproportionalto Real space 18 26 34 42
the system width.
Hybrid space 26.0 45.7 66.5 84.6
In hybrid space, the situation is more complicated. In
particular, the bond dimension of the MPO depends on
wherewecutthesystem: ifthecutisbetweentwoneigh- mension, a good indicator for the computational cost of
boringringsofthecylinder(dashedlineAinFig.1),only the DMRG is the average MPO bond dimension, which
real-space-like hopping terms (3a) are cut, and the bond is approximately twice as large in the hybrid-space rep-
dimension is again proportional to the system width. If resentation as in the real-space case (Table I). Note that
the cut separates sites of the same ring of the cylin- the averaged MPO bond dimension for the hybrid-space
4
MPO may still vary slightly for different orderings of size has a decisive influence on the performance for fixed
the momentum points of each cylinder ring within the accuracy. For two-dimensional systems to which the en-
DMRG chain. tropy area law [27] applies, the entropy is proportional
tosystemwidthL . Theresultingexponentialscalingof
y
thecomputationalcostwithsystemwidthisafundamen-
C. Real-space correlation functions in hybrid space tallimitationofthehybrid-spaceDMRGaswellasofthe
real-spaceDMRG andallknownMPS-basedalgorithms.
In this section, we at first neglect the variation in the
The equal-time spin, charge, and pair-field correlation
required bond dimension m on model parameters and
functions are defined as
lattice size and analyze the scaling of the computational
S(r,r(cid:48))= S (r)S (r(cid:48)) , cost of the real-space and hybrid-space DMRG for fixed
z z
(cid:104) (cid:105) m. Wethusestimatetheperformancegainofthehybrid-
C(r,r(cid:48))= n(r)n(r(cid:48)) n(r) n(r(cid:48)) ,
(cid:104) (cid:105)−(cid:104) (cid:105)(cid:104) (cid:105) space relative to the real-space representation. We then
Dyy(r,r(cid:48))=(cid:104)∆†y(r)∆y(r(cid:48))(cid:105) (7) compareourestimatetomeasurementsoftheactualrun-
time and memory usage for typical calculations and fi-
with nally come back to the issue of the dependence of the
relative accuracy of the hybrid-space and real-space rep-
S (x,y)=c† c c† c ,
z xy↑ xy↑− xy↓ xy↓ resentations on bond dimension m.
n(x,y)=c† c +c† c ,
xy↑ xy↑ xy↓ xy↓
1
∆†(x,y)= (c† c† c† c† ), (8) A. Estimated scaling of the computational cost
y √2 xy+1↑ xy↓− xy+1↓ xy↑
whereS (x,y)measuresthelocalspin,n(x,y)isthelocal The majority of the computational cost in the DMRG
z
charge density, and Ơ(x,y) is the creation operator for algorithm comes from applying the Hamiltonian to a
y
a pair of spin-up and -down particles on sites (x,y) and state, H Ψ ; this is the fundamental step in itera-
| (cid:105)
(x,y+1). tive eigensolvers such as the Lanczos or Davidson algo-
Measuringtwo-pointreal-spacecorrelationfunctionsin rithms [28, 29]. Therefore, we examine the scaling of the
hybrid space raises the same issues as implementing the operationsrequiredforthisoperationwiththedimension
MPO of the Hamiltonian: applying the Fourier trans- of the physical lattice sites, d, cylinder width Ly, cylin-
formation (2) to the correlation functions (7) introduces derlengthLx,andMPSbonddimensionmbasedonthe
sumsovermultiplemomenta similartoterm(3c). These structures of the MPS and MPO used. The computa-
sums can again be factorized analogously to Eq. (5): the tional costs of other operations, e.g., changing the active
pair-field correlation functions in hybrid space can be sites of the DMRG in the sweeping process, benefit from
written as the hybrid-space representation in the same way.
We assume that the bond dimension of the MPO is
Dyy(r,r(cid:48))= 21L2y (cid:88)ky eiky(y−y(cid:48)) (cid:104)Oy†(x,ky)Oy(x(cid:48),ky)(cid:105) prerporpeosretniotantailontso,Lasyisinthbeoctahserefoarl-tshpeacHeuabnbdarhdymbroidd-eslp(asecee
Sec. IIB). In our estimate, we neglect the possibility of
(9)
exploiting the conservation of spin and charge quantum
with composed operators
numbers, as it would have the exact same effect on the
O†(x,k )=2 (cid:88)cos(p k /2)c† c† , computational costs in both representations.
y y y− y xpy↑ xky−py↓ Inthereal-spacerepresentation,oneLanczosstepthen
py requires (d2L ) multiplications of m m matrices, re-
(cid:88) y
Oy(x,ky)=2 py cos(py−ky/2)cxpy↑cxky−py↓. (10) sOu(ldti2nLg2yinLOxOm(d32KL)ymop3e)raotpioenrastipoenrsDpMerRL×Ganscwzoeseps,tewpitahnda
fixed number of Lanczos steps per DMRG step, K. The
Thespinandchargecorrelationfunctionscanbetreated corresponding memory costs are (d2m2K) for the K-
analogously. Thus, the correlation functions (7) can be dimensional Krylov space and (OL m2) for the left and
y
measured without changing the scaling of the total com- O
right DMRG block in the current block-site-site-block
putational cost of the algorithm.
configuration.
Next, we analyze the scaling of the costs in the
hybrid-space representation under the assumption that
III. PERFORMANCE thehybrid-spaceMPShasthesametotalbonddimension
m. Since the Hamiltonian (5) conserves the transverse
Herewestartbyemphasizingthatthedimensionofthe momentum quantum number, every matrix of the MPS
MPS bonds required to obtain a given accuracy scales can be written in a block-diagonal form with L blocks
y
exponentially with the block entropy; thus, the scaling ofsizem(cid:48) m(cid:48) withm(cid:48) m/L ,whereeveryblockcorre-
y
× ≈
of the block entropy with system parameters and lattice sponds to one momentum sector. Therefore, the compu-
5
tational costs reduce to O(d2L−y1m3) per Lanczos step 102 102
and (d2Lxm3K) per DMRG sweep, and the memory m3 m2
cOo(smtsO2b)efocormtheeOD(Md2RLG−y1bmloc2kKs.)Afosridteh-ebyK-sriydloevcosmpapcaeriasnond [h] 101 ∝hreyablrid B] 101 ∝hreyablrid
of all costs is given in Table II. p G
Note that the above argument is only valid because ee 100 [
w
cnthounemsMberePvreSsdebcqotnourdasnsotdifteaicepospmrspouxochismeaainstteotlyhLeeyqtumoatolamlsizesepnitmnum(cid:48).orFqoupraanortttihucemler /swall10−1 RAM 100
t (a) (b)
number, there typically is a different distribution of the
quantum number sectors, with only a few large sectors 10 2 10 1
dominating the computational costs of the algorithm; in − 103 104 − 103 104
thesecases,thespeedupdependsprimarilyonthesizeof m m
the largest sectors rather than on the number of sectors. 30 5
(c)
25
4
TABLE II. Scaling of the runtimes of a single Lanczos step,
g
TmLeamncozroys,caonstdsoafssaocDiaMteRdGwistwhetehpe,KTrsywleoevp,spanacde,scMalKinryglovo,fatnhde up20 savin3 (d)
the left and right DMRG blocks of the current block-site- d15
e y
site-block configuration, Mblock, as a function of the cylinder pe or2
lmen,gtthheLloxc,atlhelatctyilcienddeirmwenidstiohnLdy,,athnedMthPeSdbimonendsidoinmeonfstiohne s10 Ly=4 mem Ly=4
Krylov subspace, K. 5 Ly=6 1 Ly=6
Ly=8 Ly=8
Real space Hybrid space 0 0
0 10k 20k 0 10k 20k
T (d2L m3) (d2L−1m3)
Lanczos O y O y m m
T (d2L2L m3K) (d2L m3K)
sweep O y x O x
M (d2m2K) (d2L−1m2K) FIG. 2. (Color online) Performance comparison between
MKrylov O (L m2) O (ym2) real-spaceandhybrid-spaceDMRG,calculatedforU/t=4.0
block O y O and n=0.875; (a) wall time per DMRG sweep and (b) peak
memoryconsumptionfora16×6cylinderasafunctionofthe
In conclusion, the computational costs of the hybrid- MPSbonddimensionm. Thedashedgraylinesdepictthem3
and m2 scaling expected in the m→∞ limit. (c) Speedup
space DMRG are expected to be independent of L for
y
and(d)memorysavingsofthehybrid-spaceDMRGcompared
fixed m, whereas for the standard real-space DMRG the
to real-space DMRG as a function of m for 16×4, 16×6, and
runtime and the memory consumption scales as L2 and
y 16×8 cylinders.
L , respectively (Table II). Estimating the bulk MPO
y
bond dimension for the Hubbard model to be approx-
imately twice as large in hybrid space as in real space
(Table I), we estimate a total speedup of roughly L2/2. for large m in both cases, while the peak memory con-
y
sumption scales with m2 [Figs. 2(a) and 2(b)]. The de-
viation from this limiting behavior for smaller m is due
B. Measured performance tothequantumnumberbookkeepingandotheroverhead
in the code and is amplified for the hybrid-space DMRG
because of the additional momentum quantum number.
Inordertoinvestigatetheactualperformance,wecom-
In agreement with the predictions above, the speedup
pare the computational costs of real-space and hybrid-
of the hybrid-space DMRG over the real-space calcula-
space calculations for Hubbard cylinders with length 16,
tions is larger for wider cylinders [Fig. 2(c)]. Because of
widths 4, 6, and 8, and periodic transverse boundary
the additional overhead, the full speedup of the hybrid-
conditions. All calculations were carried out using the
same code and using six physical cores on an Intel(cid:13)R spacecodeisonlyseenatlargem. Intheregimeinwhich
Xeon(cid:13)R X5650 CPU. In both representations, we exploit bothmethodsstillprovideresultswithinreasonabletime
in our calculations, the hybrid-space DMRG is approx-
the block-diagonal structure of the MPS and MPO ma-
imately 12 times faster for L =4 and up to 20 and 26
triceswithrespecttotheconservedchargeandtotalspin y
timesfasterforL =6andL =8. Intermsofthepeak
quantum numbers; in the hybrid space version, we fur- y y
memory consumption, we observe a weaker influence of
ther decompose the matrices using the transverse mo-
L : the memory savings for the larger m values treated
mentum quantum number. y
varies only between a factor of 4 for L =4 and 4.5 for
Figure2showsacomparisonoftheruntimeandmem- y
L =6 and L =8 [Fig. 2(d)].
ory requirements for different L as a function of m. As y y
y
expected, the CPU time per sweep is proportional to m3 In Fig. 3, we compare the CPU time per sweep and
6
15
4
L =4 L =4
y y
h] 10 Ly =6 Ly =6 )3
p[ 7.5 Ly =8 GB]10 Ly =8 S(i2
e [
/sweall 5 RAM 5 01 (a) hreyablr-sidp-ascpeace
w 0 12 24 36 48
t 2.5
(a) (b)
5
0 0
10k 20k 30k 10k 20k 30k
4
m m
)
i3
(
S
2
FIG.3. (Coloronline)(a)CPUtimepersweepand(b)peak
hybrid-space
memoryconsumptionofhybrid-spaceDMRGcalculationsfor 1 (b)
real-space
Hubbardcylinderswithlength16andwidthLy asafunction 0
oftheMPSbonddimensionm. Thecalculationswerecarried 0 24 48 72 96
out for n=0.875 and U/t=4.0.
i
FIG. 4. (Color online) Von Neumann entropy S(i) as a
thememoryconsumptionofthehybrid-spaceDMRGfor function of the MPS bond index i for (a) 8×6 cylinders at
different L . In agreement with Table II, the computa- U/t=4.0andhalf-fillingand(b)16×6cylindersatU/t=8.0
y
and n=0.875. The sites of the real-space and hybrid-space
tional costs are almost independent of L . The minor
y
lattices are mapped to the MPS sites in an x-direction-first
growth of the runtime that is still present is clearly sub-
ordering;accordingly,thegrayverticallinesindicatecutsbe-
linear in L . The measured computational costs deviate
y tween neighboring rings of the cylinders.
somewhat from the estimated costs even in the large-m
limit, in which the influence of overhead should be in-
significant. Inparticular,themeasuredabsoluteruntime
of the hybrid-space DMRG is not strictly independent as U/t is increased and the system is doped. Figure 5
of L , and thus the observed speedup over real-space illustrates the good agreement of the ground-state en-
y
DMRG grows more slowly than L2. This deviation is ergies for fixed m as well as after extrapolation to zero
y
caused by fact that the MPO dimension is not perfectly truncation error for both parameter sets in Fig. 4. The
proportional to L for small L (Table I). slight divergence in the extrapolations for the hybrid
y y
versus the real-space calculations, especially evident in
After having compared the computational costs for
Fig. 5(b), shows the limitations of the extrapolation, es-
fixed m, we now investigate how the change of basis in-
peciallyforthereal-spaceresults. Ascanbeclearlyseen,
fluences the block entropy and the convergence of both
the maximum MPS bond dimension available is signif-
methods. Since the MPS bond dimension required to
icantly larger for the hybrid-space algorithm, resulting
reach a fixed truncation error grows exponentially with
in a more well-controlled extrapolation and higher accu-
the von Neumann entropy of the DMRG blocks, S(i),
racy in the extrapolated energy. In Fig. 5, we expect
even small changes in the entropy can influence the con-
that the real-space results would converge towards the
vergence significantly. Figure 4 shows comparisons of
hybrid-space results if m could be further increased in
S(i) in real and hybrid representations for the half-filled
the real-space calculation; this is not practically possi-
system at U/t=4.0, Fig. 4(a), and for the doped sys-
ble here. For narrower systems, the energies converge
tem, n=0.875, at U/t=8.0, Fig. 4(b). In Fig. 4(a), it
more rapidly with m, and both methods yield results
can be seen that the block entropy in the hybrid-space
that agree very well. For example, for the 16 4 cylinder
representation differs only slightly from that in the real-
×
at U/t=8.0 and n=0.875, we obtain the same extrap-
space representation. Note that if we cut the system
olated energy, E /t= 0.75114(2), for each method.
between neighboring rings of the cylinder (gray lines), 0
−
the block entropy actually is the same in both cases, as Even with the significant reductions in computational
expected. This can also be seen in Fig. 4(b); however, costs, achieving good convergence for wider cylinders is
here the entropy in the hybrid representation is percep- still very expensive, making an efficient parallelization
tibly higher between these points. This illustrates that indispensable. In this respect, the hybrid-space DMRG
the nonlocal nature of the interaction within the rings hasanadvantageousproperty: theadditionalmomentum
does lead to a moderate increase of the entropy for cuts quantumnumberleadstofiner-grainedquantumnumber
within the rings in the hybrid representation, especially sectors,whichmakesforbetterloadbalancing. Inpartic-
7
atstrongcoupling. (Real-spaceDMRGactuallybecomes
hybrid-space
0.822 exactintheatomiclimitratherthanthestrong-coupling
− real-space limit,butalocalinteractionthatisstrongrelativetothe
1
E/t0 −00..882234 (a) 30.02k7.5k25.0k22.5k20.0k17.5k 1250..00kk 17.52k.5k 15.0k hAhcfouoasrtpsswpctbuihantesegstwdsgwaieesmeicntnuheesricbsnayellllodtiynchikdbenereirrSnnienrtgcirgsn.osgtpIhsyI(eIsaeaBnsesyd,sFtttmhhiegeemo.drhe4ceyal)rob,al-srtsaeiepdrnae-tsxceopecxatecmhcseseasetsetDhnlwiMtomrhdoRiiptcaG.hyt)
−
becomes smaller for weaker interaction strength. There-
0 1e-5 2e-5 3e-5 fore, the hybrid-space algorithm has slightly better rel-
0.740 ative convergence at moderate interaction strength and
E/t0 −−0.745 (b)hreyablr-s3i0dp.02-ka7s.5c2pk5e.a02k2c.5ek20.0k17.5k15.0k15.102.5k12.5k1100.0.0kk 8.0k csmthpaoneaIrncseaeitltasmheoceoechtecfohascslnuoliodpdbwla,leensinoactgyhtl,,aahawrnagnteewdtrciheatbqehlocumnutahdlola-edttdieerimmertaaehetl-enessp-spgiiinroantonc,euecrmnhamdact-terhsigtothaean,notadreten.hegdenimeprreaegairyils--,
k
fieldcorrelationfunctions. Allresultsareextrapolatedto
zero truncation error. The maximum MPS bond dimen-
0.750
− sion was at least 30000 for cylinder lengths 16 and 32,
0 2e-5 4e-5 6e-5 27500 for length 48, and 25000 for length 64.
∆ξ
FIG. 5. (Color online) Ground-state energy obtained from 0.95 (a) 32×4, PBC, λ = 8.0
hybrid-space DMRG and real-space DMRG as a function of )
x
the discarded weight per site, ∆ξ, for (a) 8×6 cylinders at (0.875
n
U/t=4.0andhalf-fillingand(b)16×6cylindersatU/t=8.0
and n=0.875. The MPS bond dimension is increased every 0.8
other sweep and is written alongside the corresponding data
points. The solid lines are linear fits through the last five 0.95 (b) 32 4, ABC, λ = 8.0
×
data points of each series and indicate the zero-truncation- )
x
error extrapolations of the ground-state energies. (0.875
n
0.8
ular, the momentum quantum numbers yield equal-sized
0.95 (c) 32 6, PBC, λ = 8.0
quantum number sectors, resulting in more equal-sized ×
chunksofwork. Inourimplementation,weapplyshared- )
x
memoryparallelizationtothecombinedloopsoverquan- n(0.875
tum number sectors and terms of the Hamiltonian, as
described for classical DMRG in Ref. [30]. In order to 0.8
further extend the applicability of the hybrid algorithm, 0.95 (d) 32 6, ABC, λ = 8.0
additional steps such as real-space parallelization [31] or ×
)
efficient distributed memory parallelization [32] could be x
(0.875
implemented. n
0.8
0 8 16 24 32
IV. RESULTS
x
We study the ground state of the doped Hubbard FIG.6. Charge-densitydistributionn(x)forlength-32Hub-
model on width-4 and width-6 cylinders for U/t=4.0 bard cylinders at U/t=4.0 and n=0.875. Panels (a)–(d)
and U/t=8.0 at filling n=N−1(cid:80)rσnrσ =0.875. We showwidth-4andwidth-6cylinderswithperiodic(PBC)and
choose these two values of U/t for the following reasons: antiperiodic (ABC) transverse boundary conditions, as indi-
(i) In general, we want to compare the physics and the cated. Theblacklinesshowthezero-truncation-errorextrap-
performance of the method at moderate coupling with olateddensities,andthegraylinesshowthenon-extrapolated
values. The gray dashed lines indicate the average filling.
thatatstrongcoupling. (ii)Itisinterestingtoinvestigate
the stability of inhomogeneous phases as the interaction
strengthischanged,inparticular,whetherastripephase Figure 6 shows the charge-density distribution in
remainsstableforweakerinteraction. (iii)Forreal-space the longitudinal direction, n(x)=L−y1(cid:80)yσnxyσ, for
DMRGmethods,thenumericalconvergence,i.e.,thebe- U/t=4.0 and different widths L and boundary con-
y
havior of the block entropy, is generally more favorable ditions. In all cases, we find a stripe pattern with
8
∆ξ 0.0,buttheλ=5.3statehaslowerenergyforany
0.95 (a) 32 4, PBC, λ = 8.0 acce→ssible fixed m. Because of this, the λ=8.0 state be-
×
) comesunstableifP issettozeroatanypointduringthe
x
(0.875 calculation due to the fact that the DMRG algorithm
n
is always driven towards the state that has the lowest
0.8 energy within the finite reduced subspace of the Hilbert
0.95 (b) 32 4, ABC, λ = 8.0 space treated, i.e., the state space with fixed MPS bond
× dimension m. The ground state within this subspace is
x) not necessarily the ground state for fixed truncation er-
(0.875
n ror ∆ξ or, indeed, for ∆ξ 0. In general, the DMRG
→
shouldconvergetothetruegroundstateatsomelevelof
0.8
accuracy(becauseitbecomesexactform ),butthe
0.95 (c) 32 6, PBC, λ = 8.0 required m may be inaccessible. Thus, a→ppl∞ying a pin-
×
ningfieldisnecessaryhereinordertoenabletheDMRG
)
x
(0.875 to track both states as the discarded weight is varied;
n
in practice, we use the fixed pinning field amplitude of
P = 0.01 for the depicted range of ∆ξ for both states
0.8
even though it would have been possible to set P =0 in
0.95 (d) 32 6, PBC, λ = 5.3
× the latter part of the λ=5.3 calculation.
)
x
(0.875 0.740
n
−
λ 5.3 P =0.01
0.8 ≈
λ 8.0 P =0.01
0 8 16 24 32 0.742 ≈
x −
2
5.
FIG.7. Charge-densitydistributionn(x)forlength-32Hub- t 0.744 20. 3 0k
bshaorwd ccyylliinnddeerrswaitthUd/iffte=re8n.t0wainddthn,b=ou0n.8d7a5r.ycPoanndeitlsion(as),–a(ndd) E/0 − 35.300.02k5.0k0k 35.0k0.0k
wavelengthλofthecharge-densitystripes,asindicated. The k
0.746
black lines show the zero-truncation-error extrapolated den-
−
sities, and the gray lines show the non-extrapolated values.
The gray dashed lines indicate the average filling.
0.748
−
0 2e-5 4e-5 6e-5 8e-5
wavelength 8; for Ly = 4, each stripe contains 4 holes ∆ξ
[Figs. 6(a) and 6(b)], and for L = 6 each stripe con-
y
tains 6 holes [Figs. 6(c) and 6(d)]. For larger interaction
FIG.8. (Coloronline)Discardedweightextrapolationofthe
strength, U/t=8.0, Fig. 7, we find two different stripe
energy for λ = 5.3 and λ = 8.0 on 16×6 cylinders at U/t =
configurations, wavelength λ=8.0 [Figs. 7(a)–7(c)] and 8.0 and n = 0.875. The energy contribution of the pinning
wavelength λ=5.3 (more exactly, λ=16/3) [Fig. 7(d)]. fieldhasbeensubtracted,andtheMPSbonddimensionmis
For both configurations, we add a pinning field, Eq. (4), indicated by the labels alongside selected data points.
withtheappropriatewavelengthtostabilizethestate. In
most cases it is sufficient to add the pinning field during The amplitude of the charge-density modulations,
the initial sweeps of the calculation and ramp its am- ∆n=1/2 max [n(x)] min [n(x)] , is plotted as a
x x
plitude P down to zero in subsequent sweeps. However, function o{f the invers−e cylinder le}ngth in Fig. 9 for
for width 6 cylinders at U/t = 8.0, the amplitude of the both values of U/t. The trend shows that in all cases
field must be held finite (we take P =0.01) during the the stripes have finite amplitude in the infinite-cylinder
entire calculation in order to stabilize the λ=8.0 stripe length limit. The fact that the stripes are enhanced
pattern. After subtracting the field energy, we find that for L =6 indicates a striped ground state in the two-
y
the λ=8.0 phase is lower in energy and is thus globally dimensional limit. As the interaction is increased from
stable; it should also be noted that the λ=5.3 phase U/t=4.0 to U/t=8.0, the amplitude of charge density
only occurs for the Ly =6 lattices with periodic bound- in the λ=8.0 stripes increases.
ary conditions. Typically, a striped charge-density distribution in the
In order to elucidate the convergence of the λ = 5.3 doped Hubbard model is accompanied by a striped stag-
and λ = 8.0 states for L = 6 at U/t = 8.0, we examine gered spin-density distribution with a wavelength that is
y
thediscarded-weightextrapolationofbothpinning-field- double that of the charge-density distribution [6]. Mea-
stabilized states in Fig. 8. As can be seen, the λ = 8.0 suring these stripes in hybrid space would require break-
state has lower energy for any fixed ∆ξ as well as for ingthetranslationinvarianceinthetransversedirection,
9
3
0.06 LLyy==44,,PABBCC,,λλ==88..00 LLyy==66,,PABBCC,,λλ==88..00 ky/π (a) 3 (b)
=1
2
n0.04 )k ±1/2 2
∆ ( 0
S
0.02 S
1
1
(a)
0.00
Ly=4,PBC,λ=8.0 Ly=6,PBC,λ=8.0 0 0 2 8 1 0 0 2 8 1
0.06 Ly=4,ABC,λ=8.0 Ly=6,PBC,λ=5.3 1/ 7/ 1/ 7/
15
n0.04 ky/π (c) (d)
∆ 3
1
0.02 10 2/3
(b) ()k ±1/3 2
S ±
0.00 S 0
5
4 8 2 6 1
1/ ∞ 1/6 1/4 1/3 1/1
1/Lx 0 0
0 2 8 1 0 2 6 1
/ / / 1
1 7 1 /
FIG. 9. (Color online) Amplitude of the charge-density 13
k /π
stripes, ∆n, for width-4 and width-6 cylinders with periodic x
k /π
(PBC) and antiperiodic (ABC) transverse boundary condi- x
tions at n=0.875 with (a) U/t=4.0 and (b) U/t=8.0 as a
function of the inverse cylinder length, 1/Lx. FIG. 11. (Color online) Spin structure factor SS(k) at
n=0.875 and U/t=8.0 for (a) Ly =4, PBC, λ=8.0, (b)
Ly =4, ABC, λ=8.0, (c) Ly =6, PBC, λ=8.0, and (d)
Ly =6, PBC, λ=5.3.
k /π (a) (b)
y
1.5
2
1 whichcouldonlybedoneatthecostofslowingdownthe
1/2 algorithm significantly. Instead, we calculate the struc-
) ±
k 1.0 0 ture factor of the equal-time spin correlations,
(
S
S 1
0.5 SS(k)= 1 (cid:88)eik(r−r(cid:48)) S(r,r(cid:48)); (11)
N
rr(cid:48)
0.0 0 for finite cylinder length we expand in the particle-in-
a-box eigenmodes, as in Ref. [7]. Figures 10 and 11
ky/π (c) (d) show SS(k) for U/t=4.0 and U/t=8.0, respectively,
1 2 at n=0.875 (the same parameter sets as in Figs. 6
4 2/3 and7). AscanbeseeninFigs.10(a)–10(d)forU/t=4.0
)k ±1/3 and in Figs. 11(a)–11(c) for U/t=8.0, SS(k) is peaked
(S ± at k ( 7/8π,π) in all cases in which period λ=8.0
S 0 1 ≈ ±
2 stripes are present. Thus, the spin correlations are anti-
ferromagnetic with a π-phase shift every eighth site, as
expected. For the λ=5.3 stripes, Fig. 11(d), the peak
0 0 shifts to k ( 13/16π,π) which corresponds to antifer-
0 /2 /8 1 0 /2 /8 1 romagnetic≈sp±in correlations with a π-phase shift every
1 7 1 7
5.3sites,alsocompatiblewiththestripestructure. Aswe
k /π k /π
x x haveseenfortheamplitudeofthecharge-densitystripes,
the amplitudes of the peaks of the spin structure factor
FIG. 10. (Color online) Spin structure factor SS(k) for (a) also increase with increasing interaction strength.
and (b) 32×4 and (c) and (d) 32×6 Hubbard cylinders with In order to test for d -pairing-induced supercon-
x2−y2
periodic[(a)and(c)]andantiperiodic[(b)and(d)]boundary ductivity, we have calculated the equal-time pair-field
conditions at n=0.875 and U/t=4.0. correlation functions and have compared their decay to
that of the spin and charge correlations. To compen-
10
sate for the stripe structure of the ground state, we take
averaged absolute values; e.g., for the spin correlation 10−1 (a) U/t=8.0
functions we define U/t=4.0
1 (Lx−(cid:88)lx)/2+3 (l)x 10−3
S(lx)= S(x,x+lx) . (12) y
8 | | Dy 10 5
x=(Lx−lx)/2−4 −
The decay of the correlations with distance along the 10 7
−
cylinder is shown in Fig. 12 for width-4 cylinders. As
10−1 (b) U/t=8.0
10−1 (a) U/t=4.0
)x 10−3
l
(
10 3 y
− Dy 10 5
−
10 5
− Dyy(lx) 10−7
S(l )
x
10−7 C(lx) 0 10 20 l 30 40 50
x
10−1 (b) FIG. 13. (Color online) Equal-time pair-field correlation
functionsfor64×4Hubbardcylindersatn=0.875fordiffer-
10 3 entU/tfor(a)periodicand(b)antiperiodicboundarycondi-
− tions.
10 5
− D (l )
yy x
S(l ) Fig. 13(b), changing U/t has virtually no effect.
x
10−7 C(lx) The obvious question to be addressed is why the pair
correlations, and, indeed, also the charge and spin cor-
0 10 20 30 40 50 relations, are strong at short length scales, but have no
lx long-range or quasi-long-range (i.e., critical power-law)
order. Here we point out that the stripe state breaks
FIG. 12. (Color online) Equal-time pair-field Dyy(lx), spin the translational symmetry and is locked into the finite
S(lx),andchargeC(lx)correlationsasafunctionofthelongi- lattice,sothatthechargeordermanifestsitselfinavari-
tudinal distance lx for 64×4 Hubbard cylinders at n=0.875 ation of the local charge density. One then expects the
andU/t=4.0for(a)periodicand(b)antiperiodicboundary chargecorrelationwiththelocalordersubtractedouttobe
conditions. short-range, i.e., exponentiallydecaying. Asimilarargu-
ment holds for the spin correlations, which are locked to
can be seen by the approximately linear behavior of the the charge correlations (but have twice the wavelength
envelopes of the correlation functions on the semiloga- and a π-phase shift between stripes). The d pair
x2−y2
rithmic scale, all three correlation functions decay ex- correlationsdoshowstrongshort-rangecorrelations, but
ponentially at moderate to long distances for both peri- are also exponentially decaying, meaning that pairing is
odic and antiperiodic boundary conditions. For periodic not even present in the one-dimensional sense, i.e., that
boundary conditions, all three correlation functions de- the correlations decay as a power law. To explain this
cay approximately at the same rate, with the spin corre- behavior, wemakethreepoints: (i)Thelocked-incharge
lations having a larger absolute value. For antiperiodic ordercouldprecludeanyotherthanshort-rangeorderin
boundary conditions, the spin correlations are dominant all correlation channels. (ii) The pairing correlations are
for larger distance. Thus, we do not find that pairing measured perpendicularly to the stripes, i.e., in a direc-
correlations are long-range or, indeed, even dominant. tioninwhichchargetransportforstaticstripeswouldnot
Weinvestigatetheeffectofinteractionstrengthonthe be expected. (Note that the cylindrical geometry locks
strength of pairing correlations by comparing the pair in transverse stripes, and that measurement of pair cor-
correlationfunctionsfortheλ=8.0stripesforU/t=4.0 relations in the transverse direction over any significant
and U/t=8.0 in Fig. 13. For periodic boundary con- length is not possible due to limitations in the treatable
ditions, Fig. 13(a), there is a slight suppression of the system width.) (iii) For this doping of the system and
exponentially decaying correlation function as U/t is in- wavelengthofthestripes,thestripesarecompletelyfilled
creased, whereas for antiperiodic boundary conditions, with holes and thus insulating, so that quasi-long-range