Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Disorder-Free Glass Transitions of Spins and Orbitals
in a Frustrated Pyrochlore Magnet
1

Kota Mitsumoto1 , Chisa Hotta2 and Hajime Yoshino1,3
Graduate School of Science, Osaka University, Toyonaka, Osaka 560-0043, Japan
Osaka University, Toyonaka 560-0043, Japan.
2
Department of Basic Science, University of Tokyo, Tokyo 153-8902, Japan
3
Cybermedia Center, Osaka University, Toyonaka, Osaka 560-0043, Japan

Abstract:

glasses obtained by supercooling molecularliquids or making densely packed soft-matterslike colloids. The most important basic question there is to clarify whether the so-called
Kauzmann transition[9], which is a putative,
idealized thermodynamic glass transition, cantake place in the supercooled liquid state.
Recently important progress was made on thetheoretical side: a mean-field theory which
combines the density functional theory of
liquids [10] and the replica method used instatistical mechanics of disordered systems[11] established that the ideal glass transition
exists (at least) in the large dimensional limit(see for a review [12]). However, the fateof ideal thermodynamic transition in finite
dimensions is largely unknown [13].
The
problem is difficult to solve precisely becauseof the glassiness: the extremely sharp increaseof the relaxation time (viscosity) going deeper
into the supercooled liquid state makes it very
difficult to measure quantities in equilibrium.
Another reason is that the super-cooled liquidstate is a meta-stable state which should
be replaced by the crystalline state in true
equilibrium. A possible way out of these
difficulties is to seek simple lattice models,
which mimic the basic phenomenology ofstructural glasses, carefully designed to avoid
crystalline transitions, and allow various
simulation techniques developed for latticemodels [14].

In this activity report, we present an
overview of our recent extensive Monte Carlosimulations performed on the ISSP supercomputer (SYSTEM B) [1][2][3] to study our effective theoretical model for the disorder-freespin-glass transition observed experimentallyin the pyrochlore oxide Y2 Mo2 O7 . By taking
into account not only the spin but also the orbital degrees of freedom of the Mo ions, wefound a new type of glass transition - spinorbital glass transition.

1

Introduction

Glass is a generic state of matters which canbe found in a diverse range of systems ranging from soft to hard condensed matters. Yet,
unlike crystalline states with long-ranged ordering of periodic structures, the very mechanism(s) of the emergence of glasses remainquite elusive. The exceptional case is the family of disordered magnets with quenched disorder, i. .e. spin-glasses. In spin-glasses, the existence of thermodynamic glass transitions isestablished by experiments, theories and simulations (see for reviews [4, 5, 6, 7]). But apparently, glasses without quenched disorder ismore ubiquitous in nature and in industrial
materials. Understanding the mechanism of
the emergence of glasses without quenched disorder remains as one of the most important
unsolved problems in physics [8].
Majority of glassy systems without
quenched disorder emerge from supercooled
liquids. Typical examples are the structural

In principle, glass transitions similar to
those found in the supercooled liquids canalso take place in spin systems. Indeed one
can show exactly that the ideal glass tran-

35

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

imental system, this observation means thatthe model is too crude to capture the real system which exhibits a spin-glass transition. In
order to explain the spin-glass transition of theY2 Mo2 O7 system, previous theoretical studies assumed the presence of some amount of
quenched disorder [28, 29]. However, from themean-field theoretical point of view, disorderfree spin-glass transitions are possible [15].
Albeit in a somewhat artificial case with ahighly non-linear interaction , where the model
is built to describe an optimization problem
like the graph coloring with continuous colors, a 2nd order disorder-free spin-glass transition involving full replica symmetry breaking
emerges without passing through the superparamagnetic phase has been found (See sec.
10 of [15]).
The purpose of the present work is to
uncover the microscopic mechanism of thedisorder-free spin-glass transition observed in
the pyrochlore magnet. Our key observationin the present work [1], which is motivatedby a recent experiment [23], is that in the pyrochlore oxides, not only the spins but also theorbitals (lattice displacements) exhibit glasstransitions simultaneously - spin-orbital glasstransition. In the following, we first explain
our theoretical model and then present our results obtained by extensive Monte Carlo simulations at ISSP.

sitions can take place in spin systems without quenched disorder, by examing supercooledparamagnetic states in the large dimensional
limit [15]. However, in reality, the metastability limit of the paramagnetic state toward
long-ranged ordering of periodic structures(corresponding to the Kirkwood instability ofliquids [16]) may preempt the ideal glass transition (see for instance [17]).
A very promising playground to explore
glass physics is geometrically frustrated magnets on corner-shared lattices such as the
Kagome lattice and pyrochlore lattice (seeFig. 1 (a)). In general, glassy states appear
when periodic long-range ordering is avoided
down to very low temperatures or high densities where strong interactions become important. Therefore frustration is a key ingredient for glass physics [18]. The most spectacular example is the pyrhoclore oxcide Y2 Mo2 O7 ,which is known experimentally to exibit clearcut spinglass transitions without appreciable
quenched disorder [19, 20, 21, 22, 23, 24, 25].
Moreover, in sharp contrast to the structural glass transition mentioned above, thespin-glass transition emerges directly out ofthe paramagnetic state as a 2nd order transition much like the conventional spin-glasstransitions[4, 5, 6, 7], without the need to go
through the metastable ’supercooled’ paramagnetic state. Actually, quite remarkably, thephenomenology of the spin-glass transition andeven the critical exponents associated with thespin-glass transition [19, 20, 25] are nearly indistinguishable from that of the conventional,
canonical spinglasses with quenched disorder.
On the theoretical side, microscopic explanation of the disorder-free spin-glass transition remained a big challenge for more than 30
years. The principal magnetic interaction in
the pyrochlore oxide Y2 Mo2 O7 is the antiferromagnetic interaction between the magnetic
moments of the Mo ions sitting on verticies ofthe pyrochlore lattice. It has been establishedthat the purely antiferromagnetic Heisenberg
spin model on the pyrochlore lattice exhibits
no transitions down to T = 0 because of thevery strong geometrical frustration [26, 27].
This is remarkable given that the system isthree dimensional. Coming back to the exper-

