Finite-temperature phase transitions in a two-dimensional boson Hubbard
  model
Finite-temperature phase transitions in a two-dimensional boson Hubbard model
Min-Chul Cha1 and Ji-Woo Lee2
Department of Applied Physics, Hanyang University, Ansan 426-791, Korea
Department of Physics, Myongji University, Yongin 449-728, Korea
We study finite-temperature phase transitions in a two-dimensional boson Hubbard model with
zero-point quantum fluctuations via Monte Carlo simulations of quantum rotor model, and construct
the corresponding phase diagram. Compressibility shows a thermally activated gapped behavior in
the insulating regime. Finite-size scaling of the superfluid stiffness clearly shows the nature of the
Kosterlitz-Thouless transition. The transition temperature, Tc, confirms a scaling relation Tc ∝ ρ
with x = 1.0. Some evidences of anomalous quantum behavior at low temperatures are presented.
PACS numbers: 73.43.Nq, 74.25.Dw, 05.30.Jp
Recently quantum phase transitions[1, 2] have drawn a
lot of attention in systems of interacting particles. Typ-
ically strong interactions suppress the itineracy of par-
ticles to induce a strongly correlated insulating phase,
whereas with weak interactions a conducting phase is
stable. The criticality of these zero-temperature phase
transitions can be investigated at low, but finite, tem-
peratures. How quantum fluctuations associated with a
quantum critical point(QCP) have influence on phases at
finite temperatures [3, 4, 5] is a theoretically interesting
and an experimentally relevant question.
At finite temperatures, it is expected that a quantum
phase transition turns into a classical one with the same
order parameter or disappears. Remnant quantum fluc-
tuations near a QCP may bring anomalous properties[3],
which can be captured by scaling relations, and lead
to crossover behaviors as temperature rises. Some pos-
sibilities such as reentrant behaviors due to the inter-
play of quantum and thermal fluctuations have been
proposed[6].
These issues can be clarified by direct investigations of
a generic quantum mechanical model. So far most of the
theoretical investigations heavily rely on the exact solu-
tion of the quantum Ising model, available strictly in one
dimension[4]. Interacting bosonic systems simulated via
Monte Carlo methods, not suffering from negative sign
problems, will be an ideal place to study these problems.
In previous works, a quantum XY model, equivalent to
hard-core bosons at half-filling, showed the Kosterlitz-
Thouless(KT) transition[7] at finite temperature in two
dimensions[8, 9]. In the model with nearest neighbor
repulsion, destruction of the solid order as well as the
superfluidity by thermal fluctuations was observed [10].
However, generic finite-temperature phase diagrams have
not been constructed.
In this work, we investigate the thermally driven phase
transitions of a two-dimensional quantum rotor model,
which is believed to share the same critical properties of
a soft-core generic boson Hubbard model[11], via Monte
Carlo simulations. The results are summarized in the
phase diagram as shown in Fig. 1. Finite-size scaling
properties of the superfluid stiffness confirm that the na-
ture of the classical phase transition associated with the
destruction of superfluidity is consistent with that of the
KT transition, and clearly support the scenario of the
universal jump at the critical point[12]. Finite tempera-
ture, T , sets the size in the temporal direction, leading to
a scaling behavior[4, 11] Tc ∝ ρ0x with x = 1.0, where Tc
is the transition temperature and ρ0 is the superfluid stiff-
ness at zero temperature. The compressibility diverges at
the transition. In the insulating regime at low tempera-
ture, thermally activated behavior of the compressibility
with a finite energy gap is observed. Some anomalous de-
pendence of energy and specific heat on T , possibly due
to quantum fluctuations, are observed for T < 0.25U .
The Hamiltonian of a boson Hubbard model reads
nj(nj − 1)− µ
nj − w
ibj + b
jbi),(1)
where bj(b
j) is the boson annihilation(creation) operator
0 0.01 0.02 0.03 0.04 0.05
µ=0.9
Normal fluid
Mott Insulator
Superfluid
KT transition
Gapped fluid
FIG. 1: (Color online) Phase diagram on the space of hop-
ping strength, t(=
n0(n0 + 1)w), and temperature, T , in
unit of U . The solid line denotes the classical phase transi-
tions, which terminates at a QCP at T = 0. The dotted line
represents crossover between gapped fluid and normal fluid.
http://arxiv.org/abs/0704.0396v1
at the j-th site, and nj is the number operator. U and
w stand for the strengths of the on-site repulsion and of
the nearest neighbor hopping, respectively, and µ is the
chemical potential.
It is convenient to put µ/U + 1/2 = n0 + n̄ with an
integer n0 and −1/2 < n̄ ≤ 1/2 so that n0 represents the
background number of bosons per site and n̄ is a charge
offset. When n̄ = 0, the density of bosons is fixed to
a commensurate filling across the transition. For non-
integer n̄, however, an integer filling in a Mott insulator
shifts to a non-integer filling in a compressible fluid. We
study the phase transition of the latter case in (2+1)-
dimensional L×L×Lτ square lattices, where L denotes
the size in a spatial dimension and Lτ in the temporal
dimension.
Since the phase transition of the model in Eq. (1)
is characterized by the establishment of phase coher-
ence, we may rewrite the Hamiltonian in terms of
the phase angle θj of bosons by replacing bj(b
j) =√
−iθj (
nj + 1e
iθj ) with nj =
. Under the as-
sumption that the nature of the transition is governed
only by the fluctuations of θj , not those of the hop-
ping strength, we replace nj → n0 so that bj(b†j) =√
−iθj (
n0 + 1e
iθj ). Then, the Hamiltonian is re-
duced to a quantum rotor model
nj(nj − 1)− µ
nj − 2t
cos(θi − θj),(2)
where t =
n0(n0 + 1)w. Here we take the number of
bosons nj ≥ 0.
Through a path integral mapping, we can construct
the corresponding classical action[14]
Jτr (J
r − 1)− ǫµJτr − ln IJxr (2ǫt)− ln IJyr (2ǫt)(3)
with the partition function
∇· ~J=0
{ ~Jr}
~J], (4)
where ǫ = β/Lτ is a lattice constant in the imaginary
time axis for an inverse temperature β, ~Jr is an integer
current at site r = (j, τ) with a spatial index j and a tem-
poral index τ , which is conserved at each site as denoted
by ∇ · ~J = 0, and Im(x) is the modified Bessel func-
tion given by the relation eK cos θ =
m=−∞ Im(K)e
In this work, we investigate the properties of the model
in Eq. 3 via Monte Carlo simulations using a recently
proposed worm algorithm [13]. In order to reduce the
systematic errors in discretizing the imaginary time axis,
we need to take ǫ
tU ≪ 1. We take Uǫ = 0.5 - 2 for
t ≪ U and set the energy unit U = 1.
The superfluid stiffness in a finite system is given by[14]
ρL = β
−1L2−d〈W 2x 〉, (5)
0.02 0.03 0.04 0.05
L=128
0 0.5 1
L exp(-b (t
) µ=0.9
=0.0409
b=1.85
0 10 20 30 40 50 60
0.005
0.015
L=128
0 0.1 0.2 0.3 0.4
L exp(-b (β
) µ=0.9
t=0.034
t=0.034
b=3.35
=28.8
FIG. 2: (Color online) Finite-size scaling behaviors of the su-
perfluid stiffness as a function of (a) hopping strength and
(b) temperature. For both cases, data collapsing onto a sin-
gle curve works fine in terms of the scaling parameter L/ξ as
shown in insets, consistent with the nature of the KT transi-
tion and the universal jump at the critical point.
where Wx = L
r and 〈...〉 denotes the averages
over the probabilites determined by the partition func-
tion of Eq. (4), and d is the spatial dimensionality. Sim-
ilarly the compressibility is
κ = βL−d[〈N2〉 − 〈N〉2], (6)
with N = L−1τ
r . The energy expectation is given
〈H〉 = L−1τ 〈
〉 , (7)
and the specific heat is CV = L
−d(∂〈H〉/∂T ).
We consider the case for µ = 0.9 so that n0 = 1 and
n̄ = 0.4. Figure 2 shows the finite-size scaling behav-
ior of the superfluid stiffness as a function of (a) t and
(b) β. Finite-size scaling properties of the transition can
be obtained by plotting the curves in terms of a scal-
ing variable L/ξ, where ξ is the correlation length. Here
we assume an essential singularity[15] ξ ∼ exp(bδ−1/2),
where δ = t − tc (or β − βc) is a tuning parameter and
b is a non-universal scaling factor. In terms of this scal-
ing variable, we obtain high-quality data collapsing onto
a single curve for different sizes, consistent with the na-
ture of the KT transition. The scaling behavior also sup-
ports the scenario of the universal jump of the superfluid
stiffness[12], (π/2)βcρ∞ = 1, at the critical point in the
thermodynamic limit. By extrapolating the single curves
to the critical point, we find that (π/2)βcρ∞ ≈ (a)1.01
and (b)1.06. These numbers are, however, sensitive to
fitting parameters b and tc(βc).
Figure 3a shows the behavior of the compressibility.
The finite-size scaling ansatz of the compressibility is
written in the form
κ = Lz−dX̃κ(L(t− tc)1/ν , β/Lz), (8)
where X̃κ is a dimensionless scaling function and z is the
dynamical critical exponent. For the generic superfluid-
insulator transition(GSIT), z = 2 is expected[11]. The
crossing behavior of the compressibility curves for differ-
ent sizes at the critical point t0c = 0.023±0.001, therefore,
represents the scaling properties near the QCP, where t0c
is the critical hopping strength at zero temperature. For
different values of µ, we have similar results with t0c just
shifted.
We find that the compressibility diverges at the transi-
tion. In the superfluid side, κ ∼ 1/(t− t0c). This strongly
supports that the longe-range density fluctuations drive
the transition. In the insulating side, the compressibility
has an activated form e−∆gap/T with a finite energy gap
∆gap. This dependence is shown in Fig. 3b for different
t, from which we can calculate ∆gap as shown in the in-
set. For small t, we need a large number of Monte Carlo
steps to obtain equilibrium and have bigger error bars in
determination of ∆gap. The gap vanishes around t = t
as expected.
Thus we have a so-called ’V-shaped’ phase diagram
(Fig.1). In the insulating regime, the Mott insulator
exists at T = 0, which turns into an activated gapped
fluid with a finite energy gap at low temperature. It
gradually disappears in a high-temperature normal fluid.
This crossover line can be specified by the condition
∆gap/T ≈ 1. The phase coherence in a superfluid at
T = 0 is destroyed by quantum fluctuations to form a
QCP or by thermal fluctuations at T > 0 to define clas-
sical phase transitions. The phase boundary in Fig. 1 is
obtained by tuning t for given T (black circles) as well
as by tuning β for a given t (red squares). Note that the
phase boundary follows a scaling relation Tc ∝ |t− t0c |zν ,
which implies that β determines the correlation length in
the temporal direction, where ν is the correlation length
0 0.01 0.02 0.03 0.04 0.05
L=12  Lτ=18
L=16  Lτ=32
L=20  Lτ=50
L=24  Lτ=72
L=28  Lτ=98
µ=0.9
0 50 100 150 200
1e-05
0.0001
0.001
0 0.01 0.02
µ=0.9
t=.005
t=.010
t=.015
t=.020
t=.021
t=.022
t=.023(b)
FIG. 3: (Color online) (a) Compressibility of the boson Hub-
bard model shows behavior of the GSIT with z = 2.0, diverg-
ing at the transition. (b) In the insulating regime, we have
thermally activated behaviors, κ ∼ e−∆gap/T , from which
∆gap can be evaluated. Inset: ∆gap as a function of t, van-
ishing at the QCP.
critical exponent. The boundary in Fig. 1 is consistent
with the expectation zν = 1[11] for the GSIT.
It is interesting to check the predicted scaling relation
[4, 11] Tc ∝ ρ0x in this model. Figure 4 shows that
the zero-temperature superfluid stiffness ρ0, denoted by
dotted line, which obtained via extrapolation of values at
T > 0, follows ρ0 ∝ |t− t0c |, implying that x = 1.0. It is
consistent with the hyperscaling argument[11] suggesting
x = z/(d+ z − 2).
We expect that this quantum criticality disappears
as temperature rises, which means quantum fluctuations
possibly leave some tracks in bulk properties at low tem-
peratures. Figure 5 shows the specific heat, CV , and the
energy expectation values, 〈H〉, as a function of T for dif-
ferent t. Sharp rises of CV in the conducting regime or
round up-rises in the insulating regime are followed by
indents, regions indicated by N, which apparently rep-
resent anomalous behavior due to quantum fluctuations
and disappear at high temperatures for T & 0.25. This
feature strongly suggests a crossover in normal fluid from
quantum mechanical to classical regime. Similarly the
curves of 〈H〉 show bumps, indicated by H, only in the
range where quantum critical fluctuations are expected
to have effects.
In summary, we have investigated the phase transi-
tions at finite temperature in a two-dimensional quan-
tum rotor model in which intrinsic zero-point fluctua-
tions are present. Finite-size scaling of the superfluid
stiffness shows an essential singularity of the KT phase
transition and the universal jump at the critical point.
The compressibility diverges at the transition. In the
insulating regime, the compressibility shows a thermally
activated behavior, κ ∼ e−∆gap/T , from which we can
successfully evaluate the gap. This indicates that the in-
sulating behavior at low temperature gradually crosses
over to the behavior of normal fluid as temperature in-
creases. The transition temperature Tc shows a scaling
behavior Tc ∝ |t − t0c |, showing that finite T limits the
length of quantum fluctuations in the temporal direction,
and a hyperscaling relation Tc ∝ ρ0. The behavior of the
specific heat and the energy suggests that, as tempera-
ture rises, quantum critical regime near a QCP crosses
over to classical regime.
MCC would like to thank Gerardo Ortiz for helpful
discussions and the hospitality of Department of Physics,
Indiana University, where parts of this work were carried
out. This work was supported by Korea Research Fund
grant No. R05-2004-000-11004-0.
0.02 0.025 0.03 0.035 0.04 0.045 0.05
β=200
β=400
µ=0.9
FIG. 4: (Color online) Superfluid stiffness for different β. As
β increases, the size dependence becomes smaller. This allows
us to extrapolate the curves to obtain zero-temperature su-
perfluid stiffness, ρ0, in the thermodynamic limit as denoted
by dotted line. It shows that ρ0 ∝ |t− t
c | with t
c ≈ 0.22.
[1] S. Sachdev, Quantum Phase Transitions (Cambridge
University Press, Cambridge, 1999).
[2] S. L. Sondhi, S. M. Girvin, J. P. Carini, and D. Sahar,
Rev. Mod. Phys. 69, 315 (1997).
[3] P. Coleman and A. J. Schofield, Nature (London) 433,
226 (2005).
[4] A. Kopp and S. Chakravarty, Nature Phys. 1, 53 (2005).
[5] S. Sachdev, Phys. Rev. B 55, 142 (1997).
[6] S. Kim and M. Y. Choi, Phys. Rev. B 41, 111 (1990).
[7] J. M. Kosterlitz and D. J. Thouless, J. Phys. C 6, 1181
(1973).
[8] H.-Q. Ding and M. S. Makivić, Phys. Rev. B, 42, R6827
(1990).
[9] K. Harada and N. Kawashima, J. Phys. Soc. Jpn. 67,
2768 (1998); A. W. Sandvik and C. J. Hamer, 60 6588
(1999).
[10] G. Schmid, S. Todo, M. Troyer, and A. Dorneich, Phys.
Rev. Lett., 88, 167208 (2002).
[11] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S.
Fisher, Phys. Rev. B 40, 546 (1989).
[12] D. R. Nelson and J. M. Kosterlitz, Phys. Rev. Lett., 39,
1201 (1977).
[13] F. Alet and E. S. Sørensen, Phys. Rev. E 67, 015701(R)
(2003); Phys. Rev. E 68, 026702 (2003).
[14] M. Wallin, E. S. Sørensen, S. M. Girvin, and A. P. Young,
Phys. Rev. B 49, 12115 (1994).
[15] J. M. Kosterlitz, J. Phys. C 7, 1046 (1974).
0 0.1 0.2 0.3 0.4 0.5
t = 0.005
t = 0.020
t = 0.040
t = 0.050
t = 0.070
t = 0.100
-0.85
t = .005
t = .010
t = .020
t = .025
t = .030
t = .040
0 0.1 0.2 0.3 0.4
t = .040
t = .050
t = .060
t = .070
t = .080
t = .090
t = .010
FIG. 5: (Color online) Specific heat, CV , as a function of T
for different t. Sharp rises in the conducting regime, signa-
ture of the superfluid transition, or round up-rise of CV in the
insulating regime are followed by indents which disappear in
high temperature region, T & 0.25. Insets: The curves of the
energy expectation values, 〈H〉, have bumps at low tempera-
tures possibly due to the effects of quantum fluctuations.
