tutorial on

Semi-classical Monte Carlo and molecular dynamics
simulations of spin coherent states

Rico Pohler.pohle@aion.t.u-tokyo.ac.jp

第19回 HPC-Phys 勉強会@ RIKEN

Aug. 31, 2023

More is different …

P. W. Anderson, Science 177, 393-396 (1972)

Taisen Deshimaru
弟子丸 泰仙Shibuya at Halloween(after COVID)

8E

6single electron
analytically exact

1023 electrons in real crystal
with complex intertactions,impurities, distortions, …

420-4-20x2

How do so many electrons
collectively behave?

42Monte Carlo …. does it work?

Neutron scattering in pyrochlore
magnet Yb2Ti2O7
with spin anisotropic interactions
MCSF

Neutron scattering in diamondspiral spin liquid MnSc2S4

Neutron scattering fortriangular magnet
YbMgGaO4

Magnetic susceptibility for thecentered pyrochlore magnet[Mn(II)(ta)2]

J. A. M. Paddison et al.,Nature Physics 13, 117-122 (2017)

R. P. Nutakki et al., PRR 5, L022018 (2023)

NSF
INS

MC

L.D.C. Jaubert et al.,PRL 115, 267208 (2015)

INS

SF

NSF

L.J. Chang et al.,Nature Commun 3, 992 (2012)

S. Gao et al.,Nature Physics 13, 157-161(2017)

3

R. P. Nutakki et al., arXiv:2203.08780v1

Monte Carlo …. does it work?

Neutron scattering of Ca10Cr7O28 compared to simulations on theBilayer- breathing Kagome lattice model

Neutron scattering of NaCaNi2F7 compared tosimulations on the Heisenberg model

INS

MD

MD

R. Pohle et al., PRB 104, 024426 (2021)

S. Zhang et al., PRL 122, 167203 (2019)4

Advantages

Disadvantages‣ versality of models
applicable for any type of lattice in any dimension(1D, 2D, 3D / Ising, XY, Heisenberg / any type of interactions /impurities / eld / …)

‣ Not an approximation‣ no quantum entanglement

numerically exact

classical evaluation of local energies‣ versality of measures
measure any physical observable‣ thermalization(magnetic order parameter, structure factor, speci cheat, susceptibilit, spin dynamics)

system may fall out of equilibrium‣ large system size
> 10,000 sites( nite size scaling is easy)‣nite temperaturescomparison to experimentbuild intuition for physical system
fifi

5fifi

Why is it a useful tool?

Monte Carlo sampling - a basic example

Book recommndation:

Monte Carlo Methods in
Statistical Physics
M. Newman, G. Barkema

area circle
πr 2 π N◯= 2 = =
area square
4r
4N□‣ evaluate condition numericallyif

22x +y ≤1N◯ = N◯ + 1
thenN□ = N□ + 1

redN□ = N□ + 1

blueelse‣ Central limit theorem
relative error ~ 1/ N

‣ especially good for multi-dimensional integration
6