2

Model

We consider two kinds of dynamical variables
associated with the Mo4+ (4d2 , S = 1) ionswhich sit on the vertices i = 1, 2, . . . , N of
the pyrochlore lattice. The first is the classical
Heisenberg spin Si = (Six , Siy , Siz ) normalized
as |S| = 1 which represents the magnetic moments of the Mo ion.
The second is the ’displacement’ of the Mo
ions σi = (σi,x , σi,y , σi,z ). The recent experiment [23] suggests that the pyrhocolore lattice is distorted by displacements of the Mo
ions. The displacement of a vertex (the Moions) take place along the lines connecting thetwo tetrahedra that share the vertex. Thus
the displacements of the vertices of a given
tetrahedron are either pointing toward or away36

CHAPTER 1. INTRODUCTION

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

(b)

?

a)

d)

b)

?
?

c)

strated lattices
with antiferromagnetic interactions on a triangle (top) and a tetrahedron (bottom). (b):simplexes.

%(')%* $)/$,$)$&(')"* $).$,$)$&(')#* $)0$,$)$&(')$* $)++,$)$&-"#% lattice.
Figure 1: Spins on deforming pyrhoclore(a): Pyrhoclore lattice. (b): O2 ions and(')%* $)/$,$)$&$
4+"
the numbers
1 − 4 are the sublattice indices ofmagnetic Mo ions around the i site, where(')"* $).$,$)$&(')#* $)0$,$)$&etrical Frustration!#(')$*$)++,$)$&represent
thepositionsof the Jahn-Teller distorted i
Mo4+ ions. Red and blue dashed circles
#
4+
2
!"
&#$
&'$
&"$
ions (i, j) interact
throughthe O ion as shown in the inset panel,ion. Spinson representative
a pair of Momagnetic
we review some basic aspects of geometrical frustration
in some$st of all let us note that their thermodynamic properties
unusual.
In Mo-O-Mo
ordinary systems,
wherearealpha
is theangle and r̂ij is the unit vector in the i → j direction. (c): Ice-type!#
ase monotonically by lowering the temperature and displacements
the system finally of
reachgroundstate
thetheMotetrahedron.
The different color bonds represent different exchange
is the lowest and the entropy becomes 0, satisfying the third law of thermodynamics. In this
!"
&#$
&'$
&"$
interactions. (taken from [1, 2])transition may take place from a disordered phase to an ordered phase, e.g. ferromagnets
the other hand, strongly frustrated systems may not exhibit any long-ranged order evenower than the energy scale of the interaction and opens various possibilities of unexpectedfrom the center of the tetrahedron. Furthertemperatures. In particular, some geometrically frustrated spin system shows the disorderthe yet.
experimentnsition whose microscopic mechanism has not beenmore,
understood
The goal ofsuggests
this thesisthat the disoblem for the case of a pyrochlore magnet Y2 Mo2 O
placements
follow the ’ice-rule’ [30] : the ’2-in7.
rustration in spin systems means a situation where 2-out’
spins instructure
the system (seecannotfind1a (d)).
config-As suggestedFig.
ly satisfies all the interactions. The simplest example
is
the
case
where
three
Ising
spins
on
by the experiment, we parametrized
the disgle triangle interact with each other antiferromagnetically (See Fig. 1.7 (a)). If we put an up
placements
as
σ
=
σ
ê
.
Here
we
introducediiantiferromagν
and a down spin on the site 2, the remaining spin on the site 3 cannot satisfy
ancan
Ising
variable
σi =to ±1
to represent ’in’
-3 and 2-3 simultaneously. This frustrated situation
be easily
generalized
simplexess, tetrahedra, e.t.c with spins put on the vertexes andand antiferromagnetic
interactions
along on êν which’out’ displacements
dependingwith
Jσi ,σj = J[1 + δ(r̂ij · σi + (−r̂ij ) · σj )] (δ > 0),(2)
The 1st term of Eq. (1) describes the exchange interaction between the spins on thenearest-neighbour sites < i, j > on the pyrhoclore lattice. The exchange interaction Eq. (2)between the spins on the i-th and j-th Mo ions,with the energey scale J(> 0), is mediated by
the O ion and depends on the angle α of the
Mo-O-Mo bond. (see Fig. 1 (b)(c)) Thus theexchange interaction depends on the displacements σi and σj of the Mo ions. This form
is constructed in a way to reproduce the values of the effective interactions derived microscopically from the perturbation process on a
Kanamori-type of Hamiltonian [1, 2]. Moreprecisely we parametrize it as Jσi ,σj = (1 +2δ̃)J (in, in),
√ J (in, out), (1 − 2δ̃)J (out, out),
where δ̃ = 6δ/3.
The 2nd term of Eq. (1) represents the elastic energy of the Mo4+ displacements. The
elastic energy is minimized if the ice-rule is satisfied.
There are three parameters in this system;T̃ = kB T /J is the dimensionless temperature,˜ = /3J the ratio of the energy scales between
the exchange interaction and the elastic energy

is a unit vector in the [111], [11̄1̄], [1̄11̄] andonsider extended networks made of such simplexes,
suchdirections
as the triangularlattice, kagomé[1̄1̄1]
respectively
for the sub-lattices
hlore lattice with antiferromagnetic coupling on the edges (See Fig. 1.7 (b)). The important
ν
=
1,
2,
3,
4
which
the
i-th
spin
belongs to (Sees that they are corner sharing: adjacent simplexes are connected to each other only throughFig. 1.(b)). The microscopic mechanism of thepeculiar displacements is presumably a JahnTeller effect: the displacements break the 3fold rotational symmetry of the trigonal crystal
field for the Mo ions resulting in the splittingof its eg orbital [3]. Thus we regard σi as thevariable which represents the configuration of
the orbital of the Mo ion.
Our effective model for the pyrhoclore magnet Y2 Mo2 O7 is given by the following effective
Hamiltonian,

H=



<ij>

Jσi ,σj Si · Sj − 



<ij>

σi · σj ( > 0)(1)

37

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo
2.5. COMPLEX ENERGY LANDSCAPE OF OUR EFFECTIVE MODEL

65 ENERGY LANDSCAPE OF OUR EFFECTIVE MODEL
2.5. COMPLEX

2.5. COMPLEX ENERGY LANDSCAPE OF OUR EFFECTIVE MODEL

65

65

√
of decoupled
the displacement,
6δ/3ouris theam- is decoupled into (1)= 0. Then our system is
into (1)and= δ̃0.=
Then
system(a)is decoupled
(a)
=
0.
Then
our
systeminto
plitude of the displacement
call(1)(a)the spin part with purely antiferromagnetic
in-the spin(hereafterpart with wepurelyantiferromagnetic
ina)
the
spin
part
with
purelyantiferromagnetic
insimplypart
as T,with, δ). At  →[26,∞, 27]
the andlattice the orbital part withteractions [26, 27] and (2) themthe orbital
teractions [26, 27]teractionsand (2) the orbital
part (2)withdistortion
becomesstatic.
pure ’spin-ice’ type interaction
[31,
32].
Both
pure
’spin-ice’
typeinteraction[31, 32]. Bothpure ’spin-ice’ type interaction [31, 32]. Bothof them have essentially flatFinallyenergylandscape
let
us
note
what
happensif
weswitchofthem
have
essentiallyflatenergylandscapeof them have essentially flat energy landscapethecoupling
between
the
spins
and
orbitalsso that the two degrees ofoff
freedomremaindisso thatofthe
two degrees
of freedom remain disso that the two degrees
freedomremain disdisplacements)
δ=
Then ourordered down to T = 0. (latticeordereddown to T=
0. 0. downorderedto Tsystem= 0.
is decoupled into (1) the spin part with purely
antiferromagnetic interactions [26, 27] and (2)b)c)
3 Simulation
3 Simulation Methods
Simulation
Methods
the orbital(lattice) 3
partMethods
withpure ’spin-ice’(c)(b)(c) (b)(b)(c)
typeIninteraction[30, 31].simulations
Both of them
have esorder
toperforminequilbrium
In order to perform simulations
inflatequilbrium
In order tosoperform
in equilbrium
sentiallyenergythat the simulations
twoourwe
combinedthelandscapefollowing methodsin
we combined the following
methods
inourwecombinedthefollowing
methods
in ourdegrees of freedom
remain disordered down to
Carlo simulations.
Monte Carlo simulations. T = Monte
Monte Carlo simulations.
0.Single spin update:
The usual single spin
Single spin update: The usual single spin Single spin update: The usual single spin
updates by the conventional Metropolis algoupdates by the conventional Metropolis algo-updates by the conventional Metropolis algorithm is employed to update the Ising variables
rithm is employed to update the Ising
variablesrithm is employed to update
theFigure2.14: Ising
Complex variables
energy landscape of our e ective modeli = ±1 (i = 1, 2, . . . , N ).
Schematic picture of the energy landscape of our e↵ective model (right) obtained by combining a flat spin
Figure2.14: Complex energy landscape of our(a):e ective
model
3
Simulation
Methods
Figure2.14: distortion
Complexenergy(leftlandscapee ective
=
±1
(i
=
1,
2,
.
.
.
,
N
).
energylandscape
(left
top)
and a zigzag but
degeneratelandscapebottom). of(b):our
Systemsize modeli= ±1picture(iWe=
1,energy
2, adopted.landscape. . , Nof).ourae↵ective model (right) obtained
Loop update (a):method:alsoi Schematicof theby combining
a flat spin
Schematicof theenergydi↵erent
landscapeof oure↵ective
the fluctuation
ofpicture
the energy
betweenenergyminimawith model (right) obtained by combining a flat spin
L dependence of V (E) which quantifies (a):energylandscape
(left
top)
and
a
zigzagbut
degeneratedistortion
landscape(leftbottom).
(b):
Systemsize
energylandscape
(left
top)
and
a
zigzagbut
degeneratedistortion
landscape (left bottom). (b): System size
Loop update method: We
also adopted
aL dependence(E) which quantifies fluctuation
a commonlattice-ice configuration.
(c): Systemsize L dependence of V
Loop
method:We
adoptedaLconfigurations.
nonlocalupdate method
calledupdate(E)thewhichloop
quantifies
the fluctuation
the also
energybetweendi↵erent
energyminima withof Vupdateof theof
energybetweendi↵erencelattice-ice
which quantifies
the fluctuation of the energy between di↵erent energy minima withdependenceof V (E)(a):Figure2:
Energylandscape.
Schematic
Figure2:
Energylandscape.
(a):
Schematic(E)
whichquantifies
fluctuationacommonlattice-ice
configuration.
(c):
Systemsize
LdependenceofV(E) which quantifies fluctuation
a common lattice-ice configuration. (c): System size L dependence of V
nonlocal update method called
the loop
updateIn orderto perform
simulationsinequilibalgorithm,
used tosimulate
spin-ice
systems
nonlocalupdatemethodcalledpicture
the loopof the energy
betweendi↵erencelattice-ice configurations.
Figure2:
Energylandscape.
(a):
Schematicof updateenergyTheessentiallyoflandscape.
the energy
between2:
di↵erencelattice-ice
configurations.
FigureEnergylandscape. (a): Schematic
pictureofenergylandscape.
Theessenalgorithm, used to simulate
spin-ice
systems. simulate
updatetheIsingvariables
rium,[33],wetocombined
thefollowingin
itoalgorithm,
usedspin-ice
systems
flat
energylandscapeof the (withspins(left/top)
pictureof methods
energylandscape.
Theenergy
essentially
pictureofenergylandscape.
The essentiallytiallyflat
landscapesmany
enMonteCarloreflection
method:Toupour MonteCarlosimulations.
.
[33], to update the Ising variablesi.
[33],flat
to updatetheIsingvariablesiandorbitals(left/bottom)
intheabsenceof of the spins (left/top)energy landscape ofergeticallythe spinsdegenerate(left/top)flat
energylandscapeminima)
of
the
spins
(i =single
1, 2,. .spin. , N ), method:date theHeisenberg
spins
Susuali CarloMonte Carlo reflectionSingle
method:To
up- Monte
reflection
To upspin
update:The orbitals
theincoupling= 0 changes
into complex
oneand(left/bottom)
theandabsenceorbitals(left/bottom)
in the absence of(left/top)orbitalsofand(left/bottom)
in the abweused
thesinglespin
updatesbythe
1, 2,the. . . conventional, N ),datethe
date the Heisenberg spinsupdates
Si (i = by
Metropolisinthepresenceofthecopling(right).
(b),(c)
(i
=
1,
2,
.
.
.
,
N
),
thecoupling
Heisenberg=algospinsSi
0 changes
complexonechanges into complex one
senceintoof thecoupling
δthe= 0couplingchange into=a 0comMetropolisbyreflection
method [34].
we used the single spin
updatesvarianceofby
the theenergybetween di↵erent energyrithm
is employedthe
forwetheIsingvariables
σi the=used
thesinglespincopling
updatesinthepresenceof(right).
(b),(c)inthepresenceofthe
one in the presence of the coupling δ > 0copling (right). (b),(c)Over relaxation method: We also used plex
minimaobtained by quench experiments (seeMetropolis reflection method±1 (i[34].
= 1, 2, . . . , N ). Metropolis
reflection
method[34].
energydi↵erent
energyvarianceof the energy
the over relaxation variancemethod ofto theupdatethebetween(right).
(b),(c)varianceof the energy
betweenbetween di↵erent energy
textfor
details).
Over relaxation method:We
also usedOver
method:Wealso usedLoop
updatemethod:Werelaxation
also
adopteda quenchminimaobtainedbyexperiments(seeminima
Heisenberg
spins.
[35]
differentenergyminimaobtainedby quenchobtainedby quench experiments (seethe over relaxation method
toupdateexchange
thethe text
over
relaxation
methodto
updatethenonlocalupdatemethodcalledthe
loop
updateReplicamethod:
On
top
of
the
for details).
experiments (taken fromtext[2])for details).
Heisenberg spins. [35]Heisenberg
spins.systems[35]
algorithm,
used to we
simulate
spin-ice
above methods
used thereplicaexchange
Replica exchange method:
top
oftothe[32], method
to On
for thevariables
σi .exchange
Replicamethod:On top of the[36]Isingacceleratetheequilbrationofabove methods we used theMontereplicaexchange
4
Results
above
methods
weusedthereplicaexchange
the system.
Carlo reflection method: For the
method [36] to accelerateHeisenberg
the equilbrationofmethodto),accelerate= 1, 2, [36]...,N
we used the equilbration of
spins Si (i1st term of the Hamiltonian Eq. (1) by zero
the system.
thesystem.
the single
spin
updatesbytheMetropolisre- 4.1 A simple analysis of the energy
4 Results
temperature, greedy Monte Carlo updates oflandscapeflection method [33].
the spins.
procedure
ob1st
term
of
the
Hamiltonian
Eq. (1)Repeatingby zero1stthis
termof the we
Hamiltonian
Eq. (1) by zero
4.1
A
simpleanalysisoftheenergyOver-relaxation
method:
We
also
used
taineda largenumberofenergyminimaofvar4 Results
4
ResultsAs
wementioned
at
the
end
of
sec.
2,
the
entemperature,
greedyMonteCarloupdatesoftemperature, greedy Monte Carlo updates oflandscapemethod to update thethe over-relaxation
ious
energies E.
Fig.
2 (b)becomes
we display
theergy
landscapeof In
oursystemessenthe spins. Repeating thisprocedure
weobthe
spins. Repeating
2i
2 this procedure we obHeisenberg
spins.
[34]
varianceoftheenergy
V
(E)
=
h
E
h
Ei
4.1 A simple analysisoftheenergyspin
Aof theenergyAs we mentioned4.1
in theendsimple
sec. number
2,analysis
the en-oftiallyflat(with ofmany
energeticallydegenertained
aoflarge
energyminimavar-tained
a large number of energy minima of varbetweendi↵erentof thespins withlandscapeReplicaexchangemethod:On
top
theFig.ate
ergy landscapeof oursystembecomesessenlandscapeminima)if we realizations
switchthecoupling
ious
energies
E.of In2 (b)we display
theiousoff
energies
E. In beFig. 2 (b) we display theacommonspin-icelike
orbitalconfiguration
tiallyflat if weo↵coupling between tween the spins
above
methods
weswitch
usedvariancethethereplica-exchange
2 i and2
orbitals(lattice
displace2
oftheenergy
V
(E)
=
h
E
h
Ei
varianceof theh Ei 2spin
As we mentioned in the end of
2, theorbital
en-As degrees
wethe
mentionedin theofsec. 2, overthe en(averageddi↵erent
realizationsof energy
such or-Vspin (E) = h E ithesec.
spinof freedom=
0. ments)
method[35] andto accelerate
equilibrationof end
,
i.
e.
δ
=
0.
Here
let
us
show
the
rebetween di↵erentof thespins
betweenthe
di↵erentof the spins withergy landscape of our system
becomesessen-aergy
bitals).
Hereh . . .iwith
represent
averagerealizations
overourrealizations
systembecomesessenHereletus show
resultlandscapeof a simpleofquench
ex- sult
the wholesystem.
of a simple
quencha experiment
which
proa
commonspin-icelike
orbitalconfiguration
spin-ice
like orbital configuration
tially flat if we switch o↵ the couplingdi↵erent realizations
the spins witha comtially
flatinsight
if we switch
the coupling
between of common
perimentsbetween
obtainsome
on howo↵the
vides
some
insights
on(averaged
howandthe .overenergyland- realizations of such orWe
considertothe
periodic
systems
ofdi↵erent
cubic realizations(averaged
over
of
such
ordi↵erentmonorbitalconfiguration.
.
represents
the spin and orbital degrees ofenergy
freedom= 0.thespin
orbital degrees
landscapechanges
inand
the presenceof the of freedom = 0.
3 unitscapechangesin thepresenceoftheh .coupling
cells
with
togeometry consisting of L
bitals).
Here
h
.
.
.i
represent
averagethe theaverageovber
di↵erent
thebitals).realizations
Here. .iof represent
the average over
Here let us show a result of a simple
quenchexcoupling.
First
wegenerated’spin-ice’configHere
let
us
show
a
resultof a simple
quenchex-over3> spins
0. First
wegenerated’spin-ice’configtally N = 16L spins, andperform
120 sta- of δtheorbitals.
In
(c)
we
show
the
correspondingdi↵erent
realizations
withacomdi↵erent
realizations
of
the
spins
with a comurationsof’srandomly.
Note
that
any
of
periments to obtain some insight on how i theperiments to obtain some insight on how theofσrepresents’s randomly.
Note between
that anydifof
tistically independent runs
fororbital
the systemsize urations
ionetomeasurethefluctuationmonconfiguration
and
.
.
.
them
minimizethe
2nd termof thechanges
Hamiltomon orbital configuration and . . . represents
energy landscape changes in thepresenceof theenergylandscapein themthe presenceof theminimizesthe
2ndterm of the
HamiltoL =nian
4, 5,
6, 8,
evaluating
theaveragesandferentrelizations
ofthe
the
=
Eq.
(1).
Second,we
performedasimpletheaverageovber
di↵erent
realizationsofthe orbitals
averageVlatticeovber(E)di↵erent
realizations of thecoupling. First we generated’spin-ice’ configcoupling.
FirstInwethe
generated’spin-ice’config2 Second, we performed a simple
nian
Eq.
(1).
mean-squared
errorsofobservable.
fol2
energydescent
simulation
withrespect
to
the
i
h
Ei
h
E
orbitals.
In
(c)
we
show
the
corresponding
In (c) we show the corresponding
urations of i ’s randomly.lowing
Noteanalysis,that anyurationsfocus
of ion’s randomly.
Note thatanysimulationof orbitals.
descent
with respect to theweofmainly
a repre- energyonetomeasurethefluctuation
betweendifonetomeasurethe fluctuation between difthem minimize the 2nd termof theHamiltothe of
Hamilto1stofterm
the Hamiltonian Eq. (1) by zerosentative
systemat δthem= 1.5,minimize
 = 0.6. the 2nd term
4
relizations
ofwe
theperformed
orbitals Valattice(E) =ferent relizations of the orbitals Vlattice (E) =nian Eq. (1). Second, we performed a simplenianferent
Eq. (1).
Second,simple222energy descent simulation with respect to theenergyh E i descenth Ei simulation with respect to the h E 2 i
h Ei
spin

lattice

spin

spin

lattice

438

lattice4Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

78

CHAPTER 3. NUMERICAL SIMULATION

temperature, greedy Monte Carlo updates ofthe spins with orbitals σi being fixed. Repeating this procedure we obtained a large numberof energy minima of various energies E. In
Fig. 2 (b) we display the variance of the energy Vspin (E) = E 2  − E2 between different
realizations of the spins with a common orbital(lattice) configuration (averaged over different
realizations of such orbitals). Here . . . represents the average over different realizations of
the spins with a common orbital (lattice) configuration and . . . the average over different realizations of the orbitals. In Fig. 2 (c) we show
the corresponding one to measure the fluctuation between different realizations of the or2
bitals (lattice) Vlattice (E) = E 2  − E . Both(b) and (c) show that the variance grows withthe system size L suggesting lugged, complex
energy landscape.

4.2

(a)

C (t)

(c)

(d)

⌧

⌧s

Cs (t)

(b)

T

Tc

T

Tc

(f)

C̃ (t/⌧ )C̃s (t/⌧s )(e)

Critical slowing down without
long-ranged ordering
Figure 3.4: Relaxation of spins and distortions for L = 6, ✏ = 0.6(a): Spin auto-correlation function. (b): Orbital auto-correlation function. (c): Typical relaxation timeof spins. (d): Typical relaxation time of distortions. The black solid lines in the figure (b) and (d) deBy lowering the temperature we observed that
note the fits by ⌧ = A(TT)
with (A, T , z⌫) = (3.99(89), 0.0695(19), 2.79(12)) and (A, T , z⌫) =Figure3.55(4))
3: Autocorrelation
functionsand plots
relax(4.27(31),by using the least squares method.
(e), (f): Scalingof auto-correlation
the dynamics of both the spins and orbitals 0.0686(5),functions using their relaxation times. The solid lines represent the exponential fits.
ations
times
of
the
spins
and
orbitals.
L=
6,(lattice displacements) slow down significantly.
In Fig. 3 (a)-(b) we display the time autocorre-  = 0.6. (taken from [1])lation functions of the spins Cs (t) and orbitals
Cσ (t) plotted against time t. Here we used ing. Some pinch points can be seen in Sσ ,’single spin updates’ for the spins and orbitals for example at (1/2, 1/2, 1/2), at low temperto measure the dynamical observables, start- atures as well known in spin-ice systems [31].
ing from the initial configurations equilibrated
To summarize the above observations sugusing all protocols listed in sec 3. As shown in gests a simultaneous, 2nd order glass transi(c)-(d), the relaxation times grow in lowering tions of the spins and orbitals (lattice displacethe temperature following a power law. The ments) at a common temperature, into a new
independent power-law fits on the two observ- glass phase - spin-orbital glass phase, whichables suggest a common critical temperature emerges directly from the high temperature
Tc ∼ 0.07 with different exponents zν ∼ 4.2 paramagnetic (liquid) phase.
(orbital) and zν ∼ 3.5 (spin). As shown in (e)(f) the auto-correlation functions follow scaling
4.3 Non-linear susceptibilities
laws in terms of scaled times. We have checkedthat finite-size effects of the auto-correlation The signatures of the spin-orbital glass can befunction is negligible within the temperature detected by observing non-linear susceptibiliand time scales shown in the Figure. We have ties of the two degrees of freedom. As shown
found that the same type of scaling holds for in Fig. 5, the non-linear susceptibilities com = 0.65, with Tc ∼ 0.086 and almost the same puted using fluctuation formulae (see [1] forexponents as obtained above [2, 3].
the details) at lower temperatures grow negaAs shown in Fig. 4, the structure factor tively rapidly when increasing the system sizeof the spins Ss and orbitals (lattice displace- L. This observation provides an interestingments) Sσ show no hints of long-range order- suggestion for experiments: by observing not
c

39

z⌫

c

s

c

3.4. SIMULATION RESULTS

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

84

87

CHAPTER 3. NUMERICAL SIMULATION

(a)

(a) 1

(b)

0.016

( )
3

0.014
0.012

3

0.5

0.01
0

0.008
0.006

-0.5

0.004
0.002

-1
1

0
0.025

(b)

T

T
(c)

0.02

0.5

(d)

0.015
0

-1
-1

| 3|

0.005
-0.5

0

0.5

1

-1

-0.5

0

0.5

1

( )
3

0.01
-0.5

0

(c)

T

T

Figure 3.12: Magnetic and dielectric nonlinear susceptibilities at ✏ = 0.6(a), (b): Size dependence of magnetic and dielectric nonlinear susceptibilities. The dashed line represents
the freezing temperature Tc ⇡ 0.07 determined by the auto-correlation function. (c), (d): logarithmic scaleof the absolute values of the nonlinear susceptibilities.

Figure 5: Non-linear susceptibilities of thespins and orbitals (lattice displacements).  =0.6. Panels (a),(c) display the non-linear susceptibilities of the spin while (b),(d)display
the non-linear susceptibilities of the lattice displacements (orbitals). (taken from [1])(d)Figure 3.9: Structure factors for L = 6, ✏ = 0.6 above Tc (T = 0.4) in the left side and below Tc(T = 0.0625) in the right side(1)(a) and (b) show the spin structure factor Ss (0, h, l) and the orbital structure factor S (0, h, l), respectively.
(2)(c) and (d) show the spin structure factor Ss (h, h, l) and the orbital structure factor S (h, h, l), respectively.

ventional type, i. e. the one due to quencheddisorder, may take place at a lower temperature.
To get some insight into this issue within
the available computational resources, we analyzed the heat capacity C and the fraction Pof the tetrahedra which satisfies the ice ruleat various values of . At  = 0.6, 0.65, forwhich we have determined the critical temperature of the spin-orbital glass transition, we
found that the heat-capacity exhibit a peak atslightly higher temperatures, similarly to thecase of the conventional spin-glasses [4] andthe experimental result in the Y2 Mo2 O7 system [36]. We have checked that the finite-size
effects of the heat capacity are very weak [1]).
As shown in Fig. 6, the heat-capacity exhibit
an additional peak at a higher temperature atlarger values of  > 1.0. The peak at thehigher temperatures appears to follow the peakof dP/dT . This observation suggests that atlarger values of , the liquid state of the lattice
distortions exhibits a smooth crossover frompurely random one to more spin-ice like one
at higher temperatures. On the other hand,the other peak of the heat capacity located at

