Thermodynamics of two-dimensional spin models with bimodal random-bond disorder

0.28 MB·English

by  Baoming Tang

Thermodynamics of two-dimensional spin models with bimodal random-bond disorder
Baoming Tang,1,2 Deepak Iyer,1 and Marcos Rigol1
1Department of Physics, The Pennsylvania State University, University Park, PA 16802, USA
2Department of Physics, Georgetown University, Washington, DC 20057, USA
We use numerical linked cluster expansions to study thermodynamic properties of the two-
dimensional classical Ising, quantum XY, and quantum Heisenberg models with bimodal random-
bond disorder on the square and honeycomb lattices. In all cases, the nearest-neighbor coupling
5 betweenthespinstakesvalues±J withequalprobability. Weobtainthedisorderaveraged(overall
1
disorder configurations) energy, entropy, specific heat, and uniform magnetic susceptibility in each
0
case. These results are compared with the corresponding ones in the clean models. Analytic ex-
2
pressionsareobtainedforlowordersintheexpansionofthesethermodynamicquantitiesininverse

n temperature.
u
J PACSnumbers: 61.43.-j,05.50.+q,75.10.Jm,05.10.-a

2

] I. INTRODUCTION Ourgoalinthisworkistousearecentlyintroducednu-
h merical linked cluster expansion (NLCE) for disordered
c
systems10 to study the thermodynamic properties of the
e
Solid-state materials deviate in various ways from the
m classical Ising (with S = 1/2), and quantum (spin-1/2)
periodic idealizations sometimes used to describe them
XY and Heisenberg models with bimodal random-bond
- theoretically. In crystals, for example, such deviations
t disorder on the square and honeycomb lattices. NL-
a can occur due to the presence of lattice distortions, im-
CEs allow us to obtain finite temperature properties of
t
s purityatoms(thatmayormaynotbemagnetic),andva- those models in the thermodynamic limit through the
. canciesinthelattice,collectivelytermedasdisorder. Dis-
t exact diagonalization of finite-size clusters. We specifi-
a order can significantly influence material properties. A
m callystudytheenergy,entropy,specificheat,anduniform
dramatic example for noninteracting electrons is Ander-
magnetic susceptibility (for the magnetization in the z-
- sonlocalization.1 Inspinsystems,the focus ofthis work,
d direction)asa functionoftemperature. Since anyglassy
quenched disorder in the spin interactions leads to frus-
n phase is only expected to emerge at zero temperature in
o tration and can generate spin glasses.2 A spin glass is a these models, if at all, we do not study the spin-glass
c remarkablestateofmatterwhere,looselyspeaking,spins order parameter. In future work, we will study quan-
[ are “frozen” in an irregular pattern, i.e., they display a tum quenches10–13 to examine the possibility of disorder
very slow dynamics under external driving. Although
2 driven localization in 2D.
v this phase does not exhibit spin order in a traditional Our results are briefly as follows: For clean systems,
0 sense (e.g., as a ferromagnet), it is still distinctly differ-
theymatchwellwithknownresultsforthesquarelattice,
9 ent from a paramagnet. Experimentally, spin glassesare
whilewereportadditionalresultsforthehoneycomblat-
9 generally associated with a cusp in the ac susceptibility
tice. For the disordered systems, we unveil some inter-
0
atacertaintemperatureabovewhichthesystembehaves
0 estingfeatures. IntheIsingmodel,theuniformsuscepti-
like a paramagnet.3 Whereas no standard order param-
. bility χ 1/T for all orders of the linked cluster expan-
1 eter captures a glassy transition, glass order parameters ∼
sion – we demonstrate this explicitly. The susceptibility
0 exist that do so (see, e.g., Refs. 2 and 4).
5 in the Heisenberg model also increases with decreasing
1 In two-dimensional (2D) lattices (the focus of this temperature (up to the lowest temperature we can ac-
: work), the effects of quenched disorder on classical cess). These two cases differ starkly from the XY model
v
spins have been extensively studied over the years.2 where the uniform susceptibility (for magnetization in
i
X However, the calculation of thermodynamic properties the z-direction, i.e., the same quantity calculated in the
r and correlation functions continues to be a computa- other two models) shows a plateau at low temperatures.
a
tional challenge.5,6 While it is generallybelieved that no At high temperatures, the clean and disordered models
spin-glass phase exists for nonzero temperatures,2 zero- behave identically as regards the energy, specific heat,
temperature phases continue to be debated for various and entropy up to third order in inverse temperature,
models ofinterest(see,e.g.,Refs.5and7). Forquantum and show identical susceptibilities up to second order –
spinmodels,thesignproblem8,9inquantumMonteCarlo weshowthis explicitlyinSec.IVviaahightemperature
simulations and the exponential growth of the Hilbert expansion.
space, relevant to full exact diagonalizationcalculations, The presentation is organized as follows. In section
represent an even greater challenge. Because of this, the II, we introduce the three models we study (the spin-
properties of disordered quantum spin systems, and of 1/2 Ising, XY, and Heisenberg) and summarize some of
quantumspinglassesinparticular,haveremainedessen- theirknownpropertiesinthesquareandhoneycomblat-
tially unexplored. The existing literature in the subject tices. Section III briefly describes NLCEs for systems
has almost exclusively dealt with classical models. with disorder. Numericalresults for the latter, andtheir
2
comparison with those for clean systems, are presented The Ising model has a discrete Z symmetry, i.e., the
2
in Sec. IV. We conclude with a brief summary in Sec. V. transformation Sz Sz for all j leaves the Hamil-
j → − j
tonian invariant. This symmetry can be spontaneously
broken at sufficiently low temperatures to create an or-
II. MODELS dered phase. For Sz = 1/2, the critical temperature is
i ±
given by T /J = 1/2log(1+√2) 0.57 for the square
c
| | ≈
We areinterestedinthermodynamicpropertiesofvar- lattice, and T /J = 1/2log[(√3+1)/(√3 1)] 0.38
c
| | − ≈
ious spin-1/2 models on the square and honeycomb lat- for the honeycomb lattice, and is the same for the ferro-
tices (see Fig. 1). In the absence of disorder, and for and antiferromagnetic cases.14 The latter is because, for
nearest neighbor interactions, those models do not ex- bipartite lattices, a unitary transformation relates both
hibit frustration on either lattice, which are both bipar- models. This can be easily seen by rewriting the Hamil-
tite. While the thermodynamic properties of the vari- tonian in Eq. (1) as H = J Sz Sz , where A
Ising i i i,A i,B
ousmodelsstudiedherearequalitativelysimilaronboth and B are the two sub-lattice indices. One can then go
lattice geometries, there are significant quantitative dif- fromJ J ,i.e.,betweenthePferro-andantiferromag-
i i
ferences, e.g., in the critical temperatures for the onset netic mo→de−ls, via the transformation: Sz Sz and
of quasi-long-range order.14 These differences have their Sz Sz . i,A →− i,A
i,B → i,B
origin in the different coordination number in both lat-
The specific heat of the clean model diverges at the
tices, with the honeycomb lattice having the smallest
critical temperature for both the square and the hon-
one. Hence, not surprisingly, for the spin-1/2 antifer-
eycomb lattices. Equivalently, the derivative of the en-
romagnetic Heisenberg model, the staggered magnetiza-
ergy diverges at the critical temperature but the en-
tion is significantly suppressed on the honeycomblattice
as compared with the square lattice.15 Disorder, on the ergy remains finite throughout. For the antiferromag-
netic model, the susceptibility is known to be finite ev-
otherhand,leadstofrustrationonbothlattices,andtoa
erywhere, with an infinite slope at the critical temper-
qualitativechangeoftheintermediatetolowtemperature
ature. The maximum value of the susceptibility occurs
propertieswithrespecttothecleansystems. Frustration
at T = 1.537T and T = 1.688T for the square and
can be easily identified by trying to assign spins to the m c m c
honeycomb models respectively,14,19 i.e., above the crit-
varioussitesinFig.1tominimizetheenergy. Ifonetakes
ical temperature. Our results for the clean systems are
theIsingHamiltoniandiscussedbelow,onefindsthatfor
consistentwiththese. However,wecannotstudythecrit-
the overwhelming majority of disorder realizations there
ical phase or the properties of the system very close to
is no single spin configurationthat minimizes the energy
criticality (see Figs. 3 and 4).
on all bonds.
Weshouldstressthatbothforthecleananddisordered
TheIsingmodelwithbimodaldisorderhasbeenexten-
cases,weexpectquantumfluctuationstostronglymodify sively studied in the past.4,5,7,20–27 It is reasonably well
the results for the spin-1/2 XY and Heisenberg models established that no glassy phase exists for T > 0.22,23
from their classical counterparts (see, e.g., Ref. 16 for
It has, however, been established that a glassy phase
examples of the effects of quantum fluctuations in frus- appears at zero temperature in this model.7 At finite
trated spin-1/2 XY and Heisenberg models on the hon-
temperature, our results for this case are described in
eycomb lattice). Our focus in this work is on systems
Sec. IVA.
with bimodal random-bond disorder, where the nearest-
neighborcouplingbetweenthespinstakesvalues J with
±
equal probability.
A. Ising model
The Hamiltonian for the spin-1/2 Ising model can be
written as
H = J SzSz (1)
Ising ij i j
Xhiji
where Sz = 1/2 is the spin at site-i, and the sum is
i ±
overnearestneighbors. Inthe absenceofdisorder(J =
ij
J for all i,j), Eq. (1) has served as the quintessential FIG. 1. (Color online) Schematic of lattice models, square
model for magnetism and was solved exactly on a 2D (left) and honeycomb (right), with bond disorder considered
squarelatticebyOnsager.17Inthepresenceofcontinuous in this work. The red and blue bonds represent Jij = ±J.
disorder, the Hamiltonian in Eq. (1), commonly known Black dots at vertices represent the spins. It is easy to see
as the Edwards-Anderson model,18 has also become a thatthebonddisordercausesfrustrationbytryingtoarrange
widely studied model for spin-glasses. thespins to minimize energy.
3
B. XY Model non-Abelian, as opposed to the XY model which has
an Abelian symmetry group, U(1) [O(2)] for the quan-
The spin-1/2 XY model can be written as tum (classical) cases. Vortices or point-defects, which
are responsible for the BKT transition,31 occur only in
Hˆ = J (SˆxSˆx+SˆySˆy) (2) thelattercase.36 Furthermore,in2D,O(N)modelswith
XY ij i j i j
N 3 or SU(N) models with N 2, i.e., with non-
Xhiji Abe≥liansymmetrygroups,canbesho≥wnviaperturbation
where Sˆx,y are spin operators at site i, proportional to theory to be asymptotically free, which for spin mod-
i els translates to a renormalization group flow towards
the Pauli matrices. We consider only the isotropic case,
paramagnetism.37,38
where the modelhas a continuousU(1) rotationsymme-
Theresultsforthe cleansystempresentedhereforthe
tryintheplane. Thepresenceofacontinuoussymmetry
precludes, via the Mermin-Wagner theorem,28 a finite- square lattice are nearly identical to the spin-1/2 results
presented in Ref. 39, which also compares with other
temperature phase-transition involving the breaking of
known results.
this continuous symmetry from occurring. For d 2 di-
≤
mensions,the fluctuations inanyputative orderedphase
appearing from breaking a continuous symmetry grow
III. NUMERICAL LINKED CLUSTER
with system size for finite temperatures, destroying any
EXPANSIONS
order (see e.g., Ref. 29). Hence, this model has an or-
dered phase only at T =0.
However, in 2D there can still be a finite temperature Numerical linked cluster expansions (NLCEs) are a
Berezinskii-Kosterlitz-Thouless (BKT) transition30,31 computational technique that can be used to calcu-
below which the system exhibits quasi-long-range spin late extensive properties (per lattice site) of translation-
order. The critical temperature for the BKT transi- ally invariant lattice systems. NLCEs, which are based
tion in the spin-1/2 XY model in the square lattice is on linked cluster expansions,40–42 were introduced in
T /J 0.34.32,33 We shouldstressthatboth classical34 Ref.43, whereitwasshownthatthe resultsobtainedfor
c
and|q|u≈antum models with continuous symmetries in 2D thermodynamic properties were exact in the thermody-
canexhibit this kind of behavior.35 We should alsomen- namic limit for systems with a finite correlation length.
tion that most studies in the literature report results for Furthermore, results could be obtained at significantly
the ferromagnetic XY model (J < 1). However, like for lowertemperaturesascomparedtohigh-temperatureex-
the Ising model, in the square and honeycomb lattice pansionsformodelsthatdeveloplong-rangeorderatzero
a unitary transformation relates the ferro- and antifer- temperature. In several subsequent works, NLCEs have
romagnetic models and the critical temperature is the beenshowntobeapowerfulcomputationaltechniquefor
same in both. Our calculations for the susceptibility of determining not only thermodynamic properties of a va-
thecleancaseonthesquarelatticeconvergedowntotem- rietyoflatticemodels,44–49butalsoforstudyingthermal-
peratures of T/J 0.4, which is compatible with the ization(or the lackthereof) atlong times after a quench
onset of quasi-lo|ng|-≈rangeorder for T /J .0.34.32,33 in isolated quantum systems.10–13 For completeness, we
c
| | provide a brief description of NLCEs. Details of how to
implement them can be found in Ref. 50.
C. Heisenberg model In NLCEs, the expectation value of an extensive ob-
servable,per-site, inatranslationallyinvariantsystem
O
can be calculated as a sum over contributions from all
The spin-1/2 Heisenberg model, also known as the
clusters c of different sizes that can be embedded on the
XXX model, can be written as
lattice
Hˆ = J Sˆ Sˆ (3)
Heis ij i· j = M(c) WO(c), (4)
Xhiji O c ×
X
where Sˆi =(Sˆix,Sˆiy,Sˆiz), andSˆix,y,z arespinoperatorsat whereM(c)isacombinatorialfactorequaltothenumber
site i, proportional to the Pauli matrices. The Heisen- of ways that a particular cluster c can be embedded on
berg model has an SU(2) symmetry, the highest of the that lattice. W (c) is the weight of cluster c for the
O
three models considered in this work. The ground state given observable, which is calculated via the inclusion-
in the clean case (Jij = J) is an ordered ferromagnet exclusion principle
orantiferromagnetdependingonthesignofthecoupling
constantJ. LikefortheXYmodel,long-rangeorderonly W (c)= (c) W (s). (5)
O O
occurs at zero temperature. However, in contrast to the O −
s⊂c
XY model, the 2D Heisenberg model does not develop X
quasi-long-rangeorderatfinite temperature. Thisis due (c) is the expectation value of the observable on the
O O
to the fact that the internal symmetry group, SU(2) specific cluster c. Within NLCEs, (c) is calculated us-
O
[O(3)] for the quantum (classical) Heisenberg model, is ing full exact diagonalization. The expansion is carried
4
disorder average of the weights is in turn given by
W (c)= (c) W (s). (7)
O O
O −
s⊂c
X
In other words, the disorder average can be carried out
orderbyorderforeachobservable. Thecalculationsthen
proceedas for the translationally invariantsystem if one
replaces (c) by (c).
O O
Thecomputationsinthepresenceofdisorderaremuch
more challenging than for translationally invariant sys-
tems because of the additional average over all possible
disorder realizations. For example, the largest clusters
FIG. 2. (Color online) Clusters up to thefourth order in the
we consider here for the quantum models in the square
sitebasedNLCEonthesquarelattice. Thetwo3-siteclusters
lattice have 13 sites. At this order, there are a total of
have the same Hamiltonian. At fourth order, in addition to
three clusters with identical Hamiltonian, two topologically 1.9 106 connected clusters, of which 5,450 are topo-
new clusters appear – the closed loop and the “⊥”. Each ∼logical×ly distinct.44 Each of these has to be fully diago-
topologically new cluster is diagonalized separately. nalized for the 2ℓ disorder configurations corresponding
to the ℓ bonds in the cluster. This has to, of course, be
carried out for all lower orders as well, each with a dif-
outorder by order,i.e., by firstconsideringclusters with ferent set of topologically distinct clusters and disorder
one-site,then two-sites,andsoon. Beyondthe bare sum configurations. For the clean systems, we report results
in Eq. (4), several resummation schemes exist that ac- for cluster with up to 15 sites for the square lattice. In
celeratetheconvergenceoftheexpansion.44 Herewewill thatcaseone hasto diagonalize42,192topologicallydis-
report results from Wynn and Euler resummation tech- tinct clusters with 15 sites.
niqueswhenevertheyallowustoextendtheconvergence NLCE calculations fail to converge when correlations
of the results to lower temperatures.44 in the thermodynamic limit extend beyond the largest
clusters considered. Therefore, NLCEs cannot be used
Examples of clusters up to fourth order in the site-
to calculate observables in phases with long-range order
based NLCE used here on the square lattice are shown
inFig.2. At thirdorder,althoughtherearetwogeomet- unless one tailors the expansion to account for those.51
rically different clusters, they are topologically identical. Since disorder usually shortens correlations at low tem-
peratures, NLCEs are particularly useful to study quan-
They have the same Hamiltonian for the models with
tum disordered systems, despite the increase of compu-
nearest interactions we consider here. At fourth order,
tational cost because of the disorder average. It will be-
there are three clusters (including the one with the four
comeapparentwhenwediscussourresultsforthevarious
sites on a line) that again have the same Hamiltonian
for the models considered here. However, two topologi- thermodynamic properties of interest in this work, that
cally new clusters appear, namely, the closed loop and NLCEsconvergetolowertemperaturesindisorderedsys-
temswhencomparedtocleansystems. Asmentionedbe-
the “ ”. They have to be individually diagonalized.
⊥ fore, quantum Monte Carlo simulations in the presence
Fromthefourthorderandbeyondthenumberofdistinct
of disorder are severely limited by the sign problem.
topologicalclustersincreasesrapidly(exponentiallywith
thenumberofsites),makingthecalculationsincreasingly
costly. References 44 and 49 provide details on the var-
IV. RESULTS
ious topological clusters on the square and honeycomb
lattices, respectively, as well as the number of such clus-
ters as a function of the order of the expansion. In this section, we discuss the results of our NLCE
basedstudy ofthe three spinmodels describedin Sec. II
Recently,inRef.10,itwasshownthatNLCEscanalso
on the square and honeycomb lattices. In what follows,
be used to study systems with disorder. As described
we set J = 1 as the energy scale. For each model and
above, NLCEs can only be used for translationally in-
lattice geometry, we report the energy (E), entropy (S),
variant systems. A priori, disorder breaks translational
specificheat(C ),anduniformsusceptibilityforthemag-
invariance. However, we are only interested in disorder v
netization along the z-axis(χ) as a function of tempera-
averaged physical quantities. If we take a disorder av-
ture. These quantities are defined as
erage over all possible disorder configurations in models
with bimodal disorder, we restore translational invari- Hˆ Hˆ2 Hˆ 2
ance and the equivalent of Eq. (4) reads E = h i, C = h i−h i
N v NT2
(8)
logZ E (Sˆz)2 Sˆz 2
= M(c) WO(c) (6) S = + , χ= h i−h i
O × N T NT
c
X
where the overline denotes a disorder average, the an-
where the overline denotes the disorder average. The gle brackets denote the thermal expectation value in the
5
grand-canonicalensemble(atzerochemicalpotential),N Letustreattheabovetermsonebyone. Inthefirstterm,
isthenumberofsites,andT isthetemperature. Asmen- for ij = kl , the trace is identically zero as shown
h i 6 h i
tionedearlier,the disorderaverageis carriedoutoverall above. For the case when say i = l, but j = k, we
6
disorderconfigurationsateachorderintheNLCE.Inall effectively haveanew Hamiltonianfor three neighboring
cases, the disorder averaged results are compared with spins. It is easy to verify explicitly that this trace also
the ones in clean systems. vanishes. The only possibility left for a nonvanishing
For all observables, we report bare NLCE results for contribution is ij = kl . The trace of second term
h i h i
thehighesttwoordersofoursite-basedNLCEexpansion, in the brackets is nonzero only for i = j. In the third
which are determined by the number of sites l in the term, againfor k =i or j, the trace is zero. For, say j =
6
largest clusters studied. Namely, we report the results k, considering only the diagonal elements of the matrix,
from Eq. (6) when the contributions of all clusters with SˆzSˆzSˆz Sˆz, we see that the trace vanishes. Taking
i j j ∝ i
upto l sitesareadded, wherel takesthe twolargestval- these into account, and writing logZ to O(β2), we have
ues in our calculations for each model in each lattice ge-
ometry. Wealsoreportresultsusingtwodifferentresum- logZ β2 Nh2
mation schemes, indicated as Wynnn and Eulern. The N =log2+ 2 2NN Tr Ji2jHˆi2j + 4 (10)
" #
subscript n denotes the order of the resummation pro- · Xhiji
cess (see Ref. 44 for details). The resummation schemes
allow us to access significantly lower temperatures than For both, the cleansystem andthe systemwith bimodal
the bare results in some cases as indicated below. disorder, Ji2j =J2. We therefore get,
Before discussing eachmodel andobservable in detail,
we review a few general observations for completeness logZ =log2+ β2 J2Tr Hˆ2+ h2 +O(β3) (11)
and pedagogy. In all models and observables discussed N 4 2 2 2
(cid:18) (cid:19)
here, the numerical results at intermediate to high tem-
peraturesinthepresenceofdisorderareclosetothoseof where Hˆ is the Hamiltonian for a two-site system and
2
the correspondingcleansystem. Whereas this is obvious the subscript “2” on the trace indicates a trace over the
for temperatures so high that the first order correction Fock space of a two-site system. The disorder average
to the infinite temperature result is negligibly small, we does not change the above expression, and therefore to
notice from our results in Figs. 3–8 that the observables this order, the clean and the disordered systems behave
in clean and disordered models are indistinguishable for identically. The energy, for instance, is given to this or-
temperaturesaslowasT =2(barringthesusceptibility). der by E = (∂logZ/∂β)/N = βJ2Tr (Hˆ2)/2, the
− − 2 2
In order to show why this is so, we expand the par- specific heat is C = β2J2Tr (Hˆ2)/2, and the uniform
v 2 2
tition function for small inverse temperature β 1/T, susceptibility is χ=β/4.
≡
Z =Tr(e−βHˆ) Tr(1 βHˆ +β2Hˆ2/2+...). The mod- In fact, it is straightforward to check that the par-
≈ −
els we consider have only nearest-neighborcoupling, i.e., tition function in the clean and disordered systems re-
Hˆ = [J Hˆ +h(Sˆz +Sˆz)/2], where we have in- mains the same at O(β3), except for terms proportional
hiji ij ij i j
cluded a magnetic field h as a source to calculate the to h2. In other words, the energy, entropy, and specific
P
uniform susceptibility in the z-direction. The most gen- heat are the same for clean and disordered systems up
eraltwo-siteHamiltonianthatdescribesallmodelsofin- to third order in β, but the uniform susceptibility devi-
terest here is given by H =γ(SˆxSˆx+SˆySˆy)+∆SˆzSˆz, ates. The following expression can be derived along the
ij i j i j i j
which becomes the Ising model for γ = 0,∆ = 1, the linesofEqs.(9)–(11)forthethirdordercorrectiontothe
XY model for γ = 1,∆ = 0, and the Heisenberg model partition function:
for γ = ∆ = 1. With this in mind, to first order in
β, the high temperature expansion for Z can be writ- logZ β2 h2
ten as Z = 2N − β hijiTr[JijHˆij + h(Sˆiz + Sˆjz)/2], N =log2+ 4 (cid:18)J2Tr2Hˆ22+ 2 (cid:19)−
where N is the number of lattice sites. Note first that
Tr(Sˆz)=0(thePaulimPatricesaretraceless),andsecond, β3h2J Tr Hˆ SˆzSˆz. (12)
i − 12 ij 2 ij i j
that Tr(Hˆ ) = 0, so that the linear correction vanishes.
ij
Therefore, to first order in β, the clean and the disor- There is no sum over i,j, which represent the two sites
dered system have the same partition function. This is in a 2-site system. For the cleanmodelone just needs to
true regardless of the type of disorder. replaceJ withJ intheaboveformula,whileforthedis-
ij
To second order, after expanding Hˆ2, we have orderedone,thedisorderaverageproducestwotermsfor
J (which cancel each other, implying that disorder ex-
β2 ±tends the paramagnetic behavior in the susceptibility to
Z =2N + Tr J J Hˆ Hˆ +
2 ij kl ij kl lower temperatures). One can the see that for h=0, all
"
hijXi,hkli thermodynamic quantities studied here, except the sus-
ceptibility,areidenticalinthecleananddisorderedcases
+h2 SˆizSˆjz +h JijHˆijSˆkz (9) up to O(β3). One canfurther verify that this changesat
#
Xi,j hXiji,k fourthorder,wheredifferencesemergebetweencleanand
6
disordered systems. An example of a term that makes a sion for all temperatures.
difference at fourth order is a square loopwith four sites Figures 3(a) and 4(a) show that, as mentioned before,
and four bonds (see Fig. 2). Even in the Ising case, the a generic feature in the presence of disorder is that the
Hamiltonian H12H23H34H41 =1/64 has a nonzero trace energy tends to plateau at the lowest temperatures ac-
with a product of four different Jij. cessibletous. Inthatregime,theentropyissignificantly
Weshouldstressthat,inallmodelsstudiedhereinthe higher than in the clean systems [Figs. 3(b) and 4(b)].
presence of disorder, we find that there is a significant Distinct to the Ising models, the sharp divergence in
amount of residual entropy (when comparing with the the specific heat in the clean case [see Figs. 3(c) and
clean systems) at the lowest temperatures we are able 4(c)],whichindicatesthephasetransition,isreplacedby
to access with NLCEs. This is a clear indication of the whatappearstobesmoothpeakinthepresenceofdisor-
lack oforderat those temperatures. The behavior ofthe der. Themaximumofthatpeakappearsattemperatures
entropy,coupledwithasaturationoftheenergyobserved lowerthan the criticaltemperature in the cleancase. At
atthelowesttemperaturesaccessibletous,confirmsthat highertemperatures,Eq.(11)givesresultsthatagreefor
there are many energy levels close to each other at low E, C , and S down to T 1. We note that NLCE re-
v
≈
energies. This is the hallmark of frustration. sults for the disordered model are well converged down
to T 0.2 to 0.3, while the NLCE results for the clean
≈
modelconvergetotemperaturesthatareclosetoT ,and
c
A. Ising model agree with the analytic results in the disordered phase.
ResultsfortheuniformsusceptibilityinFigs.3(d)and
Figures 3 and4 show the energy (a), entropy(b), spe- 4(d) show that this quantity behaves very differently in
cific heat (c), and uniform susceptibility (d) for the dis- the cleanand disorderedsystems. In the disorderedcase
ordered spin-1/2 Ising model on the square and honey- it exhibits a 1/T behavior at all temperatures, both on
comb lattices, respectively. For the square lattice, we the square and honeycomb lattices. An order by order
also plot the exact results for E, S, and C for the clean linkedclusterexpansionrevealsthatthe onlynonvanish-
v
model.17,52 Several approximate analytic estimates exist ing contribution to the susceptibility in the disordered
for the susceptibility, but there is no closed form expres- case comes from the single-site system, and is trivially
0 0
(a) E (b) S (a) E (b) S
0.6 0.6
-0.2
0.4 -0.2
Disordered Clean Disordered 0.4
l = 13 l = 14 l = 14 Clean
-0.4 l = 14 l = 15 l = 15 l = 16
Wynn Wynn 0.2 Wynn l = 17
Euler6 Euler7 -0.4 Euler6 0.2
11 12 12
2 0
(c) Cν (d) χ χ = 0.25/T (c) Cν (d) χ χ = 0.25/T
1.5 Clean
Exact 0.8
0.1
1 0.1
0.4
0.5
0 0
0.1 1 0.1 1 10 0.1 1 0.1 1 10
T T T T
FIG. 3. (Color online) Spin-1/2 Ising model on the FIG.4. (Coloronline)Spin-1/2 Ising model on the hon-
squarelattice. Panels(a)–(d)showtheenergy,entropy,spe- eycomb lattice. Panels (a)–(d) show the energy, entropy,
cificheat,anduniformsusceptibilityvsT,respectively. Solid specific heat, and uniform susceptibility vs T, respectively.
lines depict disorder averaged quantities, while dashed lines Solid lines depict disorder averaged quantities, while dashed
depict results for the clean system. Thin lines report bare linesdepictresultsforthecleansystem. Thinlinesreportbare
results for thelast two orders of the NLCE, while thick lines results for the last two orders of theNLCE, while thick lines
reporttheresultsoftworesummationtechniques. Athincon- reporttheresultsoftworesummationtechniques. Athincon-
tinuouslinefollowing theresultsoftheresummationsreports tinuouslinefollowingtheresultsoftheresummationsreports
results for a lower order of the same resummation technique results for a lower order of the same resummation technique
and is used to gauge their stability. The dotted vertical line andisusedtogaugetheirstability. Resummation resultsare
marks the position of the phase transition. The dashed dot- not presented for the clean case as they do not extend the
ted line shows exact analytic results for the clean system in convergence to lower temperatures. The dotted vertical line
thethermodynamic limit. marks theposition of thephase transition.
7
proportionalto 1/T. All higher order contributions van- disorder averaged logZ for the clusters with one, two,
ish. We show this for the first few orders of the NLCE and three sites shown in Fig. 2.
on the square lattice. Below are the expressions for the
logZ(1) =log e−β2h +eβ2h
logZ(2) = 1 l(cid:16)og 2e−β4J +(cid:17)eβ(J+44h) +eβ(J−44h) +β β
2 →−
(13)
logZ(3) = 1hlog (cid:16)eβ2h +e3β2h +eβ(J2+h) +eβ(J2−h(cid:17)) +β βi +
2 →−
1(cid:16) βh −βh −β(J+3h) −β(J−3h(cid:17)) β(J+h) β(J−h)
+ log 2e 2 +2e 2 +e 2 +e 2 +e 2 +e 2 +β β
4 →−
h (cid:16) (cid:17) i
The uniform susceptibility can be obtained from χ = multiplicities, etc.):
β−1(∂2logZ)/(∂h)2 evaluated at h = 0. We get for the
β
above three orders, W(1) =χ(1) =
χ 4
W(2) =χ(2) 2W(1) =0 (15)
β 2β 3β χ − χ
χ(1) = , χ(2) = , χ(3) = . (14) W(3) =χ(3) 2W(2) 3W(1) =0.
4 4 4 χ − χ − χ
Indeed, only the single-site cluster contributes. One can
Already, one can see that no new contributions appear checkthisathigherorders,andforthehoneycomblattice
at the higher orders. To confirm this, we calculate the as well. This is not the case for the XY and Heisenberg
weightsofthethreeclusters(seeRef.44fordetailsabout models discussed below.
Earlierstudiesfor the Isingmodelwithbimodaldisor-
0 0
(a) E (b) S (a) E (b) S
0.6 0.6
-0.2
0.4 -0.2 0.4
Disordered Clean Disordered Clean
-0.4 l = 12 l = 14 l = 13 l = 15
l = 13 l = 15 l = 14 l = 16
Wynn Wynn 0.2 Wynn Wynn 0.2
-0.6 Euler6 Euler7 -0.4 Euler6 Euler7
11 12 11 12
0 0
(c) C (d) χ (c) C (d) χ
0.6 ν χ = 0.25/T ν χ = 0.25/T
0.4
0.4
0.1
0.2 0.2 0.1
0 0
0.1 1 10 0.1 1 10 0.1 1 0.1 1 10
T T T T
FIG.5. (Coloronline)Spin-1/2XYmodelonthesquare FIG. 6. (Color online) Spin-1/2 XY model on the hon-
lattice. Panels (a)–(d) show the energy, entropy, specific eycomb lattice. Panels (a)–(d) show the energy, entropy,
heat,anduniformsusceptibilityvsT,respectively. Solidlines specific heat, and uniform susceptibility vs T, respectively.
depictdisorderaveragedquantities,whiledashedlinesdepict Solid lines depict disorder averaged quantities, while dashed
results for the clean system. Thin lines report bare results linesdepictresultsforthecleansystem. Thinlinesreportbare
for the last two orders of the NLCE, while thick lines report results for the last two orders of theNLCE, while thick lines
theresultsoftworesummationtechniques. Athincontinuous reporttheresultsoftworesummationtechniques. Athincon-
line following the results of the resummations reports results tinuouslinefollowingtheresultsoftheresummationsreports
for a lower order of the same resummation technique and is results for a lower order of the same resummation technique
used to gauge their stability. The dotted vertical line marks and is used to gauge their stability. Note that the results
theposition of the phasetransition. converge tosimilar temperatures as in thesquare lattice.
8
deronthesquarelatticefoundalow-temperaturescaling temperature for low temperatures shows that increasing
of the specific heat (i.e., the exponent α in a power law temperature does not increase the disorder in the spin
fit of the low-temperature specific heat, C T−α) that correlations in the z-direction.
v
is different from the model with continuou∼s disorder.4 We should add that the classical XY model has been
At even lower temperatures, a crossover in the scaling studied in the presence of Gaussian-random dilution53
behavior of Cv has been reported.7 Unfortunately, our and bimodal dilution54 of ferromagnetic bonds. In these
results do not converge at low enough temperatures to works, the BKT transition was seen to slowly disappear
observe such power laws. However, for T > 0.3, our re- as the dilution was increased. Here we only have consid-
sults are consistent with those in other studies4,7 (since ered the fully disordered case, i.e., an equal distribution
we consider Sz = 1/2,our temperatures are lowerby a of ferro- and antiferromagnetic bonds, so we do not ex-
factor of four from±those studies, which took Sz = 1). pect that any remnants of the BKT phase are presentin
±
our calculations.
B. XY model
C. Heisenberg model
Figures5and6showresultsforthespin-1/2XYmodel
on the square and honeycomb lattice, respectively. The Figures 7 and 8 show results for the spin-1/2 Heisen-
resultsforallquantitiesarewellconvergeddowntoabout berg model on the square and honeycomb lattices, re-
T 0.1 to 0.2. Figures 5(a),(b) and 6(a),(b) show that spectively. In Figs. 7(a) and 8(a), one can see that the
≈
the behavior of energy and the entropy is qualitatively plateau in the energy at low temperatures is the clear-
similar to the one observedin the Ising model. However, est of all models studied in this work. The onset of this
the results for the XY model in the presence of disorder plateau occurs at T 0.2 for both models. Remarkably,
≈
converge at lower temperature than those for the Ising in the honeycomb geometry, the energy in the presence
model. For the XY model, the specific heat in the pres- of disorder converges down to T 0.02. The entropy
≈
ence of disorder exhibits a peak that is well resolved by [Figs. 7(b) and 8(b)] behaves similarly to the entropy in
ourNLCE[Figs.5(c)and6(c)]. Wenotethattheenergy, the XY model, but also converges to very low temper-
entropy, and specific heat follow the second order result atures (T 0.03 to 0.04) in the honeycomb geometry.
≈
in Eq. (11) for T > 2 in the square lattice and T > 0.7 Like for the XY model, the specific heat exhibits a clear
in the honeycomb lattice. peak in the presence of disorder. The temperature at
Interestingly, Figs.5(d) and 6(d) show that in the XY which the maximum of that peak occurs is very close
model the uniform (z-)susceptibility in the presence of (but slightly larger) to the one in the clean model.
disorder exhibits a plateau for low temperatures. This In contrast to the XY model, the uniform suscepti-
is qualitatively different from the behavior observed for bility in the presence of disorder increases rapidly with
the Ising model. The fact that the response to an ex- decreasing temperature at the lowest temperatures ac-
ternal magnetic field in the z-directionis independent of cessible to us. The susceptibility therefore behaves qual-
0 0
(a) E (b) S (a) E (b) S
0.6 0.6
-0.2 Disordered Clean
-0.2 l = 13 l = 15
0.4 l = 14 l = 16 0.4
-0.4 DisollWr d==ey rn11en23d ClellWa ==ny n11n45 0.2 -0.4 WEuylenrn161 WEuylenrn172 0.2
-0.6 6 7
Euler Euler
11 12
0 0.4 0
(c) Cν (d) χ (c) Cν (d) χ 10
0.4 χ = 0.25/T
χ = 0.25/T 1 0.3
1
0.2
0.2
0.1 0.1 0.1
0 0
0.1 1 10 0.1 1 10 0.01 0.1 1 0.01 0.1 1 10
T T T T
FIG. 7. (Color online) Spin-1/2 Heisenberg model on FIG. 8. (Color online) Spin-1/2 Heisenberg model on
the square lattice. Same as Fig. 5 but for the spin-1/2 the honeycomb lattice. Same as Fig. 6 but for the spin-
Heisenberg model. 1/2 Heisenberg model.
9
itativelysimilartotheIsingmodel. Itisworthemphasiz- in spin-1/2 Ising, XY, and Heisenberg models. For all
ing that our results for all observablesin the honeycomb models we find that disorder leads to a saturationof the
lattice appear to converge at temperatures significantly energy at the lowest temperatures accessible to us, in a
belowT =0.1,andquiteclosetoT =0.01fortheenergy, regimewheretheentropyishigherthanincleansystems.
entropy, and uniform susceptibility. This makes apparent that there are many low lying en-
Our calculations for the Heisenberg model reachlower ergy states. For the disordered classical Ising model, on
temperatures and access regimes beyond what has been both the square and the honeycomb lattice, the diver-
possiblewithquantumMonteCarlosimulations. Forex- gence of the specific heat in the cleancase is replacedby
ample, Ref. 55 studied the case of diluting an antifer- a peak, and the uniform susceptibility follows an inverse
romagnetic model with ferromagnetic bonds of varying temperature law for all temperatures in the presence of
concentration. The results presented there did not reach disorder. This was explicitly verified order by order. In
theequalprobabilityJ = 1casediscussedherebecause the Heisenberg model, we also find that the susceptibil-
±
of the sign problem. For lower concentrations of ferro- ityincreaseswithdecreasingtemperature foralltemper-
magnetic bonds (below 30%), the calculations were still aturesaccessibleto usinthe presenceofdisorder. Inthe
limited to temperatures above T 0.3. XY model, on the other hand, we find that the suscep-
≈
tibility exhibits a plateau at low temperatures. On both
the XY and Heisenberg models, our NLCE calculations
V. SUMMARY were able to resolve a peak in the specific heat.
We have used numerical linked cluster expansions to
study thermodynamic properties of spin models with bi- VI. ACKNOWLEDGMENTS
modal ( J) bond disorder. The results reported are in
±
thethermodynamiclimitattemperaturesforwhichthey We acknowledge support from the National Science
are well converged. Foundation (NSF), Grant No. OCI-0904597, and the
Wehaveunveiledvariousinterestingeffectsofdisorder U.S. Office of Naval Research.
1 P. W. Anderson,Phys. Rev.109, 1492 (1958). 21 J. A. Blackman, J. R. Gonc¸alves, and J. Poulter,
2 K. Binder and A. P. Young, Phys. Rev.E 58, 1502 (1998).
Rev.Mod. Phys.58, 801 (1986). 22 A. K. Hartmann and A. P. Young,
3 M. Weissman, Rev.Mod. Phys.65, 829 (1993). Phys. Rev.B 64, 180404 (2001).
4 H. G. Katzgraber, L. W. Lee, and I. A. Campbell, 23 I. Gruzberg, N. Read, and A. Ludwig,
Phys. Rev.B 75, 014412 (2007). Phys. Rev.B 63, 104422 (2001).
5 F. Parisen Toldin, A. Pelissetto, and E. Vicari, 24 C. Amoruso and A. K. Hartmann,
Phys. Rev.E 84, 051116 (2011). Phys. Rev.B 70, 134425 (2004).
6 C. Thomas and A. Middleton, 25 J. Lukic, A. Galluccio, E. Marinari, O. C. Martin, and
Phys. Rev.E 87, 043303 (2013). G. Rinaldi, Phys. Rev.Lett. 92, 117202 (2004).
7 C. Thomas, D. Huse, and A. Middleton, 26 H. G. Katzgraber and L. W. Lee,
Phys. Rev.Lett. 107, 047203 (2011). Phys. Rev.B 71, 134404 (2005).
8 P. Henelius and A. Sandvik, 27 W. Atisattapong and J. Poulter,
Phys. Rev.B 62, 1102 (2000). New Journal of Physics 11, 063039 (2009).
9 M. Troyer and U.-J. Wiese, 28 N. D. Mermin and H. Wagner,
Phys. Rev.Lett. 94, 170201 (2005). Phys. Rev.Lett. 17, 1133 (1966).
10 B. Tang, D. Iyer, and M. Rigol, arXiv:1411.0699. 29 A. Altland and B. Simons,
11 M. Rigol, Phys. Rev.Lett. 112, 170601 (2014). Condensed Matter Field Theory (Cambridge Univ. Press,
12 B. Wouters, J. De Nardis, M. Brockmann, 2006).
D. Fioretto, M. Rigol, and J.-S. Caux, 30 V. L.Berezinskiˇi, Sov.JETP 34, 610 (1972).
Phys. Rev.Lett. 113, 117202 (2014). 31 J. M. Kosterlitz and D. J. Thouless,
13 M. Rigol, Phys. Rev.E 90, 031301(R) (2014). J. Phys. C: Solid State Phys.6, 1181 (1973).
14 M. Sykesand M. E. Fisher, Physica 28, 919 (1962). 32 K. Harada and N. Kawashima,
15 Z. Weihong, J. Oitmaa, and C. Hamer, J. Phys. Soc. Jpn. 67, 2768 (1998).
Phys. Rev.B 44, 11869 (1991). 33 J. Carrasquilla and M. Rigol,
16 A. Di Ciolo, J. Carrasquilla, F. Becca, M. Rigol, and Phys. Rev.A 86, 043629 (2012).
V. Galitski, Phys. Rev.B 89, 094413 (2014). 34 V. L.Berezinskiˇi, Sov.JETP 32, 493 (1971).
17 L. Onsager, Phys.Rev.65, 117 (1944). 35 J. Frhlich and T. Spencer,
18 S.F.EdwardsandP.W.Anderson,J.Phys.F:Met.Phys. Comm. Math. Phys. 81, 527 (1981).
5, 965 (1975). 36 R.Kenna, arXiv:cond-mat/0512356 (2005).
19 M. E. Fisher, Physica 25, 521 (1959). 37 D. Gross and F. Wilczek,
20 S.ChoandM.P.A.Fisher,Phys. Rev.B 55, 1025 (1997). Phys. Rev.Lett. 30, 1343 (1973).
10
38 D.J.Gross,Nucl.Phys. B - Proc. Suppl.74, 426 (1999), 46 E. Khatami and M. Rigol,
{QCD} 98. Phys. Rev.B 83, 134431 (2011).
39 A. Cuccoli, V. Tognetti, R. Vaia, and P. Verrucchi, 47 E. Khatami and M. Rigol,
Phys. Rev.B 56, 14456 (1997). Phys. Rev.A 84, 053611 (2011).
40 C. Domb, Advancesin Physics 9, 149 (1960), 48 B. Tang, T. Paiva, E. Khatami, and M. Rigol,
http://dx.doi.org/10.1080/00018736000101189. Phys. Rev.Lett. 109, 205301 (2012).
41 C. Domb, Advancesin Physics 9, 245 (1960), 49 B. Tang, T. Paiva, E. Khatami, and M. Rigol,
http://dx.doi.org/10.1080/00018736000101199. Phys. Rev.B 88, 125127 (2013).
42 M. F. Sykes, J. W. Essam, B. R. Heap, and B. J. Hiley, 50 B.Tang,E.Khatami, andM.Rigol,Comput.Phys.Com-
Journal of Mathematical Physics 7, 1557 (1966). mun.184, 557 (2013).
43 M. Rigol, T. Bryant, and R. R. P. Singh, 51 E. Khatami, R. R. P. Singh, and M. Rigol,
Phys. Rev.Lett. 97, 187202 (2006). Phys. Rev.B 84, 224411 (2011).
44 M. Rigol, T. Bryant, and R. R. P. Singh, 52 K.Huang,Statistical mechanics (Wiley,NewYork,1987).
Phys. Rev.E 75, 061118 (2007). 53 Y.-B. Deng and Q. Gu,
45 M. Rigol, T. Bryant, and R. R. P. Singh, Chin. Phys. Lett.31, 020504 (2014).
Phys. Rev.E 75, 061119 (2007). 54 T. Surungan and Y. Okabe,
Phys. Rev.B 71, 184438 (2005).
55 A. W.Sandvik,Phys. Rev.B 50, 15803 (1994).