Phase separation and flux quantization in the doped quantum dimer model
  on the square and triangular lattices
Phase separation and flux quantization in the doped quantum dimer model on the
square and triangular lattices
Arnaud Ralko,1 Frédéric Mila,2 and Didier Poilblanc1
1 Laboratoire de Physique Théorique, CNRS & Université Paul Sabatier, F-31062 Toulouse, France
2 Institute of Theoretical Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland
(Dated: April 5, 2007)
The doped two-dimensional quantum dimer model is investigated by numerical techniques on the
square and triangular lattices, with significantly different results. On the square lattice, at small
enough doping, there is always a phase separation between an insulating valence-bond solid and
a uniform superfluid phase, whereas on the triangular lattice, doping leads directly to a uniform
superfluid in a large portion of the RVB phase. Under an applied Aharonov-Bohm flux, the superfluid
exhibits quantization in terms of half-flux quanta, consistent with Q = 2e elementary charge quanta
in transport properties.
PACS numbers: 75.10.Jm, 05.50.+q, 05.30.-d
Understanding electron pairing in high temperature
superconductors is a major challenge in strongly corre-
lated systems. In his milestone paper, Anderson pro-
posed a simple connection between high temperature su-
perconductors and Mott insulators [1]. Electron pairs
”hidden” in the strongly correlated insulating parent
state as Valence Bond (VB) singlets lead, once fried to
move at finite doping, to a superconducting behavior.
A very good candidate of the insulating parent state is
the resonating VB state (RVB), a state with only expo-
nentially decaying correlations and no lattice symmetry
breaking. A simple realization of RVB has been pro-
posed by Rokhsar and Kivelson (RK) in the framework
of an effective quantum dimer model (QDM) with only
local processes and orthogonal dimer coverings [2]. Even
though the relevance of these models for the description
of SU(2) Heisenberg models is still debated, this approach
is expected to capture the physics of systems that nat-
urally possess singlet ground states (GS). For instance,
specific quantum dimer models have recently been de-
rived from a spin-orbital model describing LiNiO2 [3], or
from the trimerized kagome antiferromagnet [4]. In a re-
cent work, a family of doped QDMs (at T=0) generalizing
the so-called RK point of Ref.[2] has been constructed
and investigated[5], taking advantage of a mapping to
classical dimer models [6] that extends the mapping of
the RK model onto a classical model at infinite temper-
ature, with evidence of phase separation at low doping.
However, the soluble models of Ref.[5] are ’ad hoc’ con-
structions, and this call for the investigation of similars
issue in the context of more realistic models. In that re-
spect, a natural minimal model to describe the motion of
charge carriers in a sea of dimers is the two-dimensional
quantum hard-core dimer-gas Hamiltonian:
H = v
Nc|c〉〈c| − J
(c,c′)
|c′〉〈c| − t
(c,c′′)
|c′′〉〈c|
where the sum on (c) runs over all configurations of the
Hilbert space, Nc is the number of flippable plaquettes,
the sum on (c′, c) runs over all configurations |c〉 and |c′〉
that differ by a single plaquette dimer flip, and the sum
on (c′′, c) runs over all configurations |c〉 and |c′′〉 that
differ by a single hole hopping between nearest neigh-
bors (triangular) or (diagonal) next-nearest neighbors
(square). Throughout the energy scale is set by J = 1.
A schematic phase diagram for the two lattices is de-
picted in Fig.1 in the undoped case. Remarkably, these
lattices lead to quite different insulating states. Indeed,
an ordered plaquette phase appears on the square lattice
immediately away from the special RK point, whereas a
RVB liquid phase is present in the triangular lattice
PSfrag replacements
columnar
columnar staggered
staggered
liquid
plaquette
FIG. 1: (color online) Schematic phase diagrams for the tri-
angular and the square lattice.
In this Letter, we investigate in details the properties
of model (1) on the square and triangular lattices at fi-
nite doping. Building on the differences between the two
lattices in the undoped case, we investigate to which ex-
tent the properties of the doped system are governed by
the nature of the insulating parent state. This investi-
gation is based on exact Diagonalisations and extensive
Green’s Function Monte-Carlo (GFMC) simulations [7]
essentially free of the usual finite-size limitations [8].
Phase separation: At small t, it is expected that holes
http://arxiv.org/abs/0704.0715v2
experience an effective attractive potential. It is there-
fore natural to first address the issue of phase separation
(PS), i.e. the possibility for the system to spontaneously
undergo a macroscopic segregation into two phases with
different hole concentrations. We analyze the problem as
a function of the hopping parameter t and hole concen-
tration x = nh/N , where nh is the number of holes in the
system and N the number of sites. In order to perform
a Maxwell construction we define:
s(x) =
e(x)− e(0)
where e(x) is the energy per site at doping x. This quan-
tity corresponds to the slope of the line passing through
e(0) and e(x). If the system exhibits PS, the energy will
present a change of curvature implying s(x) to have a
minimum at a critical doping xc. The fact that the local
curvature of e(x) at x = 0 is negative then implies that
the two seperated phases will have x = xc and x = 0 (the
undoped insulator). In Fig.2, typical results are shown
for both square and triangular lattices and for different
sizes. Interestingly, PS appears in both cases, but with
0 0.05 0.1 0.15
-0.025
-0.02
3x6x6
3x8x8
3x10x10
0 0.05 0.1 0.15
-0.138
-0.136
-0.134
12x12
14x14
16x16
PSfrag replacements
Triangular
v=0.85, t=0.05
Square
v=0.90, t=0.10
FIG. 2: (color online) Slope of energy density (Eq.(2)) vs
doping for different sizes. (a) Triangular lattice. (b) Square
lattice.
noticeable differences. While for the square lattice (lower
panel) the critical hole concentration xc is roughly size
independent, there is a strong size dependence for the
triangular lattice (upper panel). This size effect can be
traced back to the nature of the parent undoped state.
On the square lattice, the crystalline phase (for v < 1) at
zero doping is very robust and for increasing size, its lo-
cal order changes only weakly. On the triangular lattice,
it has been shown that size effects play an important role
[8], especially in the RVB liquid phase for 0.8 . v ≤ 1.
Periodic boundary conditions (BC) tend to stabilize the
so-called
12 phase on small clusters, and clusters
with more than 192 sites are necessary to significantly re-
duce finite-size effects, in particular, as in Fig.2, close to
the transition point with the crystalline phase. Hence the
PS observed around x = 0.075 for the 3× 6× 6 cluster is
not representative of the thermodynamic limit.
To obtain the phase diagram in the (v, x) plane, we
have performed a systematic size-scaling analysis at fixed
t and for various v’s depicted in Fig.3. In agreement with
0.003 0.004 0.005 0.006 0.007 0.008
v=0.95
v=0.90
v=0.85
0 0.002 0.004 0.006 0.008 0.01
v=0.90
v=0.85
v=0.80PSfrag replacements
Triangular
Square
Triangular
t=0.05
Square
t=0.15
FIG. 3: (color online) Scaling of the critical doping xc de-
fined by Maxwell construction with the inverse total number
of sites. (a) Triangular lattice. (b) Square lattice.
the previous discussion, a significant size dependence is
only present for the triangular lattice, in which case PS
disappears for large clusters in the RVB phase in the
vicinity of the RK point [9]. In Fig.4, we report the ther-
modynamic limit of xc for the two lattices as a function
of v, and for different values of t. For the square lattice,
calculations have been done from the RK point down
to the expected phase transition between the plaquette
phase and the columnar phase, namely v ≃ 0.6 [10]. For
the triangular lattice, the range between the RK point
down to the RVB-
12 transition point at v ≃ 0.8
has been covered [8]. These results clearly demonstrate
the difference between the square and triangular lattices.
In the first case, as soon as v 6= 1, PS occurs for x < xc.
Moreover, upon decreasing v, crystalline order strength-
ens and, for fixed t, it is necessary to consider a higher
concentration of holes to reach a stable conducting phase.
Similarly, the bigger t, the lower xc. On the triangular
lattice, a finite size-scaling analysis shows that no phase
separation appears down to a critical value v ∼ 0.9, well
above the critical value v ∼ 0.8 below which plaquette
order sets in. Although numerical limitations prevent
computations for smaller v and t, our results up to the
3 × 12 × 12-site cluster provide clear evidence for a re-
gion of PS inside the RVB region, between v ∼ 0.8 and
v ∼ 0.9 [11].
Dimer ordering on the square lattice: Next we investi-
gate how dimer order, known to exist at x = 0, evolves
under finite doping. Two scenarii are a priori possible: i)
the dimer order vanishes in the stable conducting phase
0.6 0.7 0.8 0.9 1
t=0.05
t=0.10
t=0.15
t=0.20
PSfrag replacements
Square
Triangular
FIG. 4: (color online) Phase separation boundaries for the
square and triangular lattices in the thermodynamic limit.
The dashed lines correspond to the approximate location of
the phase transition between plaquette and columnar phases
[10] for the square lattice and between plaquette and RVB
phases [8, 12] for the triangular lattice.
immediately at xc; ii) dimer order survives above xc in
a narrow region of the conducting phase. To solve this
problem, we have calculated the squared order parameter
D2(~k) in the GS |Ψ0〉 defined by:
D2(~k) =
〈Ψ0|d(−~k)d(~k)|Ψ0〉
〈Ψ0|Ψ0〉
along the path of Fig.5(a) the first Brillouin zone of the
square lattice, where d(~k) is the Fourier transform of the
dimer operator defined on the horizontal bonds. Note
that this calculation has not been tried for the triangular
lattice since no Bragg peak is present in the RVB phase,
and the algorithm is losing efficiency for v . 0.8 [9]. In
the pure plaquette phase on the square lattice, a Bragg
peak develops at point ~kM = (π, 0), the middle of the side
of the BZ. We show in Fig.5(b) a typical result for the
squared order parameter on the 196-site cluster for differ-
ent values of x. Clearly, the Bragg peak disappears upon
doping. A finite size scaling of the order parameter can
be performed thanks to the linear behaviour at low con-
centration. Within our data, D2(~k = ~kM ) ≡ D2M (L, x)
behaves like aLx + bL in the linear region. In this
case, one can determine rather precisely x+∞ such that
D2M (+∞, x+∞) = 0 i.e. x+∞ is the concentration in
the thermodynamic limit where the Bragg peak vanishes.
By definition, x+∞ = −b+∞/a+∞ ≃ 0.05(8). This value,
and the linear behaviour in the thermodynamic limit, are
displayed in Fig.5(c). If we compare x+∞ to the corre-
sponding xc from Fig.4, we can conclude that the Bragg
peak indeed vanishes at xc ≃ 0.067(5) within error bars.
Note that numerical errors increase for larger clusters,
and we are not able to use the same analysis for clusters
larger than 256 sites (the results for the 18 × 18 cluster
have not been used). Although the determination of x+∞
0 0.03 0.06 0.09 0.12 0.15
0.005
0.015
12x12
14x14
16x16
18x18
TD limit
0.005
x=0.00
x=0.03
x=0.06
x=0.09
x=0.12
PSfrag replacements
v=0.90, t=0.15
FIG. 5: (color online) (a) First Brillouin zone. (b) Momen-
tum dependence of the squared order parameter for the 196-
site cluster (14x14). (c) Squared order parameter at the M
point as a function of the hole concentration x and for differ-
ent cluster sizes. The thermodynamic limit is depicted as a
dashed line, with the corresponding error bar (see main text).
is delicate, the linear behavior of D2 vs x is consistent
with the physical behavior expected for the binary sys-
tem of Fig.4. No dimer order is present above xc, showing
that the system is simply “conducting ”in this case, with
D2(~kM ) decreasing as N
−1 (critical behaviour).
Flux quantization: Finally, let us turn to a better char-
acterization of the “conducting phase”. Since holes are
bosonic one expects the conducting phase to be a su-
perfluid (through Bose condensation)[13]. However, ex-
tra complexity results from the fact that these bosons
move in a dimer environment. To investigate this is-
sue, we pierce the torus by an Aharonov-Bohm flux
of strength φ = ξφ0 with 0 ≤ ξ ≤ 1 and where
φ0 = hc/e is the magnetic flux quantum. To reduce
finite size effects, we also consider arbitrary BC in the
second direction (y). All this is implemented by the
Peierls substitution, changing the hole hopping term into
t′ = exp(±i2πξa/Lx± i2πηa/Ly)t, where the ± depends
on the directions ±x± y, while a is the lattice parameter
and Lx and Ly the linear sizes of the system. Obviously,
the whole spectrum should be periodic in ξ with period
1. We show in Fig.6 the spectrum of the 4 × 4 cluster
on the square lattice, with 4 holes. It turns out that (i)
the GS energy exhibits well-defined minima and (ii) is
rigorously periodic with period ξ = 1/2, which means
that there is flux quantization in units of half the flux
quantum (red curve). Property (i) is typical of a super-
fluid [14]: It is the precursor of a finite barrier in the ther-
modynamic limit. [15]. Property (ii) was suggested quite
some time ago in the context of a more general QDM
by Kivelson[13], who also predicted that, in the cylinder
geometry, one should be able to tunnel between the two
branches of Fig.6, thus lifting the degeneracy at the level
crossing. This degeneracy is not lifted in our case, neither
in the torus geometry, nor in the cylinder geometry, due
0 0.2 0.4 0.6 0.8 1
-0.54
-0.51
-0.48
-0.45
k(0,0)
k(π,0)
-0.48
-0.46
PSfrag replacements
FIG. 6: (color online) Energy spectrum vs (reduced)
Aharonov-Bohm flux ξ, for a 4 × 4 cluster with 4 holes at
v = 0.70 and t = 1.00. (a) Torus geometry. Arbitrary BC are
used in the transverse direction to obtain a continuous set of
momenta leading to a continuum (colored area). (b) The two
low-energy branches for the case of a cylinder (dashed line)
and including a small bond disorder (full red line).
to the translational symmetry, which puts the two states
that are degenerate at ξ = 1/4 into different symmetry
sectors. However, getting rid of the translational sym-
metry by changing the amplitude of a local dimer flip
indeed removes the degeneracy (upper panel of Fig.6),
leading to a detectable flux quantization in units of half
the flux quantum in an experiment in which the flux is
sweeped. Thus, in our model, the ground-state energy
has periodicity hc/2e, consistent with mobile elementary
particles of charge Q = 2e in the system. Unlike what
was recently found in a bosonic model with correlated
hopping[16], these particles are not boson pairs: ¿From
the bosonic point of view, it is the statistical flux of the
dimer background that leads to the half-flux quantiza-
tion. If dimers are interpreted as SU(2) electron singlets,
these singlets are the physical pairs that lead to half-flux
quantization. This scenario is fundamentally different
from the usual mechanism related to real space pairing
of the charge carriers found e.g. in the extended Hubbard
chain[17], in the 2-leg ladders [18] or, more generally, in
Luther-Emery liquids [19], as can be inferred from the
exact degeneracy between ξ = 0 and ξ = 0.5 for finite
systems in the present case, to be contrasted with the
significant finite-size effects of the other cases.
Summary and conclusions: The numerical investiga-
tion with Green’s function QuantumMonte Carlo and ex-
act diagonalizations of the doped two-dimensional quan-
tum hard-core dimer model on the square and triangular
lattices has led to a number of interesting conclusions re-
garding hole motion in a dimer background. Phase sep-
aration is often present at low doping, as suggested by
earlier investigations, but our results indicate that it is
related to the presence of valence bond order [20]: In the
RVB phase of the triangular lattice, PS only occurs close
to the plaquette phase, where short-range dimer correla-
tions are already strong enough. Close to the RK point,
doping the RVB phase leads directly to a superfluid phase
as shown from its response to an Aharonov-Bohm flux.
Moreover, we observed that the flux quantization is in
units of half a flux quantum, consistent with the idea
that the dimer background leads to effective particles of
charge 2e. All these results are in qualitative agreement
with the gauge theories of high Tc superconductivity in
strongly correlated systems [21].
We acknowledge useful discussions with Federico
Becca. This work was supported by the Swiss National
Fund, by MaNEP, and by the Agence Nationale de la
Recherche (France).
[1] P.W. Anderson, Science 235, 1196 (1987).
[2] D.S. Rokhsar and S.A. Kivelson, Phys. Rev. Lett. 61,
2376 (1988).
[3] F. Vernay, A. Ralko, F. Becca and F. Mila, Phys. Rev.
B 74, 054402 (2006).
[4] M. E. Zhitomirsky, Phys. Rev. B 71, 214413 (2005).
[5] D. Poilblanc, F. Alet, F. Becca, A. Ralko, F. Trousselet
and F. Mila, Phys. Rev. B 74, 014437 (2006).
[6] C. Castelnovo, C. Chamon, C. Mudry and P. Pujol, Ann.
Phys. 322, 903 (2007).
[7] N. Trivedi and D.M. Ceperley, Phys. Rev. B 41, 4552
(1990); M. Calandra and S. Sorella, Phys. Rev. B 57,
11446 (1998).
[8] A. Ralko, M. Ferrero, F. Becca, D. Ivanov, and F. Mila,
Phys. Rev. B 74, 134301 (2006) and references therein.
[9] For small values of t, the GFMC algorithm is no longer
ergodic due to hole localization.
[10] O.F. Syljuasen, Phys. Rev. B 73, 245105 (2006).
[11] Long-range Coulomb repulsion is expected to have a ma-
jor role in the PS region and might stabilize stripes.
[12] R. Moessner and S. L. Sondhi, Phys. Rev. B 63, 224401
(2001).
[13] S. Kivelson, Phys. Rev. B 39, 259 (1989).
[14] N. Byers and C.N. Yang, Phys. Rev. Lett. 7, 46 (1961).
[15] By contrast, for non-interacting fermions, signs of a flat
energy curve already appear on such small clusters pro-
vided one also uses arbitrary BC. Similar arguments were
used in D. Poilblanc, Phys. Rev. B 44, 9562 (1991) for
the 2D t-J model.
[16] R. Bendjama, B. Kumar, F. Mila, Phys. Rev. Lett. 95,
110406 (2005).
[17] K. Penc and F. Mila, Phys. Rev. B 49, 9670 (1994).
[18] C. A. Hayward, D. Poilblanc, R. M. Noack, D. J.
Scalapino and W. Hanke, Phys. Rev. Lett. 75, 926
(1995).
[19] A. Seidel and D. H. Lee, Phys. Rev. B 71, 045113 (2005).
[20] Our results give an explicit characterization of the
confinement-deconfinement discussed in O.F. Syljuasen,
Phys. Rev. B 71, 020401(R)(2005).
[21] T. Senthil and P.A. Lee, Phys. Rev. B 71, 174515 (2005).
and references therein.