Figure 4: Heat maps of the structure factorsof the spins and lattice displacements. L = 6,
 = 0.6 above Tc (T = 0.4) in the left side andbelow Tc (T = 0.0625) in the right side. (takenfrom [1])only the usual magnetic non-linear susceptibility, which is used often to study conventional spin-glasses [4, 6, 7], but also dielectric non-linear susceptibility, the simultaneous
spin-orbital glass transitions suggested by thepresent work may be detected experimentally.

4.4

Static or dynamic distortions?

In our study we regard the lattice distortions (orbitals) as dynamical degrees of freedom rather than as frozen-in, static quencheddisorder. But by increasing the energy scale ofthe lattice displacements given by the parameter  in Eq. (1), one would wonder whether thelattice displacements can freeze independently
from the spins at higher temperatures. If thishappens, the spin-glass transitions of a con-

40

86

CHAPTER 3. NUMERICAL SIMULATION

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo(a) 1.7

1.5

C 1.4

1

(b)

¡=0.60¡=0.65¡=0.70¡=0.80¡=1.00¡=1.20

1.6

P

0.8

1.3

0.7

1.2

0.6

1.1

0.5
0.1

(d)

1.6
1.5
1.4

0.15

0.2

0.25

¡=0.60¡=0.65¡=0.70¡=0.80¡=1.00¡=1.200.1