Essence of Monte Carlo simulations
Master Equation:dωμ(t)dt=ων(t)P(ν → μ) − ωμ(t)P(μ → ν)][∑ωμ … occupation probabilityP(μ → ν) … transition probability

νP(μ → ν)‣ Markov process:system has no memory‣ Ergodicity:

access whole phase space in nite amount of time

‣ Detailed balance:systems in equilibriumpνP(ν → μ) = pμP(μ → ν)fi

7

is time independent

Finite-temperature Monte Carlo simulations‣ Detailed balance:

56

II Methods

Initial spin
configurationpνP(ν → μ) = pμP(μ → ν)MC timet = 0
MC moven = 0make connection between statistical mechanics and thermodynamics‣ Metropolis algorithm
N. Metropolis and S. Ulam, JASA 44 (247): 335–341 (1949)

canonical ensemble
coupled to external heat-bathSpin σt,n chosen randomly
Compute ∆E = Ef − Ei
Generate random r ∈ [0 : 1]n = n+1

t= t+1

N. Metropolis et al., J. Chem. Phys. 21, 1087–1092 (1953)
citation count > 50800

Boltzmann weightYes

r < e−β∆E

Spin flip of σt,nZ … partition function

NoStore the new
configuration

β = 1/kBT

ratio of transition probabilities

Yes

Yes

n<NNot < Teq(detailed balance)

No
Compute observables
Outputdata

most ef cient transitionfi8

Average Noobservables

L. D. C. Jaubert, PhD thesist < Ttot

YesStore dataproblem dependent spin updates

56‣ local spin updates:‣ non-local (cluster updates updates)• single spin ip• Swenden-Wang / Wolf algorithm

acceptence ratio usually T dependent

R. H. Swendsen and J. S. Wang, PRL 58, 86 (1987)
U. Wolf, PRL 62, 361 (1989)

G. Marsaglia, Ann. Math. Stat. 43, 645 (1972).

II Methods

Initial spin
configuration

MC timet = 0
MC moven = 0Spin σt,n chosen randomly
Compute ∆E = Ef − Ei
Generate random r ∈ [0 : 1]

A. Barbu, Computer Vision (pp 783 - 785), Springer• heat bath updaten = n+1

t= t+1Yes

Miyatake et al, J. Phys. C:Solid State Phys. 19, 2539-2546 (1986)

rejection free

Spin flip of σt,nNo

reduce critical scaling near phase transition• over relaxationr < e−β∆EStore the new
configuration• Loop / Worm updates

Yes

Yesn<NNot < Teq

No
Compute observables

energy conserving, decorrelate systemchange between topological sectors

flM. Creutz, PRD 36, 515 (1987)L. D. C. Jaubert et al, PRX 3, 011014 (2013)
H. Shinaoka et al, JPCS 320 (2011) 012009

Outputdata

9Average Noobservables

L. D. C. Jaubert, PhD thesist < Ttot

YesStore dataFinite-temperature Monte Carlo simulations

56

II Methods

Initial spin
configuration‣ minimize local bond energies for any type of lattice

MC timet = 0
MC moven = 0Spin σt,n chosen randomly
Compute ∆E = Ef − Ei
Generate random r ∈ [0 : 1]….
n = n+1

t= t+1

Yesr < e−β∆E

Spin flip of σt,nNNoN

Store the new
configuration

H=Siα JijαβSjβ −Siαhiα∑ ∑∑ ∑i,j α,β=xyz
i α=xyz

Yes

Yesn<N

spin con guration which respects Boltzmann distribution

Not < Teq

No
Compute observables
Output
data

fi

10

Average Noobservables

L. D. C. Jaubert, PhD thesist < Ttot

Yes

Store databrute force thermalization

56
Initial spin
configuration‣ Parallel tempering

MC timet = 0

replica exchange method

MC moven = 0hotcoldOpen MPIII Methods

Spin σt,n chosen randomly

CPU 0 CPU 1 CPU 2 CPU 3 CPU 4 CPU 5 CPU 6 CPU 7 CPU 8

T0T1

T2

T3T4T5T6T7

Compute ∆E = Ef − Ei
Generate random r ∈ [0 : 1]

T8Monte Carlo timet0n = n+1

t= t+1

Yest1t2

Spin flip of σt,nNo

t3

Store the new
configuration

t4

Yes

probability
distribution functionr < e−β∆E

Yes

n<NNo

t < TeqNooverlap

Compute observables
Outputdata

11

Average No
observables

L. D. C. Jaubert, PhD thesist < Ttot

YesStore dataspin order parameter space‣ local degrees of freedom for spinsdescrete
Isingcontinuous

q-state clock/Potts

2π
ϕ=q

σ =±1cosϕS ⃗=( sin ϕ )2π 4π 6π
ϕ = 0,,,,⋯q q q

XY=O(2)O(3) = CP1

CP2

CP(N-1)

SU(2)

SU(3)

SU(N)

ϕcos
ϕS ⃗=( sin ϕ )sin θ cos ϕ
S ⃗ = sin θ sin ϕcos θ0 ≤ ϕ < 2π0 ≤ ϕ < 2π0≤θ≤π12retemaraprofgn1 ??
inePpCpanahhstitreagraWhlecapsspin order parameter space‣ generalization to SU(N)spin coherent state| Ωi⟩ =N∑
α=1cα(Ωi) | α⟩

Ω = {ξ1, ⋯, ξN−1, ϕ1, ⋯, ϕN−1}
ϕ0 = 0ξN = 0R. Iwazaki et al., arXiv:2301.09824
CP1CP2CP(N-1)

O(3) = SU(2)

SU(3)

SU(N)

1+= c† 1+22⟩⟩i

1−= c† 1−22⟩⟩iiicoef cients

Degree of freedom: 2(N-1)

observables

dipolesquadrupoles

Siα = ⟨Ω | Sα̂ | Ω⟩i

Qiα,β = ⟨Ω | Sα̂ Sβ̂ | Ω⟩i

1 0 1x̂S =2 (1 0)

1 0 −iŷS =2 (i 0 )1 1 0
S =2 (0 −1)𝖼ẑ

𝖼𝖺𝖺𝖼𝖼𝖺𝖺𝗏

𝗏

𝖼𝖺𝗏
fi𝗏

𝗏

α,β,γ
T= ⟨Ω | Sα̂ Sβ̂ Sγ̂ | Ω⟩ioctupoles ic1(Ω) = cos ξ1c2(Ω) = e iϕ1 sin ξ113†+1⟩ = c+1i0⟩ = c0†i†−1⟩ = c−1i

⟩i
⟩i

⟩i

c1(Ω) = cos ξ1c2(Ω) = e iϕ1 cos ξ2 sin ξ1c3(Ω) = e iϕ2 sin ξ1 sin ξ20 1 01 0 12 (0 1 0)

1x̂

S =ŷ

S =i

1+ (N − 1)⟩2i⋮
1− (N − 1)
⟩20 −1 0
1 0 −12 (0 10)

1 0 0Sẑ = 0 0 0(0 0 −1)

icα(Ω) = e iϕα−1 cos ξαaj =x̂

S =120a1Sẑ =∏β=12(S + 1)j − j( j + 1)

a10⋱

0 −ia1
ia1 0ŷ
S =
⋱Sα−1

⋱
⋱ aN−1aN−1 0⋱⋱ −iaN−1
iaN−10S−1⋱−Ssin ξβ

spin order parameter space‣ generalization to SU(N)H=

HamiltonianN∑∑i,j α,βLandau-Lifschitz dynamics
Molecular dynamics

for SU(2) case, we recoverNα αβi Jijd αi∂H= fαβγ βdt∂ iβ−j∑∑i
αα αi hî α | Ω⟩i

̂ α, ̂ β = i f̂γαβγ i[ i i]γidSi∂H= − Si ×dt∂Siαi = ⟨Ω |with fαβγ = ϵαβγ

generalized spin components, fromgenerators of corresponding SU(N)dipoles, qadrupoles, octupoles,…

structure constants fαβγ

recent papers:
H. Zhang and C. Batista, PRB 104, 104409 (2021)
D. Dahlbom et al., PRB 106, 054423 (2022)
D. Dahlbom et al., PRB 106, 235154 (2022)‣ solve by numerical integration

K. Remund et al., PRR 4, 033106 (2022)
R. Iwazaki et al., arXiv:2301.09824

𝒪𝒪𝒪𝒪𝒪𝒪𝒪𝒪𝒪𝒪𝒪D. Dahlbom et al., arXiv: 2304.03874

14

numerical integration

Book recommndation:

numerically solve differential equations:dy= f(t, y)dty(t0) = y0

Euler method

mid-point method

Runge-Kutta 1st order

Runge-Kutta 2nd orderk1 = f(tn, yn)hhk2 = f tn + , yn + k1
(22 )yn+1 = yn + h f(tn, yn)yn+1 = yn + h k2

tn+1 = tn + h

tn+1 = tn + hexact solution, y(t)y(t)
tangent line

yn+1
yn

f(tn, yn)

tn

tn+1t

Solving Ordinary Differnetial
Equations I: Nonstiff Problems
E. Hairer, S.P. Norsett, G. Wanner

Runge-Kutta 4th orderk1 = f(tn, yn)hhk2 = f tn + , yn + k1
(22 )hh
k3 = f tn + , yn + k2
(22 )k4 = f (tn + h, yn + k3)

hyn+1 = yn + (k1 + 2k2 + 2k3 + k4)
6tn+1 = tn + hexact solution, y(t)

y(t)

yn+1

tangent line

hhf tn + , yn + k1
(22 )ynerror estimate:
local truncation error:O(h5)total accumulated errror : O(h4)tntn + h2
15

tn+1tnumerical integration
t/2

y(t) = e sin(5t)one example:

dy 1 t/2= e sin(5t) + e t/2 cos(5t)dt
2Euler method

mid-point method

Runge-Kutta 1st order

Runge-Kutta 2nd ordery(0) = 0
Runge-Kutta 4th orderk1 = f(tn, yn)

hhk2 = f tn + , yn + k1
(22 )hh
k3 = f tn + , yn + k2(22 )k1 = f(tn, yn)

yn+1 = yn + h k2

tn+1 = tn + h

tn+1 = tn + h

h=0.5

1

23
t4

5

6

15
10
5
0-5
-10
-15
-20
0

k4 = f (tn + h, yn + k3)

hyn+1 = yn + (k1 + 2k2 + 2k3 + k4)
6tn+1 = tn + h

h=0.5

y

15
10
50
-5
-10
-15
-20
0

yn+1 = yn + h f(tn, yn)yyhhk2 = f tn + , yn + k1(22 )

123t
16

45615
1050-5-10-15-200h=0.5

123t456

Molecular Dynamics simulations - some examples
Case for SU(2) spin coherent states on the
S=1/2, Bilayer Breathing Kagome lattice

Dynamics for Ca10Cr7O28

https://www.youtube.com/watch?v=GGug4w4MiKQ
RP, H. Yan, N. Shannon, PRB 104, 024426 (2021)
17

Molecular Dynamics simulations - some examples
Case for SU(3) spin coherent states on the
S=1, Bilinear Biquadratic model

K. Remund, R.P., Y. Akagi, J. Romhányi, N. Shannon,Phys.
Rev.
Research4,u3MC033106(2022)(a) Finite-Tphase
diagram fromsimulation

4Figure 3. Finite–temperature phase diagram of the spin–1 bilinear–biquadratic (BBQ) model on the triangular lattice, obtained from Monte
Carlo simulation of HBBQ [Eq. (9)] in the space of u(3) matrices (u3MC). Points show the location of peaks in heat capacity; phases arelabelled according to their dominant correlations (cf. Fig. 2 and Fig. 8). Details of simulation are provided in Section III A 2.
gravitational waves in spin nematics?!

to which quantum results for dynamics can be inferred from(semi–)classical simulations.
In order to keep the paper self–contained, we provide a
detailed account of derivations, and review all of the necessary mathematical formalism. However the paper is also constructed in such a way that readers uninterested in technicaldevelopment of the method can skip directly to the results provided in Section VI, Section VII and Section VIII. Additional
technical details are provided in a series of Appendices.
The paper is structured as follows:
Section II reviews the mathematical formalism needed toanalyze spin–1 magnets in terms of a u(3) algebra. A single spin–1 moment is analysed within a basis of time–reversal
invariant states, where the most general possible spin configuration can be expressed in terms of a complex vector d.
The group U (3) is shown to provide a convenient basis for
all possible operations on spin–1 moments. Expressions are
given for both the dipole
18 and quadrupole associated with aspin–1 moment, in terms of the vector d, and the 3 ⇥ 3 ma-spin–1 BBQ model on a triangular lattice, starting from theu(3) representation introduced in Section II. In Section IV A,
small fluctuations about FQ order are recast in terms of thethe matrices Â↵ . In Section IV B these are shown to provide a natural basis for a classical, low–temperature expansion, within which it is possible to calculate thermodynamicproperties. Finally, in Section IV C we use this theory to makeexplicit predictions for thermodynamic properties, including
structure factors, for later comparison with MC simulation.
In Section V, we develop an equivalent, zero–temperature,quantum theory of fluctuations about a FQ ground state. First,in Section V A we quantize the fluctuations introduced in Section IV A, and show that the resulting multiple–Boson expansion is equivalent to earlier “flavor–wave” theory. Then, inSection V B, we use this quantum theory to make explicit predictions for dynamical structure factors, for later comparison
with MD simulation.
In Section VI, the numerical methods developed in Sec-other interesting examples
Case for SU(3) spin coherent states:
S=1, CP2 skyrmions

Y. Amari et al., PRB 106, L100406 (2022)

H. Zhang et al., Nat Comm 14, 3626 (2023)
19other interesting examples

R.P., N. Shannon, Y. Motome, PRB 107, L140403 (2023)

SU(3) spins: CP2

Case for SU(3) spin coherent states on the
S=1, Bilinear Biquadratic model with
Kitaev interactions

L=120, N=28,800
FMchiral

NC 2y

SU(3)

FQ

SU(3)

q1D 1

AFM

AFMzz

FMSO
SU(3)

SU(3)

ystripyTC
q1D 2xFM Kitaev SLTC

AFM
chiralNC 1xzigzag
AFM
Kitaev SLO(3) spins: CP1
FM CSLCPl 2

NC 1

FM Kitaev SL

NC 2

CPl 1

α=x,y,z ⟨ij⟩α

1D 1⟨i,j⟩∑ ∑SiαSjα

AFM CSL∑2JS⋅
S+J(S⋅S)
2ij ] +K
[ 1 i jAFMℋ=FM
chiral

1D 2

L=120, N=28,800

stripyFM

AFM
AFM
chiralzigzag
AFM
Kitaev SL20