0.3

T¡=0.60¡=0.80¡=1.00¡=1.20

(c)

0.9

0.15

(e)

0.2

0.25

0.3

T

1

T
(f)

0.9
0.8

C 1.3

P 0.7

1.2

0.6

1.1

0.5

1

0.4

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8¡=0.60¡=0.80¡=1.00¡=1.200.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

T

T

T

Figure 3.11: ✏ dependence of heat capacity and fraction of 2-in-2-out structures for L = 5(a): Heat capacity in the narrow range of T = 0.0625 to T = 0.3. (b): Fraction of 2-in-2-out structures in
the narrow range of T = 0.0625 to T = 0.3. (c): Incremental di↵erences of fractions of 2-in-2-out structuresin the narrow range of T = 0.0625 to T = 0.3. (d): Heat capacity in the wide range of T = 0.0625 to T = 0.8.
(e): Fraction of 2-in-2-out structures in the wide range of T = 0.0625 to T = 0.8. (f): Incremental di↵erencesof fractions of 2-in-2-out structures for ✏ = 1.0, 1.2 in the wide range of T = 0.0625 to T = 0.8.

Figure 6: Heat capcity C and fraction P of tetrahedra which satisfy the ice-rule. L = 5. (a)-(c)for 0 < T < 0.3 (d)-(f) for a wider temperature range 0 < T < 0.8. (taken from [1])the lower temperatures tend to saturate to a
finite temperature T ∼ 0.16 by increasing .
Very importantly we confirmed that the orbital
degrees of freedom remain dynamically fluctuating between the intermediate temperatures
between the higher and lower peak temperatures of the heat-capacity [2, 3]. This meansthat the orbitals do not produce the putative
static quenched disorder for spins but remain
dynamical. We consider that the spin and orbital degrees of freedom interact cooperatively
at the energy scale corresponding to the lower
peak temperature and exhibit the simultaneous glass transition at that energy scale.

The criticality of the spin-orbital glass transition and its universality, consistency with theexperimental results [19, 20, 25], should be examined by approaching more closer to the critical temperatures.
Acknowledgments This work was supported by KAKENHI (No.
19H01812,
17K05533,
17K05497,18H01173,17H02916)
from MEXT, Japan. The computation in thiswork has been done using system B of theSupercomputer Center, the Institute for Solid
State Physics, the University of Tokyo.

5[1] Kota Mitsumoto, Chisa Hotta, and Hajime Yoshino. Physical Review Letters,
124(8):087201, 2020.

References

Summary and Outlook

To summarize we constructed a theoretical model to describe the spontaneous glasstransition observed in the pyrochlore oxide
Y2 Mo2 O7 and performed extensive Monte
Carlo simulations of the model. Our results
strongly suggest a new type of glass transition- spin-orbital glass transition.
There are numerous possibilities for further
works. Exploration of the parameter space,especially along the -axis, should be done
to clarify to what extent our present scenario
holds. Construction of a mean-field theory [37]would provide a useful guideline in this respect.

[2] Kota Mitsumoto, doctoral thesis, Osaka
Univ., 2020.
[3] Kota Mitsumoto, Chisa Hotta, and Hajime Yoshino. in preparation.
[4] John A Mydosh. Spin glasses: an experimental introduction. CRC Press, 2014.
[5] M Mézard, G Parisi, and M Virasoro. volume 9. World Scientific Publishing Company, 1987.

41

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[6] K Binder and AP Young. Reviews of Modern physics, 58(4):801, 1986.

[21] JS Gardner and et al. Physical review letters, 83(1):211, 1999.

[7] H Kawamura and T Taniguchi.
ume 24. Elsevier, 2015.

[22] HD Zhou, CR Wiebe, A Harter, NS Dalal,
and JS Gardner. Journal of Physics: Condensed Matter, 20(32):325201, 2008.

vol-[8] PW Anderson. Science, 267:1615, 1995.
[23] PMM Thygesen and et al. Physical review
letters, 118(6):067201, 2017.

[9] W. Kauzmann. Chem Rev, 43:219–256,
1948.

[24] CH Booth and et al. Physical Review B,
62(2):R755, 2000.

[10] Jean-Pierre Hansen and Ian R McDonald.
Theory of simple liquids. Elsevier, 1990.

[25] N Hanasaki and et al. Physical review letters, 99(8):086401, 2007.

[11] M. Mézard, G. Parisi, and M. A. Virasoro.
Spin glass theory and beyond. World Scientific, Singapore, 1987.

[26] JN Reimers, AJ Berlinsky, and AC Shi.
Physical Review B, 43(1):865, 1991.

[12] Giorgio Parisi, Pierfrancesco Urbani, andFrancesco Zamponi. Theory of Simple
Glasses: Exact Solutions in Infinite Dimensions. Cambridge University Press,
2020.

[27] R Moessner and JT Chalker. Physical review letters, 80(13):2929, 1998.
[28] TE Saunders and JT Chalker. Physicalreview letters, 98(15):157201, 2007.
[29] H Shinaoka, Y Tomita, and Y Motome.
Physical review letters, 107(4):047204,
2011.

[13] Ludovic Berthier and Giulio Biroli. Reviews of Modern Physics, 83(2):587, 2011.
[14] Yoshihiko
NishikawaandKoji
Hukushima.
arXiv
preprint
arXiv:2003.02872, 2020.

[30] L Pauling. volume 260. Cornell university
press Ithaca, NY, 1960.

SciPost Phys, 4:040,[31] ST Bramwell and MJP Gingras. Science,
294(5546):1495–1501, 2001.

[16] J. G. Kirkwood. in in Phase Transitions
in Solids, edited by R. Smoluchowski, J.
E. Mayer, and A. Weyl (Wiley, New
York, 1951).

[32] Roger G Melko, Byron C den Hertog, andMichel JP Gingras. Physical review letters, 87(6):067203, 2001.

[15] Hajime Yoshino.
2018.

[33] S Große Pawig and K Pinn. International
Journal of Modern Physics C, 9(05):727–
736, 1998.

[17] Leticia F Cugliandolo, Laura Foini, andMarco Tarzia.
Physical Review B,
101(14):144413, 2020.

[34] JL Alonso and et al. Physical Review B,
53:2537–2545, 1996.

[18] Gilles Tarjus, Steven A Kivelson, Z Nussinov, and Pascal Viot. Journal of Physics:
Condensed Matter, 17(50):R1143, 2005.

[35] K Hukushima and K Nemoto. Journal of
the Physical Society of Japan, 65(6):1604–
1608, 1996.

[19] MJP Gingras, CV Stager, BD Gaulin,
NP Raju, and JE Greedan. Journal of
applied physics, 79(8):6170–6172, 1996.

[36] NP Raju, E Gmelin, and RK Kremer.
Physical Review B, 46(9):5405, 1992.
[37] Kota Mitsumoto and Hajime Yoshino. in
preparation.

[20] MJP Gingras, CV Stager, NP Raju,
BD Gaulin, and JE Greedan. Physical review letters, 78(5):947, 1997.

42