3.3

Strongly Correlated

Quantum Systems

196

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

Mechanism of quantum spin liquidand high-Tc superconductivity
studied by excitation spectra and nonequilibrium dynamics
calculation by variational Monte Carlo methods
Masatoshi IMADA
Toyota Physical and Chemical Research Institute41-1 Yokomichi, Nagakute, Aichi, 480-1192
andResearch Institute for Sicence and Engineering, Waseda University
3-4-1 Okubo, Shinjuku-ku, Tokyo, 169-8555
Understanding physics of high-Tc cupratemodel, namely Hubbard model on the squaresuperconductors remains one of the importantlattice

problems in materials science. In this project, wesuperconducting state is severely competing

have

phase

with the stripe ordered state, while most of the

diagram of HgBa2CuO4+y [1] by solving its abphase diagram in the parameter space of theinitio low-energy effective Hamiltonian without

carrier doping concentration is dominated by the

adjustable parameters [2].

It shows astripe order and the superconducting statesuperconducting phase in a wide range of holeremains an excited state [3,4]. Therefore, thedensity δ, and its competition with charge

Hubbard model does not properly account forperiod-4 plus spin period-8 stripe order near δ ∼the experimental phase diagram of the cuprate

0.1, in agreement with experimental results.

superconductors. However, it has also been

reproduced

theexperimental

hadclarified

that

the

d-waveA crucial role of off-site interactions inshown that the superconducting phase in the abstabilizing the superconductivity is elucidated

initio Hamiltonian is adiabatically connectedwith emphasis on charge fluctuations. It alsowith that found in the superconducting excitedclarifies

mainlystate of the Hubbard model. Therefore, it is

Coulomb

intriguing to understand the distinction of theinteraction. The present achievement will enablesuperconducting state in between the ab initio

deeper, predictable understanding on openand Hubbard cases. The method for dynamicalsuperconducting

properties newly developed offer useful insights

mechanism and promote ab initio studies on

into this issue as we describe in the next section.

strongly

Calculations of dynamical quantities are athecontributed

issuesofcondensation
from

the

the

on-sitehigh-Tc

correlated

energyelectrons

beyondparametrized model studies.

challenging subject for strongly correlated

In former studies on a simplified effective

electron systems even in the linear response

197

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

regime. Recently we have developed a method

method. The application of the method to larger-

to calculate the dynamical correlation functions

sized Hubbard model on the square lattice

based on the variational Monte Carlo method

correctly[3,4]. We first developed a method to calculate

behavior at half filling and gap structures of d -the charge dynamical structure factors for thewave superconducting state of the hole doped

ground states of correlated electron systems. Our

Hubbard model in the ground state optimized bybenchmarks for the one- and two-dimensional

enforcing the charge uniformity, evidencing a

Hubbard models show that inclusion ofwide applicability to strongly correlated electroncomposite-fermion excitations in the basis set

systems.

reproduces

the

Mott

insulating

greatly improves the accuracy, in reference to

From the obtained d-wave superconducting

the exact charge dynamical structure factors forgap of the charge uniform state, we find that theclusters. Together with examination for larger

gap amplitude at the antinodal point is several

systems beyond tractable sizes by the exact

times larger than the experimental value, when

diagonalization, our results indicate that thewe employ a realistic parameter as a model ofvariational Monte Carlo method is a promising

the cuprate superconductors. The effective

way for studies on the nature of charge dynamics

attractive interaction of carriers in the d-wave

in correlated materials such as the copper oxide

superconducting state inferred for an optimized

superconductors

composite-fermion

state of the Hubbard model is as large as theexcitations are properly included in the restricted

order of the nearest-neighbor transfer, which isHilbert space of intermediate states in the linear

far beyond the former expectation in the

response theory. Our results on the charge

cuprates. The nature of the superconducting

dynamical structure factor are consistent withstate of the Hubbard model in terms of thethe particle-hole excitations inferred from theoverestimate of the gap and the attractive

single-particle spectral function A(k,ω) in theinteraction in comparison to the cuprates are

literature. The importance of incorporating

clarified [4].

nonlocal composite fermion for a more accurate

Finding quantum spin liquids is anotherdescription is also pointed out [3].

important subject of the present research project.

This method can be extended to calculate theTopological properties of the strongly correlated

one-body Green's function by improving thesystems with strong spin-orbit interaction arevariational Monte Carlo method. We benchmark

one of such subjects. We have studied

against the exact diagonalization for the one-

Heisenberg model with spin-orbit interaction

and two-dimensional Hubbard models of 16 site

(Dyaloshinskii-Moriya

lattices, which proves high accuracy of the

pyrochlore lattice by using the variational Monte

if

the

198

interaction)

on

the

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

Carlo method. This system has been proposed to

conductors

be a candidate of the quantum spin liquid in

Hamiltonian”.

earlier studies. We will report the result on thisand “Studies on Quantum Spin Liquids in

based

on

first-principles

Materials with Strong Spin-Orbit Interaction”

subject in the future.

This series of work has been done in

Another candidate of the quantum spin liquid

collaboration with T. Ohgoe, T. Misawa, M.

is found in a quantum Heisenberg model with

Hirayama, K. Ido, Y. Nomura, Y. Yamaji, R.

the geometrically frustrated square lattice, called

Pohle and M. Charlebois. The calculations were

J1-J2 model. By utilizing the machine learning

performed on computers at the Supercomputer

method, we have found a firm evidence for theCenter, ISSP, Univ. of Tokyo. The work is also

quantum spin liquid and clarified its excitation

supported by JSPS Kakenhi 16H06345, thespectra, which indicate the fractionalization of

RIKEN Center for Computational Science underspins into spinons. The subject has a closethe HPCI project (hp190145) and the MEXTHPCI Strategic Programs, CDMSI. This project

connection to the experimentally observed spin

has used the software HPhi and mVMC.

liquid state in molecular conductors. Studies on

References

the ab inito low-energy effective Hamiltonian[1] Motoaki Hirayama, Takahiro Misawa,
for a dmit compound are underway and will be
Takahiro Ohgoe, Youhei Yamaji, Masatoshi
reported elsewhere.

Imada, Phys. Rev. B 99 (2019) 245155.

This is a combined report for E project[2] Takahiro“Mechanism of quantum spin liquid and high-Tc

Ohgoe,

Motoaki

Hirayama,

Takahiro Misawa, Kota Ido, Youhei Yamaji,

superconductivity studied by excitation spectraand nonequilibrium dynamics calculation by

Masatoshi Imada, Phys. Rev. B 101 (2020)

variational

045124.

Monte

Carlo

methods”,“Development of methodology for analysis of[3] Kota Ido, Masatoshi Imada, Takahiro

nonequilibrium superconductivity in strongly

Misawa, Phys. Rev. B 101 (2020) 075124.

correlated systems via integration of electronic[4] Maxime Charlebois, and Masatoshi Imada,state theory and time-resolved experimental

arXiv:1912.09960.

data”, and shared project for post-K project“Studies on quantum spin liquids in molecular

199

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

Numerical studies on energy current correlations infrustrated quantum spin systems
Youhei YAMAJI
Department of Applied Physics, The University of Tokyo
Hongo 7-3-1, Bunkyo-ku, Tokyo 113-8656

Topological states of matters are often characterized by quantized transport properties.
As a prototypical example, two dimensional
Chern insulators are characterized by thequantized charge Hall eﬀect. In the Mott insulators, however, the quantization of thermal
Hall transports characterizes the topologicalstates, instead of the charge transports [1, 2].
The quantum spin liquid state of the Kitaev
model is a typical example of the topologicalstate that shows half-integer quantized thermal Hall transports [1]. The observation of thehalf-integer quantized thermal Hall coeﬃcientin a spin-orbit coupled Mott insulator α-RuCl3
has attracted broad attention [3].
We studied numerical methods to characterize the topological spin liquids. The thermal transport is calculated by correlations ofthermal currents and energy magnetization [4].
These quantities are formulated by energy current operators, which are given by the continuity equation of the energy density.
For the lattice hamiltonian, the energy density is represented by the bond energy as follows. Here, we show the continuity equation of
the energy density of the Kitaev-type magnets
on the honeycomb lattice as an example. The
nearest-neighbor bonds on the honeycomb lattice have three diﬀerent directions. When thethree bonds are labeled as x, y, and z, an effective hamiltonian of α-RuCl3 , the Kitaev-Γ

hamiltonian,
Ĥ =

∑

∑

(γ)

Ĥℓm ,(1)

γ=x,y,z ⟨ℓ,m⟩∈γ

is defined by the exchange coupling for the γ(= x, y, z) bond,
(
)
(γ)
γ
β
α
Ĥℓm = K Ŝℓγ Ŝm, +Γ Ŝℓα Ŝm, +Ŝℓβ Ŝm,(2)where (α, β, γ) is a permutation of (x, y, z) andK (Γ) is the Kitaev (oﬀ-diagonal symmetric)
exchange coupling. The bond energy of thehamiltonian is, then, introduced as the Fouriercomponents,
∑ ∑
Ĥq =
e+i(q /2)·(Rℓ +Rm )γ ⟨ℓ,m⟩γ

[
(
)]
γ
β
α
× K Ŝℓγ Ŝm+ Γ Ŝℓα Ŝm, +Ŝℓβ Ŝm, ,(3)
where ϵαβγ (α, β, γ = x, y, z). By using thesmall q limit of the continuity equation of thebond energy,dĤq∼ −iq · j E ,dt(4)

we determine the energy current operator j E .
When, for simplicity, K = 0 is chosen, theHeisenberg equation of motion,
]
dĤq
i[
=
Ĥ, Ĥq ,dtℏ
defines j E with the commutator,

200

(5)

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo]i[
Ĥ, Ĥq =
ℏ

(
)
i 2∑[
y y
y
x x
z
z y
y
Γ
q · δ xz,ℓ Ŝm
Ŝ
Ŝ
+
Ŝ
Ŝ
Ŝ
−
Ŝ
Ŝ
Ŝ
m
m
m
m
m
ℓ
ℓ
x
z
x
z
x ℓ
z
ℏ
ℓ
(
)
x
y
x
x
+q · δ zy,ℓ ŜmŜ z Ŝ z + Ŝm
Ŝ y Ŝ x − ŜmŜ x Ŝmz ℓ myz ℓ my
z ℓ
y
(
)]
y y
z
x
x z
z
z z
+ q · δ yx,ℓ Ŝm
Ŝ
Ŝ
+
Ŝ
Ŝ
Ŝ
−
Ŝ
Ŝ
Ŝ
my ℓ mx
my ℓ mxy ℓ mx
[
(
)
∑
i
x x
z y
y y
q · δ xz,m Ŝℓyx Ŝm
Ŝℓz + Ŝℓzx Ŝm
Ŝℓz − Ŝℓyx Ŝm
Ŝℓz
+ Γ2
ℏ
m
(
)
z z
y x
x x
+q · δ zy,m Ŝℓxz ŜmŜℓy + Ŝℓyz ŜmŜℓy − Ŝℓxz Ŝm
Ŝℓy
(
)]
y y
x z
z z
+ q · δ yx,m Ŝℓzy ŜmŜℓx + Ŝℓxy ŜmŜℓx − Ŝℓzy ŜmŜℓx ,(6)

where δ αβ,ℓ is the vector from the center of the α-bond terminated by ℓth site and to the centerof the β-bond terminated also by ℓth site.

The energy-current–energy-current correlations and energy magnetizaion for the KitaevΓ hamiltonian at finite temperatures are calculated by combining the Krylov subspace methods and typicalit approaches [5] and will bereported elsewhere [6].

References

The topological nature of the Mott insulating states will be destroyed by spontaneous
symmetry breakings. Therefore, to clarify the
topology, we need to examine whether the spin
liquid candidates show spontaneous symmetry
breakings or not. We have modified mVMC [7]to calculate higher-order spin correlations andcalculated the singlet-singlet correlatons of atypical frustrated quantum spin systems, theHeisenberg model on pyrochlore lattices. So
far, the inversion symmetry is spontaneously
broken in the lowest variational state.

[3] Y. Kasahara, et al.: Nature 559 (2018)
227.

[1] A. Kitaev: Ann. Phys. 321 (2006) 2.
[2] C. Wang, A. C. Potter, and T. Senthil:
Science 343 (2014) 6171.

[4] T. Qin, Q. Niu, and J. Shi: Phys. Rev.
Lett. 107 (2011) 236601.
[5] Y. Yamaji, T. Suzuki, and M. Kawamura:
arXiv:1802.02854.
[6] Y. Yamaji, et al.: in preparation.
[7] T. Misawa, et al.: Compt. Phys. Commun. 235 (2019) 447.

201

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

Theoretical study of spin-orbit coupled phenomena
in correlated electron systems
Yukitoshi MOTOME
Department of Applied Physics,
The University of Tokyo, Bunkyo, Tokyo 113-8656

We have theoretically studied a variety ofthis topic in review articles [11,12,13].

intriguing phenomena in correlated electron(ii) Chiral magnets: We theoretically unveiled

systems with spin-orbit coupling, ranging fromthat a nonreciprocal spin transport, which is

Mott insulators to metals. During this fiscal

switchable by a magnetic field, can occur inyear, we have achieved substantial progress on

monoaxial magnets [14]. We also clarified thatthe following topics (project numbers: H31-

interplay between spin, charge, and orbital

Ca-0052 and H31-Cb-0032). We summarize

degrees of freedom plays a crucial role inthe main results below.

stabilizing the magnetic hedgehog lattices(i) Kitaev quantum spin liquids: By using ab

recently found in B20 chiral magnets [15]. We

initio based calculations, we theoreticallyalso showed that the magnetic monopoles andproposed new candidates for the Kitaev

antimonopoles associated with the magnetic

magnets with unconventional antiferromagnetic

hedgehogs annihilate in pair by magnetic fields

Kitaevand cause topological phase transitions [15,16].

interactions:Pr-basedf-electroncompounds [1,2] and polar asymmetric systems(iii)[3]. We also clarified real-time dynamics ofproposed that the successive phase transitions

fractional excitations in the Kitaev quantum

in a spin-orbit coupled metal Cd2Re2O7 arespin liquids [4], Majorana-magnon crossover indriven by electric toroidal quadrupoles [17]. We

an applied magnetic field [5], and exotic phasealso clarified that nonmagnetic impurities cantransitions in three-dimensional extensions of

induce magnetic and toroidal multipoles in

the Kitaev models [6,7]. In addition, we haveferromagnets [18].

collaborated with several experimental groups(iv) Spin current generation: We discovered afor identifying the signatures of Majorana

new mechanism of spin current generation,

excitations in the Kitaev magnets: nuclear

which can work even in the absence of thequadrupole and magnetic resonances [8],spin-orbit coupling [19]. We also developed theRaman scattering [9], and thermal Hall effectunderstanding of this phenomenon based on the[10]. We have summarized the recent results oncluster multipole picture [20].

202

Multipole

physics:We

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

References

preprint (arXiv:2001.01899).

[1] S.-H. Jang, R. Sano, Y. Kato, and Y.

[11] Y. Motome and J. Nasu, J. Phys. Soc.

Motome, Phys. Rev. B 99, 241106(R) (2019).

Jpn. 89, 012002 (2020), INVITED REVIEW[2] S.-H. Jang, R. Sano, Y. Kato, and Y.

PAPER.

Motome, preprint (arXiv:1912.03422).

[12] 求 幸年，那須 譲治, 固体物理 54, 217[3] Y. Sugita, Y. Kato, and Y. Motome, Phys.

(2019).

Rev. B 101, 100410(R) (2020).

[13] Y. Motome, R. Sano, S.-H. Jang, Y. Sugita,[4] J. Nasu and Y. Motome, Phys. Rev.

and Y. Kato, preprint (arXiv:2001.03731), in

Research 1, 033007 (2019).

press in J. Phys.: Condens. Matter, Special[5] J. Yoshitake, J. Nasu, Y. Kato, and Y.

Issue Article.

Motome, Phys. Rev. B 101, 100408(R) (2020).

[14] S. Okumura, H. Ishizuka, Y. Kato, J. Ohe,

[6] T. Eschmann, P. A. Mishchenko, T. A.

and Y. Motome, Appl. Phys. Lett. 115, 012401

Bojesen, Y. Kato, M. Hermanns, Y. Motome,(2019).

and[15] S. Okumura, S. Hayami, Y. Kato, and Y.

S.

Trebst,

Phys.

Rev.

Research 1,

032011(R) (2019).

Motome, preprint (arXiv:1908.05044), in press[7] P. A. Mishchenko, Y. Kato, K. O'Brien, T. A.

in Phys. Rev. B.

Bojesen, T. Eschmann, M. Hermanns, S. Trebst,[16] S. Okumura, S. Hayami, Y. Kato, and Y.

and Y. Motome, Phys. Rev. B 101, 045118

Motome, JPS Conf. Proc. 30, 011010 (2020).

(2020), Editors' Suggestion.

[17] S. Hayami, Y. Yanagi, H. Kusunose, and Y.

[8] Y. Nagai, T. Jinno, Y. Yoshitake, J. Nasu, Y.

Motome, Phys. Rev. Lett. 122, 147602 (2019),

Motome, M. Itoh, and Y. Shimizu, Phys. Rev.

selected as a Cover Image.

B 101, 020414(R) (2020).

[18] S. Hayami, H. Kusunose, and Y. Motome,[9] Y. Wang, G. B. Osterhoudt, Y. Tian, P.

J. Phys. Soc. Jpn. 88, 063702 (2019).

Lampen-Kelley, A. Banerjee, T. Goldstein, J.

[19] M. Naka, S. Hayami, H. Kusunose, Y.

Yan, J. Knolle, H. Ji, J. Cava, J. Nasu, Y.

Yanagi, Y. Motome, and H. Seo, Nat.

Motome, S. Nagler, D. Mandrus, and K. S.

Commun. 10, 4305 (2019).

Burch, npj Quant. Mat. 5, 14 (2020).

[20] S. Hayami, Y. Yanagi, M. Naka, H. Seo, Y.

[10] T. Yokoi, S. Ma, Y. Kasahara, S. Kasahara,

Motome, and H. Kusunose, JPS Conf. Proc. 30,

T. Shibauchi, N. Kurita, H. Tanaka, J. Nasu, Y.

011149 (2020).

Motome, C. Hickey, S. Trebst, and Y. Matsuda,

203

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

Numerical Studies on the Role of the Incipient Band Played in Two-band
Unconventional Superconductivity
KAZUHIKO KUROKI
Department of Physics, Osaka University
1-1 Machikaneyama, Toyonaka, Osaka, 560-0043, Japanthe two-leg ladder Hubbard model, and made acomparison between the two models.
In Fig.2, we show for t′ = 0 the eigenvalue λplotted against t⊥ for the bilayer and tr for the
two-leg ladder. For each t⊥ or tr , the maximum
eigenvalue is searched by varying the band ﬁlling n.
We perform the calculation for U = 6t and U = 3t.
For U = 6t in the bilayer and for both values of
U in the ladder, λ is maximized when the bonding
band touches the Fermi level, i.e., when the bandis incipient. This is not the case for U = 3t in the(a)1.4
bonding band below EF

1.2

n=1.05

U=6

1
0.8

λ

We have been studying spin-ﬂuctuationmediated superconductivity in two-band systems
with one of the bands intersecting the Fermi level
while the other band (“incipient band”) sits inthe vicinity of the Fermi level. In ref.[1, 2], we
studied the Hubbard model on a two-leg ladderlattice, which is characterized by the nearest
neighbor hoppings in the leg (t) and the rung(tr ) directions, and also the diagonal next nearest
neighbor hoppings t′ (Fig.1).
In 2019 ﬁscal year, we have extended the study
to two dimensional (2D) systems. A simplest extension of the two-leg ladder Hubbard model to twodimensions is the bilayer square lattice Hubbard
model, where two Hubbard models on a square lattice are coupled by vertical hoppings t⊥ along withdiagonal hoppings t′ (Fig.1). Here, t⊥ controls therelative energy between the two bands like tr in
the two-leg ladder. We adopted the ﬂuctuationexchange (FLEX) approximation, and solved thelinearized Eliashberg equation to obtain the eigenvalue at a ﬁxed temperature (T = 0.05t), which
gives a measure for the superconducting transition
temperature. We investigated a large parameter
space of (t⊥ , t′ , n, U ), where n is the band ﬁllingand U is the on-site Hubbard interaction.[3] Simultaneously, we also performed similar calculation for

n=1.15

0.6

U=3

0.4
n=1.15

0.2
0

1

1.5

n=1.072

2.5

3

bonding band
below EF

3.5

4

t
(b)

bilayer

1.4

2 leg ladder

1.2

t

t
tr´

1

n=1.1

bonding band below EF

U=6

n=1.05bonding band below EF

U=3t
unit cell

λ

0.8

t
´

n=1.050.6
0.4

t

n=1.05

0.2
0
1

EF

1.2

1.4

1.6

1.8

2

2.2

tr
t~t , t’~0

t~t , t’ >0

t<t , t’~0

FIG. 1. Upper panel : the bilayer lattice (left) and thetwo-leg ladder lattice (right). Lower panels : schematic
images of the bonding and antibonding bands of thebilayer lattice. Left : t ∼ t⊥ , t′ ∼ 0, middle : t ∼ t⊥ ,t′ > 0, right : t < t⊥ , t′ ∼ 0. In the middle and rightpanels, the bonding band is incipient.

204

FIG. 2. (a) Maximized λ at T = 0.05 of the bilayer
Hubbard model with U = 6t or 3t and t′ = 0 plottedagainst t⊥ . λ is maximized at each t⊥ by varying theband filling n. At some points, n that maximizes λis denoted by arrows. (b) Similar plot for the twoleg Hubbard ladder model, where t⊥ is replaced by tr .
From Ref.[3].

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo(a)comes small, the overall tendency does not changemuch in the ladder (except that larger tr is needed
for the bonding band to be incipient because theband width renormalization is weak), but for thebilayer model, superconductivity is strongly suppressed and it is no longer optimized in the incipient band situation.

(b)FIG. 3. mVMC calculation results of the pairing correlation for the two models with t′ = 0. (a) t⊥ /t dependence of the pairing correlation P s± (interlayer nearest
neighbor pairing)for the bilayer model with various values of U/t. The band filling is n = 1.06. (b) tr dependence of the pairing correlation P d (interchain nearest
neighbor pairing) for the two-leg ladder model with various values of U/t. The band filling is n = 1.03. From
Ref.[6].

bilayer, where λ is signiﬁcantly reduced.
We also performed multivariable variational
Monte Carlo calculation (mVMC) [4, 5] study forthe two models, and calculated the pairing correlation function against t⊥ (bilayer) or tr (two-leg
ladder)[6]. The results shown in Fig.3 resemble the
FLEX results of the eigenvalue of the Eliashberg
equation shown in Fig.2 (note that FLEX is performed only in the t⊥ ≥ t or tr ≥ t regime).
The above FLEX and mVMC results suggest
that, when U is suﬃciently large, superconductivity is strongly enhanced in both the bilayer andtwo-leg ladder models when the bonding band isincipient. As t⊥ or tr becomes small and the bonding band intersects the Fermi level, superconductivity is signiﬁcantly degraded for the bilayer model,
but not so strongly in the ladder. When U be-205

It is known that the spin ﬂuctuations within acertain moderate energy range work eﬀectively as
pairing glue, while those in the very low energy
range are pair breaking[7]. This, combined with thediﬀerence in the density of states (DOS) for the twomodels due to the diﬀerence in the dimensionality,
enables us to understand the obtained results[3, 6].
Namely, when the bonding band is incipient, thelow-energy spin ﬂuctuations is suppressed because
only the states around the bonding band edge cancontribute to their development, while there is signiﬁcant contribution from the bonding band to themoderate-energy spin-ﬂuctuations. This is the reason why superconductivity is strongly enhanced in
the incipient band situation. As t⊥ or tr is reduced
and the bonding band intersects the Fermi level,
the low energy spin ﬂuctuation rapidly develops in
the bilayer lattice because the DOS of the bonding
band at the Fermi level increases, while the development of the low energy spin ﬂuctuation in theladder is moderate because the DOS of the bonding band decreases as the Fermi level moves toward
the middle of the band. This explains the diﬀerence
between the two models. As for the eﬀect of U ,the band width is reduced by the renormalization,
so that the bonding band becomes incipient for asmaller t⊥ or tr for larger U . This eﬀect is essential for superconductivity to be strongly enhancedin the bilayer model because in this way the portionof the DOS that contributes to superconductivity
is squeezed into a narrower energy range. This effect is not so important in the ladder, where the
DOS at the band edge is diverging.

[1] K. Kuroki, T. Higashida, and R. Arita, Phys. Rev.
B 72, 212509 (2005).
[2] K. Matsumoto, D. Ogura, and K. Kuroki, Phys.
Rev. B 97, 014516 (2018).
[3] K. Matsumoto, D. Ogura, and K. Kuroki, J. Phys.
Soc. Jpn. 89, 044709 (2020).
[4] D. Tahara and M. Imada, J. Phys. Soc. Jpn. 77,
114701 (2008).
[5] T. Misawa et al., Comp. Phys. Comm. 235, 447(2019).
[6] D. Kato and K. Kuroki, Phys. Rev. Research, toappear (2020) (arXiv:1912.11983).
[7] A.J. Millis, S. Sachdev and C.M. Varma, Phys. Rev.
B 37, 4975 (1988).

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

Study of magnetism and transport phenomena in strongly
correlated quantum systems using computational approaches
including non-equilibrium methods
Norio KAWAKAMI
Department of Physics,
Kyoto University, Kitashirakawa, Sakyo-ku, Kyoto 606-8502
We

have

studiedstrongly

correlated

Brillouin zone at which the effective non-quantum systems with a focus on magnetism

Hermitian Hamiltonian cannot be diagonalized,

and transport phenomena, using computational

occur at the Kondo temperature, at which theapproaches, including non-equilibrium methods.

magnetic moments are screened, and the wellknown

We have analyzed the effect of correlations

crossoverbetweenlocalized

anditinerant f electrons takes place. [2]on the emergence of quantum oscillation in
narrow-gap topological Kondo insulators byperforming

dynamical

theoryHubbard model on the honeycomb lattice by(DMFT) calculations. We have demonstrated

DMFT as well as by numerically exact

that the interplay between correlations andquantum Monte Carlo simulations. We havenonlocal hybridizations can lead to observableshown that the DMFT approximation is very

quantum oscillations even without a Fermi

accurate in the Dirac semi-metallic phase,surface and that quantum oscillations in a three-where local moment formation is present, anddimensional system can be understood from thethe spin correlation length small. Moreover, we

physics of the two-dimensional planes in thehave shown that although the quantum critical

momentum space for which the hybridization in

point can only be understood by lattice

the direction of the magnetic field vanishes.

Quantum Monte Carlo techniques, DMFT givesThisa reasonable estimate of the critical interaction

scenario

shouldmean-field

Furthermore, we have studied the one-bandberelevantfor

understanding the observations in SmB6 andstrength. [3]

YbB12. [1]References[1] R. Peters, T. Yoshida, and N. Kawakami,

We have studied the impact of nonHermiticity due to strong correlations in f-

Phys. Rev. B 100, 085124 (2019).

electron materials, by analyzing the periodic-[2] Y. Michishita, T. Yoshida, and R. Peters,

Anderson model with local and nonlocal

Phys. Rev. B 101, 085122 (2020).

hybridization in the insulating and metallic
regimes by DMFT combined with numerical[3] M. Raczkowski, R. Peters, T. T. Phung, N.

renormalization group (NRG). We have shown

Takemori, F. F. Assaad, A. Honecker, and J.

that exceptional points, momenta in the

Vahedi, Phys. Rev. B 101, 125103 (2020).

206

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

Study on metal-insulator transition in nickel oxides
Yusuke NOMURA
RIKEN Center for Emergent Matter Science, 2-1 Hirosawa, Wako, 351-0198
Rare-earth nickelates LuNiO3 shows anenergy Hamiltonians for LuNiO3 eg manifold.

interesting phase diagram. Starting from theWe explicitly include the phonon degrees ofroom temperature, by cooling, it exhibits two-

freedom using newly-developed constrained

step transitions. The first one is the metal-

density functional perturbation theory (cDFPT)

insulator transition, and the second one is the[1]. We calculate the magnitude of phonon-

magnetic transition. At the first transition,

mediated electron-electron interaction as wellinterestingly, bond disproportionation occurs,

as the effective Coulomb interactions.

where one set of octahedra is compressed, andthe other is expanded. It means that thedisplacement of breathing phonon mode is
induced, and the Jahn-teller mode becomesinactive despite the (eg)1 configuration.

Fig. 2: Band structure of LuNiO3.

As a result, we find that coupling to the JahnTeller phonon mode indeed exists, which favors
Jahn-Teller distortion. However, at the sametime, non-local Coulomb interaction, whichfavors bond disproportionation, is substantial.

Fig. 1: Crystal structure of LuNiO3 at highTherefore, the coupling to Jahn-Teller modetemperature.

and the non-local Coulomb interaction competewith each other. In the real material, the latter
wins, and the bond disproportionation takes

In this study, we derive the effective low-207

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

place at the metal-insulator transition.

References[1] YN et al., Phys. Rev. Lett. 112, 027002(2014); YN and R. Arita, Phys. Rev. B 92,
245108 (2015)

208

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

Study on frustrated magnetism using machine learning
Yusuke NOMURA
RIKEN Center for Emergent Matter Science, 2-1 Hirosawa, Wako, 351-0198
Quantum spin liquid (QSL) is an exotic

the benchmark to check the usefulness ofstate of the matter, in which long-range

machine learning, and interesting applications

magnetic order is suppressed even at zero

to explore many-body physics have not yet

temperature. The QSL has attracted much

performed.

attention not only because of its uniqueIn the study, by using a newly-developed

property such as the fractionalization of themachine learning method (RBM+PP), which

excitation

combines the restricted Boltzmann machine

but

also

becauseofpossibleand(RBM) and pair-product (PP) state [2], we

memory devices. However, to show theobtain strong support for the existence of theexistence of the QSL in the realistic quantum-

QSL phase in the thermodynamic limit of thespin Hamiltonians, except for rare cases such as

spin-1/2 Heisenberg spin system on the square

the Kitaev model where the exact ground statelattice with the nearest-neighbor and next-

can be constructed by hand, we need highly

nearest-neighboraccurate numerical methods.

respectively, called the J1-J2 Heisenberg model.

applications

toquantum

computingTo achieve high accuracy, we focus on

The

QSL

exchanges,

phase

shows

J1

and

J2,

unconventional

machine learning. Machine learning is good at

excitation spectrum, in which both the singlet

extracting essential patterns in the data. By

and triplet excitations become gapless at high-applying

body

symmetry momenta (0,0), (p,0), (0,p), (p,p), in

accurate

support of the nodal Dirac dispersion of the

problems,

it

to

the

quantum-many

we

can

obtainanrepresentation of the quantum many-body wave

fractionalized spinon at (±p/2,±p/2).

functions (vector with exponentially large
dimensions) with a finite number of parameters(data compression). Starting from the seminal

References

work by Carleo and Troyer [1], now, there is[1] G. Carleo and M. Troyer, Science 355, 602

growing evidence for the powerfulness of the(2017).

machine learning method to study many-body[2] Y. Nomura, A. S. Darmawan, Y. Yamaji, andproblems. However, studies so far are limited to

M. Imada, Phys. Rev. B 96, 205152 (2017).

209

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

Variational Monte Carlo calculation of
two-dimensional Wigner crystal under latticecommensurability
Takeo KATO
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581
4

In this work, we studied the Wigner crystallization in spinless Fermion systems on asquare lattice at finite fillings by highly accurate simulation based on the many-variablevariational Monte Carlo (mVMC) method [5].
By performing the calculations up to the orderof 103 sites, we determine the transition points
for n = 2, 4, 8, 12, 16 and draw the phase diagram as functions of rs and 1/n as shown in
Fig. 1.
In contrast to the previous works by pathintegral renormalization group (PIRG) [4], wefound that the transition points at finite fillings
becomes significantly small. Nevertheless, thetransition points seem to continuously connect
to its continuum limits, which is consistent
with the previous works. We also examined thenature of the associated metal-insulator transition by calculation the momentum distribution

30

27

mVMC
HF
PIRG
DQMC

20
rs

Wigner crystallization, which breaks thetranslational symmetry due to the long-range
Coulomb interactions, is a fundamental problem in strongly correlated electron systems.
By the diﬀusive quantum Monte Carlo method(DQMC), the critical electron density for crystallization is obtained as rs = r0 /a∗B ≃ 150
for a three-dimensional electron gas (r0 =(πne )−1/2 , a∗B : the eﬀective Bohr radius) [1].
The critical value of rs is lowered down to rs ≃31 for a two-dimensional electron gas [2, 3].
In addition to the Coulomb interaction, thelattice commensurability reduces the critical
Coulomb interaction for crystallization [4].

28

1/8

10

1/4
0

0

0.1

0.2

0.3

0.4

0.5

n
Figure 5. Caption

Figure 1: Phase diagram in the plane of theelectron-density and rs . We also show the result of the Hartree-Fock (HF) approximation.
and the charge gap [6].

References[1] D. M. Ceperley and B. J. Alder, Phys. Rev.
Lett. 45, 566 (1980).
[2] M. Imada and M. Takahashi, J. Phys. Soc.
Jpn. 53, 3770 (1984).
[3] N. D. Drummond and R. J. Needs, Phys.
Rev. Lett. 102, 126402 (2009).
[4] Y. Noda and M. Imada, Phys. Rev. Lett.
89, 176803 (2002).
[5] T. Misawa, et al., Comput. Phys. Commun.
235, 447 (2019).
[6] M. K. Ritter, K. Ido, T. Misawa, and T.
Kato, in preparation.

210

put. Phys. Commun. 235, 447 (2019).
P. P. Ewald, Ann. Phys. 369, 253 (1921).
A. Grzybowski, E. Gwóźdź, and A. Bródka, Phys. Rev. B
61, 6706 (2000).

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

Magnetic and transport properties indouble-exchange systems with quasiperiodicity
Akihisa Koga
Department of Physics, Tokyo Institute of Technology, Meguro, Tokyo 152-8551

Quasicrystals (QCs) have no translational
symmetry but possess long-range order andspecial rotational symmetry. Since its first
discovery in 80’s, intensive efforts have been
made to search various type of QC materials
and to reveal characteristic features. In thisdecade, there have been several important experimental findings in QC materials, which include quantum critical behaviors in Au-Al-Yb
alloy and the emergence of superconductivity
in Al- Zn-Mg QC alloy. These experimentssuggest important roles of electronic correlations in these materials, and stimulate theoretical studies of the correlation and inducedordered phases in QC systems.
Motivated by this situation, we have studied possibility of the excitonic insulating (EI)
state in the two-band Hubbard model [1] andthe magnetism of the double exchange model[2] on the Penrose lattice. In the former
study [1], we considered the vertex model andthe center model on the Penrose lattice to determine their ground state phase diagram interms of the real-space mean-field theory. We
have revealed the existence of the stable EI
phase induced by small interband interactions[Fig. 1(a)]. We pointed out that this originates from the electron-hole paring for thecompletely or nearly degenerated states, which
is characteristic of the QC systems. In addition, we also discussed the spatial structure of
the order in terms of the so-called perpendicular space [Fig. 1(b)]. In the second study [2],we investigated the double exchange model on
the Penrose lattice in terms of a Monte Carlo
method. We showed the existence of the ferromagnetism and revealed that its general fea-

211

Figure 1: (a) EI order parameter of the twoband Hubbard model of the center model on
the Penrose lattice. Red dotted circles indicates the stable EI states with small interactions. (b) Analysis in the perpendicular space
for the stable EI phases.
tures such as temperature dependence resemble those of the system on the normal squarelattice. By means of the analysis in the parpendicular space, we discussed the characteristic
real-space structure of the magnetism in thisQC system.

References[1] K. Inayoshi, Y. Murakami, and A. Koga:
arXiv:2002.05870 (2020).
[2] Y. Takeuchi, Master thesis, Tokyo Institute of Technology (2020. Mar).

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

Numerical analysis for nonequilibrium dynamics incorrelated electron systems with orbital degrees of
freedomAkihisa Koga
Department of Physics, Tokyo Institute of Technology, Meguro, Tokyo 152-8551
LETTERS

J. Phys. Soc. Jpn.

2

t/JK1

Strongly correlated systems with orbital de!"#(a)grees of freedom can lead to intriguing spin
⇥10 3 2(a)R
150
systems. One of the important examples is1
the Kitaev quantum spin model, which consists
100
L
B
R
0
of anisotropic exchange couplings on the honlocalized
50
1
h (t)
eycomb lattice and shows the quantum spinSz hi�nerant
02liquid (QSL) state at low temperatures. Al⇥10 3 2
12 13
16 17
20
(b)!$# (b)though the Kitaev model was originally proL
M
R
14 15
18 19
150
11
6 7
2 3
posed theoretically,itisconsidered
to
be
re1
10
FIG. 1. Kitaev model on a honeycomb lattice with zigzag edges.
4 5
8 9
100
1
Red, blue, and,green
linesotherrepresentstrongly
x, y, and z bonds,
respectively.
alized in α−RuClcorre3 and0
Solid (open) circles represent spin-1/2 in the A (B) sublattice. In
50(A along= Na,Li) arelated systems
such
A2 IrO
this figure,fourasz (green)bonds
the a directionin the R
3 exist1
region, namely, L = 4.
Fig. 3. (a) Schematic picture of the Kitaev system on the honeycomb
latticeSzalso candidates of theR Kitaev QSL. So far, thewith two0 edges. The lattice is composed of three regions. The static magnetic
2
region, where
field h is0applied in the
10right (R) 20
30the magnetization
40 appears.
static properties in the Kitaev QSL have been
In the bulk (B) region, no magnetic field is applied and the QSL state is reduces the phase transition to the spin-polarized state aroundx
alized without the magnetization. Time-dependent pulsed magnetic field isthe main focus
thischarac-we re0.042[30]field
within and
the MFmany
theory. Therefore,
hc /JK ⇠ofintroduced in the left (L) region. (b) 20-site cluster used in the exact diagoThe numbers
represent
the Kitaev
index ofspin
the latticesite.with L = 50
strict ourselves
the casehR < hc to discuss the spin
FIG.nalization.
2. Real-time
dynamicsof thesystemteristic features
suchto as
thewithfractionalizationatransport inherent in the Kitaev QSL.
and Lb = 300 after the magnetic-field pulse at t = 0. The contourzof specific heats
andthethermalHall
effectst) onSchematicthe plane of the timet and the space
x in theplot of 1:S (x,We study the time evolution of the system upon stimuli ofFigure(a)picturesof oursetup
system (a) without and (b) with the M (white) region where the QSL
magnetic pulsein the Laregion(see Fig.
1) [44, 45]. To
have been thereported.
Thus,
simplequestionprocess
involved
realizedwithnoisspinpolarization
betweenL (left gray) and spin
region,which
explicitly
given inas thetheanomalous
this end, we introduce the time-dependent Majorana MF the-andis the
R (right gray) regions (see the top of the
The dashed lines" panels).
#
naturally arises
characteristic
2 result of the spinory. Thewhat
details ofisthethe
formulations
are shown in dyRef. [46].transport.
(b)themean-fieldAtrepresent x = vt withvelocity v (seeexp, text).
(9)hL (t)the= Majoranap
The
Hamiltonian
Eq.
(1)
is
obtainedasafermionmodel
by2
namical properties of the Kitaev QSL.
2
2⇡ 2 model.
dynamics of the Kitaev
L

t/JK1

R

R

applying the Jordan-Wigner transformation to the spin opera-

where A andare magnitude and width of the Gaussian
tors [47–49].
Furthermore,
introducing thetwo kindsof MajoMotivatedby this,we bystudied
spin
z
z
B
2hSpulse.
2hS
i(x + xdwe
/2, fixt), ⇥= hi¯=rA 2/J.
In=thefollowing,
A 5=(x,1 t)
andr i, ⇥6 (x, t) =A i(x, t)Fig. 4. Real-time evolution of the change in the local magnetization in therana fermions, and ¯ , from a complex fermion at each site,
A B
¯
i,
where
x
is
the
horizontal
coordinateofr and= system with h /J = 0.01 after the introduction of the pulsed magnetic fieldIn
thiscalculation,
weexaminethe
20-sitecluster
withxdtwohitransport through
the
QSL
by
investigatingrrp
Eq. (1) is rewritten as
the1/(2
other
edge
under
the
staticwith A = 1 and = 2/J shown as the dashed line (o↵set for clarity). Dottededges,
periodic
condition3). where
This theMF
theory boundaryis exact
when magnetic
hL is=imposed
hR = for0.field,
the real-space dynamics
X Kithe direction
perpendicularto MF
the edge,as shown
in Fig.as3(b).
JK X ⇣ A of B theA S B ⌘= JK1/2
we believethat
ourresultsarereliable
far lines represent the results for the pulse with A = 1 and = 2/J.
A B A BtheTherefore,hybridization
betweenthelocal
Majorana
H(t) =ii r a+b r + i r+b r
r i¯ r ¯ rIn thecluster,theareR and
L regions
includesites.theThereexas thesmall
fieldsapplied.
In theMFfourtheory,many4 r r[2]taev system with4 ar zigzag
edge structure
fermions
yieldsfinitespin
momentum.
ist wave
twelvesites inthetheB region,where
no magnetic
fieldisoneap- We
bodyfunctionisexpressed
asadirectproductofXX⇣⌘
hL (t) afterhR
A[Fig.1(a)]. We followed
plied.
Althoughcluster we treat issmall, thebyspintransi A ¯ thei rBdynamics¯ rB +i B ¯ B , the (2)alsobody
states,whose
is toodescribed
theto
MFshowstudiedthethetime-evolutionS= 1 Kitaev
model
2 r2R r r
2 r2L r r
port characteristicof the Kitaev
systemexpected
be capHamiltonian.
Determining
the MFsat is
eachtime tofromthe polarization even in the S = 1 Kitaev model. Namely, we exmagnetic field pulse applied to one of the edges
pect that in the S = 1 Kitaev model without a static magnetic
tured since
thereofexist
fourweplaquettes
withthetime-evolution
local Z2 sym- spinthemany-bodyexistencethesimilaranomalous
wave
function,computethewhererindicates
theposition
of
the
z
bond
(the
centerofthemetry,
which iscrucialforthetheextended
peculiar spintransport,
in[50–
the field, the spin degree of freedom is fractionalized into two;with the time-dependentA Majorana
mean-fieldoftheone-body
stateswithEuler
methodA
B
B[3].
z bond; see Fig. 1). r and ¯ r ( r and ¯ r ) are the Majoranatransport
B region.
In fact, we have confirmed that, in the initial state itinerant and localized quasiparticles owing to the existence55]. The time-dependent
field is explicitly given as,
theory andfermion
theoperators
exactconnected
diagonalization.
i magnetic
with the z bond in theWesublattice(t ! A1), the hsite-dependent
magnetization mzi appears only of the local Z2 symmetry. The pulsed magnetic field in the L
t
p
h
(t)
=
exp
,
where
A
and
arez constantsfor the region creates the itinerant and localized quasiparticle exciz
z
z
L
=
h
(t)
=
0,
[H,
⌘
]
=
0
A
(B),
as
shown
in
Fig.
1.
When
h
2
R
L
in the R2⇡ region (m9 = m19 = 0.217 and m10 = m20 = 0.220).
showed that, although no spin
moments
are in-r
Gaussianfollowing,
weoffixthethe
systemsizespin
as tations, while only the formers propagate in the whole sysat each r, and ⌘r (= i¯ rA ¯ rB ) is the Z2 local conserved quanFigure pulse.
4 showsInthethe
timedependencechangein thezz
duced in the
QSL
region,the asspin
mo- References
La moment= 50 and m
Lbz [== m
300,
themstatic
field asthehR pulsed= 0.01J
tem. In the R region, the itinerant and localized quasiparticles
K , andtity. Kitaev
In the case,the model
is solvable
the Hamiltonian
magneticii (t)i ( 1)] after1
pulse
parameters
as A The= 1 magnetic
and = 2J
field
is introduced.
moments
at the sites 1 and are hybridized by the static magnetic field due to the lack ofbilinear ininterms, and its low-energy
dispersion isK .
ments are isexcitedtheofopposite
edge region
2 in thediscussing
L region areby thethroughpulsed magnetic
field. the local Z2 symmetry, leading to the finite spin oscillations.
given asp"k ' v|k kK | around the K point with the velocBefore
the induced
spin transport
the M region,[1]weA.
Kitaev:
Ann.
Phys.
321
(2006)[Fig.1(b)]. ityWethat this
On
the
other
nowithoutmagnetic
momentinthe
the static
B 2.
re- The scenario for the spin transport is schematically shown in
3JK /4. This indicatesthat anomalousand ¯ are regarded asv = revealed
examinethe hand,systemthisregion,appears
namely,gion, which
thesites
existence
local withcon- Fig. 3(a).
the itinerant
and localized
Majorana
respectively.
magnetic
field hisR consistentis appliedwith
to thein theofRthe
regionspin transport
originates
from
thefermions,fractionalserved
quantity.
Oncase,
theY.
other
hand,
after
some
delay,
thequanspin andThe spin fractionalization in the S = 1 Kitaev model has
Murakami,
A. Koga= LMinakawa,
1. In thisthere
are nolocal
conserved
Since the magnetic field hybridizes two kinds of the [2]LRT.
asuch as aoscillation
the sitesand 10 in the
R region.
ization of Majorana
the spinsin thetheHamiltonian
Kitaev isQSL
intotities,leading istoinduced
nonzeroat local
spin9 moments.
Figure2(a) been suggested in the thermodynamic properties
fermions,
no longer
exactlyJ.
Nasu:
arXiv:1912.10599(2019).
double-peak structure in the specific heat.21) We note that
This
means
that
the
wave-packet
triggeredbythemagnetic
shows the contour plot of the changes in the spin moments
solvable.
Here, we apply
the Hartree-Focktype decouthe itinerant
and localized
Majorana
fermions.
the L region
the t)iR regionthrough
re- the higher temperature peak is closely related to the nearest(x, t),inwhereO(x, t)reaches= hO(x,hO(x,1)i.theAsB exS zpulse
pling to the interaction on the z bond as i rA rB i¯ rA ¯ rB ⇠gion without spin oscillations. Since mix (t) = myi (t) = 0 for neighbor spin-spin correlations Ci j . It is known that the higher
A
B
A
B
A
A
Namely, after
excitation, only the itiner[3]pected,
A.
Koga,
T.
Minakawa,
Y.
Murakamii r r ⇥thewe find thatthewavepacket createdbythe magnetic- and
1 (x, t)+⇥2 (x, t)i¯ r ¯ r ⇥1 (x, t)⇥2 (x, t) i r ¯ r ⇥3 (x, t)
any sites, we can say that the spin moment plays no role in temperature peak in the S = 1/2 case corresponds to the mot)i rB ¯ rB +⇥3 (x,cant)⇥4 (x,propagatet) i rA ¯ rB ⇥5 (x, t)through⇥6 (x, t)i¯ rA rB +⇥4 (x, fermionsfield pulse flows to the right edge. Note that its velocity al- tion of the itinerant Majorana fermions. Therefore, one canant Majorana
the spintransportJ.
in theS = 1 Kitaev
model.
This remarkable
J.
Nasu:
Phys.
Soc.
Jpn.
89 (2020)⇥5 (x, t)⇥6 (x, t), where we have introduced the six kinds ofmost coincides with the Majorana velocity of the genuine Kiphenomenon is similar to that in the S = 1/2 Kitaev model,30) expect that the flow of Ci j is regarded as the motion of the itinA absent
B
the QSL, hencethemagnetization
is033701.
x- and t-dependent MFs as ⇥1 (x, t) = hi¯ r ¯ r i ⌘ h⌘i(x, t),taev
model,
whichis shown
as the dashedline in Fig.
2(a). erant quasiparticles in the S = 1 Kitaev model. Figures 5(a),wherethev,spintransport
is mediatedby the itinerant
Majorana
Thisfermions.
indicatesTherefore,that the propagation
is attributed
to the gapless= itinerant
hi rA rB i ⌘ h⇠i(x,t), ⇥3 (x, t) reaches= hi rB ¯ rB i =⇥2 (x, t)there. Whenthe
Majorana
ourresultssuggesttheexistenceof
the 5(b) and 5(c) show the real-time evolution of the change in the
z
z
A A
R

2

2

2hS B i(x, t) =2hS i(xxd /2, t), ⇥4 (x, t) = hi r ¯ r i =Majorana
in thewhichbulk withina small static
itinerantexcitation
quasiparticles,
is not associatedwithfield.
the spinnearest-neighbor spin correlations on the x-, y-, and z-bonds
except for the L region. Oscillations appear in the spin-spin

3

212

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

Novel topological phenomena induced by strongly
correlated dynamics
Tsuneya YOSHIDA
Department of Physics, University of Tsukuba
Tsukuba, Ibaraki 305-8571

Focusing on non-Hermitian topological phenomena, we have analyzed correlated systems.
Speciﬁcally, we have analyzed exceptionalband touching for quasi-particle spectrum of
correlated systems in equilibrium [1, 2, 3].
Our works [1, 3] have discovered novel exceptional band touching which we dub symmetryprotected exceptional rings (surfaces) in two(three) dimensions, respectively. In addition,
we have also elucidated that for Kondo insulators, such non-Hermitian band touching
emerges around Kondo temperatures. Furthermore, focusing on dynamics of open quantum systems with two-body loss, we have discovered a non-Hermitian fractional quantum
Hall state (see Fig. 1), which extends thenon-Hermitian physics to long range entangledstates. The above set up is relevant to the cold
atoms with two-body loss.

ImE

(a)

0

(b) 0

-0.05

-0.05

-10.1

-10.05
ReE

-10

nθ=0
nθ=1
nθ=2
nθ=3
nθ=4
nθ=5

ticular, we have found that mirror symmetry
and non-Hermitian topology induces novel skineﬀect protected by mirror symmetry which
may be realized for an electric circuit composed of negative impedance converters. Besides that we have elucidated the other topological phases protected by crystalline symmetry [7, 8, 9]. For instance, we have extended
the notion of higher order topology to correlated systems, which discovered a new topological state, a higher-order topological Mott
insulating state [9]. We also have proposed aneﬃcient method to compute Z4 topological invariant for glide symmetric systems [7].

References[1] Kazuhiro Kimura, Tsuneya Yoshida, andNorio Kawakami: Phys. Rev. B 100,
115124 (2019).
[2] Yoshihiro Michishita, Tsuneya Yoshida,
and Robert Peters: Phys. Rev. B 101,
085122.

-0.85 -0.8 -0.75
100+10× ReE

Figure 1: (Color Online). Energy spectrum of
non-Hermitian Hamiltonian. Panel (b) is themagniﬁed version of the range in −10.085 <ReE < −10.075.
In addition to be above works on nonHermitian systems, we have addressed topology with crystalline eﬀects [5, 6]. In par-213

[3] Tsuneya Yoshida, Robert Peters, Norio Kawakami,
Yasuhiro Hatsugai:
arXiv:2002.11265.
[4] Tsuneya Yoshida, Koji Kudo, and Yasuhiro Hatsugai: Sci. Rep. 9, 16895(2019).
[5] Tsuneya Yoshida and Yasuhiro Hatsugai:
Phys. Rev. B 100, 054109 (2019).

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[6] Tsuneya
Yoshida,Tomonari
Mizoguchi,
andYasuhiro
Hatsugai:
arXiv:1912.12022.
[7] Tsuneya Yoshida, Akito Daido, Norio
Kawakami, and Youichi Yanase: Phys.
Rev. B 99, 235105 (2019).
[8] Kazuhiro Kimura, Tsuneya Yoshida, andNorio Kawakami: JPS Conf. Proc. 30,
011012 (2020).
[9] Koji Kudo, Tsuneya Yoshida, and Yasuhiro Hatsugai: Phys. Rev. Lett. 123,
196402 (2019).

214

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

Two-Channel Kondo Eﬀect Emerging in a
Multiorbital Anderson Model Hybridized with Γ7 andΓ8 Conduction Bands
Takashi HOTTA
Department of Physics, Tokyo Metropolitan University
1-1 Minami-Osawa, Hachioji, Tokyo 192-0397

To explore the emergence of a quantum critical point (QCP) between non-Fermi-liquid andFermi-liquid phases, we investigate the statesin the vicinity of two-channel Kondo phase.
For the purpose, we analyze a seven-orbital impurity Anderson model hybridized with three
conduction bands (one Γ7 and two Γ8 ) byemploying a numerical renormalization group(NRG) method.
First we consider one f -electron state, whichis the eigenstate of spin-orbit and crystalline
electric field (CEF) potential terms. Under
the cubic CEF potential, we obtain Γ7 doublet
and Γ8 quartet from j = 5/2 sextet, whereas
we find Γ6 doublet, Γ7 doublet, and Γ8 quartet from j = 7/2 octet. By using those oneelectron state as bases, we describe the local
f -electron Hamiltonian as
Hloc =∑†(λj + Bj,µ + Ef )fjµτ
fjµτ

j,µ,τ

+

∑ ∑ ∑

j1∼j4 µ1∼µ4 τ1∼τ4

Iµj11jτ21,jµ32jτ42 ,µ3 τ3 µ4 τ4 (1)× fj†1 µ1 τ1fj†2 µ2 τ2fj3 µ3 τ3fj4 µ4 τ4 ,where fjµτ denotes the annihilation operator ofa localized f electron in the bases of (j, µ, τ ),j is the total angular momentum, j = 5/2 and7/2 are denoted by “a” and “b”, respectively,µ distinguishes the cubic irreducible representation, Γ8 states are distinguished by µ = α
and β, while Γ7 and Γ6 states are labeled byµ = γ and δ, respectively, τ is the pseudo-spin
which distinguishes the degeneracy concerning

the time-reversal symmetry, and Ef is the f electron level to control the local f -electron
number at an impurity site.
As for the spin-orbit coupling term, we obtain λa = −2λ and λb = (3/2)λ, where λis the spin-orbit coupling of f electron. For
Pr ion, we set λ = 0.1eV in this research.
Concerning the CEF potential term for j =5/2, we obtain Ba,α = Ba,β = 1320B40 /7
and Ba,γ = −2640B40 /7, where B40 denotes
the fourth-order CEF parameter for the angular momentum ℓ = 3. Note that the sixthorder CEF potential term B60 does not appear
for j = 5/2, since the maximum size of thechange of the total angular momentum is lessthan six. On the other hand, for j = 7/2,
we obtain Bb,α = Bb,β = 360B40 /7 + 2880B60 ,Bb,γ = −3240B40 /7 − 2160B60 , and Bb,δ =360B40 − 3600B60 /7. Note also that B60 term
appears in this case. In the present calculations, we treat B40 and B60 as parameters.
The matrix element of the Coulomb interaction is expressed by I. Here we do notshow the list for the explicit forms of I, butthey are expressed with the use of four SlaterCondon parameters, F 0 , F 2 , F 4 , and F 6 . Although the Slater-Condon parameters of a material should be determined from experimental results, here we simply set the ratio as
F 0 /10 = F 2 /5 = F 4 /3 = F 6 = U , where U is
the Hund’s rule interaction among f orbitals.
it is reasonable to set U as 1 eV in this research.
In Fig. 1, we show the local CEF ground-215

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

0.01-0.005

G5

-0.01
-0.1

-0.05

B60 [meV]

0.0

U=1.0
l=0.1
W=0.001

G3

x=0.0

Entropy

0.005

x=1.0
G1

0.0
B40 [meV]

0.05

0.1

Figure 1: Local CEF ground-state phase diagram on the plane of (B40 , B60 ) for n = 2.
The line is given by B40 = W x/15 and B60 =W (1 − |x|)/180 with W = 10−3 .
state phase diagram for n = 2, obtained fromthe diagonalization of Hloc . For n = 2, theground-state multiplet for B40 = B60 = 0
is characterized by total angular momentum
J = 4. Under the cubic CEF potentials, thenonet of J = 4 is split into four groups as Γ1singlet, Γ3 doublet, Γ4 triplet, and Γ5 triplet.
Among them, Γ4 triplet does not appear as asolo ground state under the cubic CEF potential with Oh symmetry.
Now we include the Γ7 and Γ8 conduction
electron bands. Here we consider only the hybridization between conduction and j = 5/2
electrons. The Hamiltonian is given byH =∑εkc†kµτckµτ +k,µ,τ

+ Hloc ,∑Vµ(c†kµτ faγτ +h.c.)

k,µ,τ

(2)

where εk is the dispersion of conduction electron with wave vector k, ckγτ is an annihilation operator of conduction electrons, and Vµdenotes the hybridization between f electronin the µ orbital and conduction electron of theµ band. Here we set Vα = Vβ = Vγ = V .
In this research, we analyze this model by
employing the NRG method [1]. We introduce
a cut-oﬀ Λ for the logarithmic discretizationof the conduction band. Due to the limitationof computer resources, we keep M low-energy

2
x=0.0
U=1.0
1.8
x=0.18715
1.6
l=0.1
x=1.0
1.4
W=0.001
1.2
V=0.6
L=8.0
1
M=5000
0.8
0.6
log f
0.4
0.5log 2
0.2
0 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0
10 10 10 10 10 10 10 10 10 10 10
Temperature T [eV]

Figure 2: Entropies vs. temperature for x = 0,0.18715, and 1 with W = 10−3 .
states. Here we use Λ = 8 and M = 5, 000.
Note that the temperature T is defined as T =DΛ−(N −1)/2 in the NRG calculation, where N
is the number of the renormalization step andD is a half of conduction band width. Here we
set D = 1 eV.
In Fig. 2, we show the calculation resultsfor entropies, when we change the local CEF
states from Γ3 doublet (x = 0) to Γ1 singlet(x = 1) on the line in Fig. 1. For x = 0,
we find an entropy of 0.5 log 2 at low temperatures, indicating the two-channel Kondo state(non-Fermi-liquid phase), while for x = 1,
we obtain the CEF singlet state (Fermi-liquid
phase). After repeating the NRG calculations,at x = 0.18715, we find an entropy plateau of√log ϕ with ϕ = ( 5 + 1)/2. This is the sameas the residual entropy for the three-channel
Kondo eﬀect. We consider that the entropy
log ϕ suggests the emergence of QCP between
non-Fermi-liquid and Fermi-liquid phases. Details will be discussed elsewhere in future [2].

References[1] H. R. Krishna-murthy, J. W. Wilkins, andK. G. Wilson, Phys. Rev. B 21, 1003
(1980).
[2] T. Hotta, in preparation.

216

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

Research of Two-Channel Kondo Eﬀect
in Transuranium Systems
Dai MATSUI and Takashi HOTTA
Department of Physics, Tokyo Metropolitan University
1-1 Minami-Osawa, Hachioji, Tokyo 192-0397

On the basis of the electron-hole relation
between f 2 and f 4 states in a j-j coupling
scheme, we hit upon an idea that two-channel
Kondo phenomena can be observed even in cubic transuranium compounds including Np and
Pu ions with local Γ3 non-Kramers doublet
ground state. To clarify this point, we analyze a seven-orbital impurity Anderson modelhybridized with degenerate Γ8 conduction electron bands by employing a numerical renormalization group (NRG) method [1].
The model Hamiltonian is given byH =∑εk c†kµτ ckµτ +

k,µ,τ

+

∑

∑

V (c†kµτ faµτ + h.c.)

k,µ,τ

†
fjµτ
(λj + Bj,µ + Ef )fjµτ

j,µ,τ

+

∑ ∑ ∑

j1∼j4 µ1∼µ4 τ1∼τ4

Iµj11jτ21,jµ32jτ42 ,µ3 τ3 µ4 τ4(1)

× fj†1 µ1 τ1fj†2 µ2 τ2fj3 µ3 τ3fj4 µ4 τ4 ,where εk denotes the energy of Γ8 conduction electron with wave vector k, ckµτ indicates the annihilation operator of Γ8 conduction electron, µ distinguishes the cubic irreducible representation, Γ8 states are distinguished by µ = α and β, while Γ7 and Γ6 states
are labeled by µ = γ and δ, respectively, τ isthe pseudo-spin which distinguishes the degeneracy concerning the time-reversal symmetry,
V is the hybridization between conduction andlocalized electrons, fjµτ indicates the annihilation operator of a localized f electron in thebases of (j, µ, τ ), j is the total angular momentum, j = 5/2 and 7/2 are denoted by “a”
and “b”, respectively, and Ef is the f -electronlevel to control the local f -electron number at
an impurity site. In the present case, we consider the hybridization between Γ8 conduction
electrons and the Γ8 quartet of j = 5/2.
As for the spin-orbit coupling term, we obtain λa = −2λ and λb = (3/2)λ, where λis the spin-orbit coupling of f electron. The
magnitude of λ depends on the kind of actinide atoms, but in this research, we set
λ = 0.3 eV, since we consider Np and Pu
ions. Concerning the CEF potential term for
j = 5/2, we obtain Ba,α = Ba,β = 1320B40 /7
and Ba,γ = −2640B40 /7, where B40 denotes
the fourth-order CEF parameter for the angular momentum ℓ = 3. Note that the sixthorder CEF potential term B60 does not appear
for j = 5/2, since the maximum size of thechange of the total angular momentum is lessthan six. On the other hand, for j = 7/2,
we obtain Bb,α = Bb,β = 360B40 /7 + 2880B60 ,Bb,γ = −3240B40 /7 − 2160B60 , and Bb,δ =360B40 − 3600B60 /7. Note also that B60 terms
appear in this case. In the present calculations,
we treat B40 and B60 as parameters.
The matrix element of the Coulomb interaction is expressed by I. To save space, we donot show the explicit forms of I here, but they
are expressed by four Slater-Condon parameters (F 0 , F 2 , F 4 , F 6 ) and Gaunt coeﬃcients.
As for the magnitudes of Slater-Condon parameters, first we set F 0 = 10 eV by hand.
Others are determined so as to reproduce excitation spectra of U4+ ion with two 5f electrons. The results are F 2 = 6.4 eV, F 4 = 5.6

217

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

Figure 1: Color contour map of the entropy
for n = 4 on the plane of (B40 , B60 ) for V =1.0
and T = 2.6 × 10−6 . White curves denote theboundaries among local CEF ground states.
eV, and F 6 = 4.1 eV [2].
In this research, we analyze the model by
employing the NRG method [3]. We introduce
a cut-oﬀ Λ for the logarithmic discretizationof the conduction band. Due to the limitationof computer resources, we keep M low-energy
states. Here we use Λ = 5 and M = 2, 500.
Note that the temperature T is defined as T =DΛ−(N −1)/2 in the NRG calculation, where N
is the number of the renormalization step andD is a half of conduction band width. We set
D = 1 eV in this research.
In Fig. 1, we show the contour color mapof the entropy on the plane of (B40 , B60 ) for
V = 1.0 and T = 2.6 × 10−6 . To visualize
precisely the behavior of entropy, we define
the color of the entropy between 0 and 1.5,
as shown in the right color bar. Note that thegreen and yellow regions indicate the entropyof log 2 and 0.5 log 2, respectively. We immediately notice that a region with an entropyof 0.5 log 2 (yellow region) almost corresponds
to that of the Γ3 CEF ground state. Notethat for large B60 , the color becomes green,
but when we decrease the temperature, wefind the entropy of 0.5 log 2 even for large B60 .
The results strongly suggest the emergence of
quadrupole two-channel Kondo eﬀect for thecase of n = 4. Note that in the present model,
we have observed the discontinuous change inentropy behavior between yellow (Γ3 ) and red(Γ1 or Γ5 ) regions. Quantum critical behavior which may occur between non-Fermi-liquid
and Fermi-liquid phases will be discussed elsewhere in future.
Here we remark that there exists a blue
region corresponding to log 3 in the local Γ5state. This is quite natural, since the local Γ5state is triply degenerate. The Γ5 moment isscreened by Γ8 conduction electrons and thus,
it is expected that the conventional Kondo effect occurs in the Γ5 region, although the magnitude of the Kondo temperature significantly
depends on the hybridization and excitationenergy.
Note also that in the region of B60 < 0,
we observe some blurry yellow spots along theboundary curve between Γ1 singlet and Γ5
triplet local ground states. The QCP is known
to appear between the CEF and Kondo singlet
states for n = 2 [4]. Thus, we deduce thatthose spots form a QCP curve, although we
could not obtain enough amounts of numericalresults to depict the smooth curve here.
Finally, we provide a brief comment on
actual materials to observe the two-channel
Kondo eﬀect for the case of n = 4. In rareearth ions, the case of n = 4 corresponds toPm3+ , but unfortunately, there exist no stable
isotopes for Pm. Thus, we turn our attention
to actinide ions such as Np3+ and Pu4+ with
5f 4 configurations. It may be diﬃcult to synthesize new Np and Pu compounds, but we
expect that Np 1-2-20 compound will be synthesized in future.

References[1] D. Matsui and T. Hotta, JPS Conf. Proc.
30, 011125 (2020).
[2] T. Hotta, Phys. Rev. B 80, 024408 (2009).
[3] H. R. Krishna-murthy, J. W. Wilkins, andK. G. Wilson, Phys. Rev. B 21, 1003
(1980).
[4] T. Hotta, Physica B 536, 203 (2018).

218

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

First-principles and quantum many-body calculations
for electronic states and superconductivity in themulti-band systems
Yoshiaki Ōno
Department of Physics, Niigata University
Ikarashi, Nishi-ku, Niigata, 950-2181

Tungsten bronzes Ax WO3 (A=Na, K, Rb,
Cs) are known to show the superconductivity for 0.15 < x < 0.35, where the superconducting transition temperature Tc is found toincrease rapidly with decreasing x up to 7K.
In addition, at further low doping x ∼ 0.05,
remarkable high-Tc superconductivities are reported at Tc =91K for A=Na and Tc =120K for
A=H. The systems also show rich structuralphases such as cubic for x > 0.4, tetragonal
for 0.1 < x < 0.4, orthorhombic for 0.03 <x < 0.1 and monoclinic for x < 0.03, which are
well accounted for by the first-principles calculations. However, the superconductivity inwhich the structural phases and those transitions are considered to play important roles fordetermining Tc was not discussed there. The
purpose of this work is to study the superconductivity in Ax WO3 on the basis of the firstprinciples calculation and to discuss the eﬀectof the ferro-orbital fluctuation which is foundto be largely enhanced around the structuralphase transition and is expected to induce thepairing interaction as previously discussed in
the iron-based superconductivity.
We calculate Tc on the basis of thefirst-principles calculation by using Quantum
ESPRESSO and find that the x dependenceof Tc is qualitatively consistent with the experiments. The significant x dependence of Tc
is mainly due to that of the electron-phonon
coupling constant λ which decreases with increasing x because of the screening eﬀect as

discussed by several authors. However, thecalculated value of Tc is considerably smaller
than the experimental values, and hence thepairing interaction mediated by the electronic
fluctuation such as the orbital fluctuation is expected to be responsible for the pairing mechanism in addition to the BCS phonon mechanism which is exclusively taken into account in
the present first-principles calculation. Then,
we also calculate the pairing interaction on thebasis of the eﬀective three band model consistsof W-t2g orbitals by taking into account of the
Jahn-Teller electron-phonon coupling within
the random phase approximation and find thatthe pairing interaction is largely enhanced dueto the ferro-orbital fluctuation which becomeslarge around the structural phase transition[1].
Recently, the Pr 1-2-20 systems PrT2 X20
have been intensively investigated concerning
their specific features including quadrupoleorders, superconductivity and non-Fermi liquid behaviors. The crystalline electric field
ground states of the Pr ions are the nonKramers doublets Γ23 which have the electric
quadrupole moments with no magnetic moment. In fact, PrTi2 Al20 and PrV2 Al20 exhibit
the ferro-quadrupole (FQ) and the antiferroquadrupole (AFQ) orders, respectively. They
also show remarkable superconductivities coexisting with the quadrupole orders. To discuss the quadrupole order and the superconductivity, the details of the energy bands andthe Fermi surfaces are important. In the Pr

219

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

1-2-20 systems, the de Haas-van Alphen experiments revealed that the Fermi surfaces
are well accounted for by the first-principles
band calculations for the corresponding La 12-20 systems, indicating that the Pr 4f electrons are suﬃciently localized. Therefore, thequadrupole orders are considered due to the
Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction between the quadrupole moments of
the localized Pr 4f electrons. The purpose ofthis study is to evaluate the RKKY interaction
on the basis of the realistic energy band structure extracted from the first-principles calculation and to discuss the quadrupole orders in
PrTi2 Al20 and PrV2 Al20 .
We perform the first-principles calculations
for LaT2 Al20 (T =Ti, V) by using WIEN2k
code where we also employ the GGA+U
method with U = 60eV to exclude the La4f components near the Fermi level. The obtained Fermi surfaces of LaTi2 Al20 are in good
agreement with those in the previous study.
Then, we construct the tight-binding modelfor the conduction electrons so as to reproduce the electronic structure near the Fermi
level by using the maximally localized Wannierfunctions which consist of 196 orbitals: La-d(5 orbitals × 2 sites), La-s (1 orbital times 2
sites), T-d (5 orbitals × 4 sites), T-s (1 orbital × 4 sites), Al-p (3 orbitals × 40 sites
and Al-s (1 orbital × 40 sites) in the conventional unit cell. Based on the eﬀective model,
we calculate the quadrupole susceptibility of
the conduction electrons and then evaluate the
RKKY interaction mediated by the conductions. The obtained results indicate that thewave vectors of the expected quadrupole orders are Q = (0, 0, 0) (FQ) for PrTi2 Al20 and
Q = (π/a, π/a, 0) (AFQ) for PrV2 Al20 , respectively, as consistent with the experiments[2].
Quasicrystals first discovered by Shechtman
et al. in Al-Mg alloy have been extensively
investigated as they show various anomalous properties due to their specific structures
which are not periodic but have kinds of trans-

lational orders with the 5-fold (8-, 10- and 12fold) rotational symmetry. Recently, the bulk
superconductivity at a transition temperature
Tc ∼ 0.05K has been discovered by Kamiya
et al. in the Al-Mg-Zn quasicrystal and its
approximant, where the basic cluster is considered to be the Bergman-type cluster, andhas attracted much attention for its symmetry
and pairing mechanism.
As for the theoretical study, the Al30 Mg40 Zn30 1/1 approximant was previously studied on the basis of
the first-principles calculation and the originof the pseudo gap near the Fermi level was discussed. However, the electronic state with thecomposition ratio corresponding to the superconducting compound was not discussed there.
Then, we perform the first-principles calculations (OpenMX code) for the Al-Mg-Zn 1/1
approximants with various composition ratios
including the superconducting one and findthat the over-all structure of the density ofstates seems to be almost independent of thecomposition ratio but the energy band dispersions near the Fermi level and the Fermi
surfaces are largely dependent on the composition ratio. On the basis of the obtained
electronic states, we also discuss the possiblesuperconductivity[3].

References[1] Takuya Sekikawa, Rai Watabe, Jun
Ishizuka, Yoshihiro Nitta, Kazuhiro Sanoand Yoshiaki Ōno, JPS Conf. Proc. 30,
011043 (2020).
[2] Yuto Iizuka, Takemi Yamada, Katsurou
Hanzawa and Yoshiaki Ōno, JPS Conf.
Proc. 30, 011152 (2020).
[3] Masaki Saito, Takuya Sekikawa andYoshiaki Ōno, 14th International Conference on the Structure of Non-Crystalline
Materials (NCM14).

220

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

Search for new types of skyrmion crystals based on
point group symmetry
Satoru Hayami
Department of Applied Physics, The University of Tokyo, Tokyo 113-8656

Multiple-Q magnetic states, which are characterized as a superposition of diﬀerent wave
vectors, open a new path toward applications
to next-generation electronics and spintronics
device in condensed matter physics. Amongst
them, the superposition of spiral states gives
rise to unconventional noncollinear and noncoplanar magnetic textures, such as the magnetic vortex crystal and magnetic skyrmion
crystal, which become a source of unconventional topological excitations and transport
properties [1]. Especially, since the discoveryof the magnetic skyrmions in noncentrosymmetric magnets, their studies have been extensively done in both theory and experiments.
Meanwhile, recent theoretical studies pointout the emergence of the skyrmion crystals incentrosymmetric Mott insulators. They arestabilized by taking into account various factors, such as frustrated exchange interactions
and thermal fluctuations [2, 3, 4, 5]. Moreover,another interesting situation is brought by considering the interplay between itinerant electrons and localized spins. There, eﬀective longranged interactions mediated through itinerant electrons play an important role in stabilizing the skyrmion and vortex crystals [6, 7, 8],
which might provide a deep understanding of
the microscopic mechanism of the multiple-Q
skyrmion, vortex, and hedgehog crystals foundin d- and f -electron materials.
In this project, we investigate the skyrmion
physic in itinerant magnets by considering
an anisotropic magnetic interaction [9]. We
start from the Kondo lattice model with221

the anisotropic xxz-type spin-charge coupling.
The Hamiltonian is represented by
H = −

∑

tij c†iσ cjσ + Jz

i,j,σ

+ Jxy

∑

i,σ,σ ′∑
i,σ,σ ′zz †σσσ
′ Si ciσ ciσ ′

y
y †
x
x
(σσσ
′ Si + σ
σσ ′ Si )ciσ ciσ ′ , (1)where c†iσ and ciσ are creation and annihilation operators of an itinerant electron at sitei and spin σ, respectively. The first term represents the hopping of itinerant electrons andthe second and third terms denote the spincharge coupling where |Si | = 1 and σ is thePauli matrix. Jxy and Jz are the anisotropic
coupling constants for xy and z spin components, respectively.
We derive an eﬀective spin model in Eq. (1)by assuming (i) the weak-coupling regime for
Jxy , Jz ≪ W (W is the bandwidth) and (ii)the multiple peak structures of the bare magnetic susceptibility. Then, the spin model isobtained as
H = −2∑

˜ Q + K̃Γ2Q ),
(JΓ
ν
ν

(2)

ν=1−3

√
x Sx
where ΓQν=
1 − A2 (SQ
+
ν −Qν
y
y
z
z
SQν S−Qν ) + ASQν S−Qν .
We adopt a
triangular lattice with six peaks in thebare susceptibility:±Q1 = ±(π/3, 0),√±Q2 = ±(−π/6, 3π/6), and ±Q3 =√±(−π/6, − 3π/6). The first term represents
the bilinear interaction with J˜ = 1 andthe second term represents the biquadraticinteraction with the positive coupling constant K̃ = K/N > 0 (N is the numberof spins). The former originates from theActivity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo
(a)

(b)

(c)

(d)

Figure 1: Snapshots of the spin configurations (see the text). The arrows represent the xy
component and the contour shows the z component of the spin moment.

second-order perturbation, while the latter
arises from the fourth-order perturbation.
The magnetic anisotropy is represented by√
A; 0 ≤ A < 1/ 2 for the easy-axis interac√tion and 1/ 2 < A ≤ 1 for the easy-planeinteraction.
The magnetic phase diagram of the eﬀective
spin model in Eq. (2) is investigated by performing simulated annealing from high temperature. We apply the Metropolis local updates under periodic boundary conditions inboth directions for the 962 spins. To find the
optimal spin configurations, we gradually reduce the temperature with the rate Tn+1 =αTn , where Tn is the temperature in the nth
step. We start from the initial temperature
T0 = 10−1 -100 . We typically set α ∼ 0.99999
and the final temperature T = 10−4 -10−3 .
As a result, we find several magnetic ordered states with a net scalar chirality including the skyrmion and vortex crystals by constructing the ground-state phase diagram ina wide range of model parameters. We showthat two types of the skyrmion crystal phaseswith the topological number of two are realized for large K: one is identified as the superposition of the magnetic vortex structuresin the xy spin component and the sinusoidal
wave in the z spin component [Fig. 1(a)] andthe other is denoted as the bubble-type magnetic configurations satisfying the sixfold rotational symmetry in both xy and z spin components [Fig. 1(b)]. We also find a triple-Q

222

state with a finite net scalar chirality, which
is diﬀerent from the skyrmion crystal by increasing the easy-plane anisotropic interaction[Fig. 1(c)]. Moreover, we find that the strong
easy-axis anisotropic interaction gives rise to
the triple-Q magnetic bubble crystal with thecollinear z spin component [Fig. 1(d)].

References[1] N. Nagaosa, Y. Tokura, Nat. Nanotech.
8 (12) (2013) 899–911.
[2] T. Okubo, S. Chung, H. Kawamura, Phys.
Rev. Lett. 108 (2012) 017206.
[3] A. O. Leonov, M. Mostovoy, Nat. Commun. 6 (2015) 8275.
[4] S.-Z. Lin, S. Hayami, Phys. Rev. B 93
(2016) 064430.
[5] S. Hayami, S.-Z. Lin, C. D. Batista, Phys.
Rev. B 93 (2016) 184413.
[6] Y. Akagi, M. Udagawa, Y. Motome, Phys.
Rev. Lett. 108 (2012) 096401.
[7] R. Ozawa, S. Hayami, K. Barros, G.W. Chern, Y. Motome, C. D. Batista, J.
Phys. Soc. Jpn. 85 (10) (2016) 103703.
[8] S. Hayami, R. Ozawa, Y. Motome, Phys.
Rev. B 95 (2017) 224424.
[9] S. Hayami, submitted.

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

Multipole and anisotropic superconductivity fromfirst-principles
Jun Ishizuka
Department of Physics, Graduate School of Science, Kyoto University,
Kyoto 606-8502, Japan

Recently discovered heavy fermion superconductor UTe2 [1], as a member of the ferromagnetic

and three two-dimensional Z2 numbers by Fermi
surface formulas [9, 10, 11]:superconductors [2, 3], provides a new platform toω=study the spin-triplet superconductivity. The spin-

1X
n(Ki )
2(mod 2).

(1)

Ki

triplet state in UTe2 is supported by various exper-

The GGA+U Fermi surfaces indicate topologicalimental results: large upper critical field Hc2 alongsuperconductivity for the reasonable values of U ,the three principal axes [1, 4], reentrant supercon-

and accompanied surface states are provided. We

ductivity under magnetic field along b axis [5], andalso predict the superconducting gap node at zero

the temperature dependence of the Knight shift inmagnetic field and under magnetic field along b

superconducting state [6].

axis. Based on these results, we proposed multiple
study

superconducting phases under the magnetic field.

[7].

Our another aim of this project is to study the

Theoretically, a band structure has been investi-

exotic phenomena due to the spatial inversion sym-

gated for UTe2 from first-principles [4]. However,metry breaking. Recently, it was found that thepreviously obtained result shows insulating state

noncentrosymmetric Nb/V/Ta superlattice showswith an extremely small gap, incompatible withsuperconducting transition at 4.4 K [12] and ex-

superconducting instability.

To overcome thishibits nonreciprocal charge transport in supercon-

program, we perform a GGA+U calculation using

ducting state. To identify the Rashba spin-orbit

WIEN 2k package [8] and find the insulator-metalinteraction in the artificial superlattice [Nb/V/Ta]ntransition by Coulomb interaction. The insulator-

without a center of inversion, we carried out den-metal transition is characterized by multi-orbitalsity functional theory calculations using WIEN2k

effect in Kondo lattice system.

The insulator-

package [8]. We created a slab [Nb/V/Ta]5 con-metal transition is a peculiar property of UTe2 ,taining 30 atoms, which corresponds to the fivewhich was not observed in other uranium-based

times stacking of two layers Nb, V, and Ta of asuperconductors nor in Kondo insulators. We canbcc structure. We obtained the Rashba splitting at

identify the three-dimensional winding number ω

the Fermi level, and the magnitude of the Rashba

In

this

anisotropic

project,

wetheoreticallysuperconductivity

in

UTe2

223

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

splitting is around 10 meV, which originate from

[8] P. Blaha, K. Schwarz, G. K. H. Mad-V atoms and is the origin of the nonreciprocity.

sen, D. Kvasnicka, J. Luitz, R. Laskowski,
F. Tran, and L. D. Marks, WIEN2k, An Augmented Plane Wave+Local Orbitals Program

References

for Calculating Crystal Properties (Karl[1] S. Ran, C. Eckberg, Q.-P. Ding, Y. Furukawa,heinz Schwarz, Techn. Universität, Wien,

T. Metz, S. R. Saha, I.-L. Liu, M. Zic,
H. Kim, J. Paglione, and N. P. Butch, Science
365, 684 (2019).

Austria, 2018).
[9] M. Sato, Phys. Rev. B 79, 214526 (2009).
[10] M. Sato, Phys. Rev. B 81, 220504(R) (2010).

[2] D. Aoki and J. Flouquet, J. Phys. Soc. Jpn.
83, 061011 (2014).

[11] L. Fu and E. Berg, Phys. Rev. Lett. 105,
097001 (2010).

[3] D. Aoki, K. Ishida, and J. Flouquet, J. Phys.
Soc. Jpn. 88, 022001 (2019).

[12] F. Ando, D. Kan, Y. Shiota, T. Moriyama,
Y. Shimakawa, T. Ono J. Magn. Soc. Japan[4] D. Aoki, A. Nakamura, F. Honda, D. Li,

43, 17-20 (2019).

Y. Homma, Y. Shimizu, Y. J. Sato, G. Knebel,
J.-P. Brison, A. Pourret, D. Braithwaite,
G. Lapertot, Q. Niu, M. Vališka, H. Harima,
and J. Flouquet, J. Phys. Soc. Jpn. 88, 043702(2019).
[5] G. Knebel, W. Knafo, A. Pourret, Q. Niu,
M. Vališka, D. Braithwaite, G. Lapertot, M. Nardone, A. Zitouni, S. Mishra,
I. Sheikin, G. Seyfarth, J.-P. Brison, D. Aoki,
and J. Flouquet, J. Phys. Soc. Jpn. 88, 063707(2019).
[6] G. Nakamine, S. Kitagawa, K. Ishida,
Y. Tokunaga, H. Sakai, S. Kambe, A. Nakamura, Y. Shimizu, Y. Homma, D. Li,
F. Honda, and D. Aoki, J. Phys. Soc. Jpn. 88,
113703 (2019).
[7] J. Ishizuka, S. Sumita, A. Daido, andY. Yanase: Phys. Rev. Lett. 123, 217001(2019).

224

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

Development and applications of continuous-timequantum Monte Carlo method in the thermodynamic limit
Takahiro OHGOE
Research Institute for Science and Engineering,
Waseda University, Okubo, Shinjuku-ku, Tokyo 169-8555

Quantum Monte Carlo methods are power

reach the 11th order in shorter CPU time bytools for investigating quantum many-bodyusing the CDet code. Furthermore, we alsosystems. However, they often offer from theobtained some results of density and Tan’snotorious negative-sign problem for fermioniccontact for imbalanced cases where thesystems and frustrated spin systems when theconventional quantum Monte Carlo method

system-size is increased. Recently, a novel

suffers from the negative-sign problem. The

determinant quantum Monte Carlo method wasdeveloped code will be useful for revealing

proposed,thefascinating physics such as FFLO or p-wave

thermodynamic limit [1]. The method is calledsuperconducting instability which are expectedconnected determinant Monte Carlo method

in imbalanced populations. This work was done

(CDet), and its applications to the Hubbard

in collaboration with Riccardo Rossi (Flatiron

model revealed non-trivial analytic structure of

institute), Felix Werner and Kris Van Houcke

the self-energy in strongly-correlated regime [2].

(ENS, Paris).

which

is

formulatedin

Therefore, it is a promising method for tacklingchallenging problems of quantum many-bodyReferences

systems.

[1] R. Rossi: Phys. Rev. Lett. 119 (2017) 045701

Here, we developed the original CDet code[2] F. Simkovic IV. and E. Kozik, Phys. Rev. B

which is applicable to the two-component

100 (2019) 121102

unitary Fermi gas. To show its advantage, the[3] R. Rossi, T. Ohgoe, K. Van Houcke, and F.

method was compared with the diagrammatic

Werner: Phys. Rev. Lett. 121 (2018) 130405

Monte Carlo method in the case of balanced[4] R. Rossi, T. Ohgoe, E. Kozik, N. Prokof’ev

population [3,4]. In the diagrammatic Monteand B. Svistunov, K. Van Houcke, and F.

Carlo method, the maximum expansion order

Werner: Phys. Rev. Lett. 121 (2018) 130406

which we could reach was 9, whereas we could

225

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

Kitaev spin liquid under time-dependent magnetic
field
Joji NASU
Department of Physics, Yokohama National University
Hodogaya, Yokohama 240-8501

Nonequilibrium phenomena in insulating
magnets such as real-time dynamics of magnetization and transport of spin excitations haveattracted considerable attention in condensed
matter physics. In the insulating magnets,
spin excitations are carried by magnons in thepresence of a magnetic order. On the other
hand, it was suggested that the spin transport
can be also governed by spinons, which are elementary excitations fractionalized from spins
in one-dimensional Heisenberg systems. Thiswas confirmed in the cuprate Sr2 CuO3 by thespin Seebeck experiments under the magnetic
field.
The Kitaev quantum spin model is anothercandidate of the magnetic systems associatedwith fractional quasiparticles. In this model,
the ground state is exactly shown to be a quantum spin liquid and the elementary excitations
are described by two kinds of quasiparticles:
itinerant and localized Majorana fermions. In
this work, we study spin dynamics under thetime dependent magnetic field to clarify therole of the fractional quasiparticles for thenonequilibrium phenomena.
Here, we examine the the real-time dynamics of the spin moments in the Kitaev model
under the time-dependent magnetic field. In
particular, we investigate the following twocases by using the Majorana mean-field theory.
First, 1) we study the eﬀect of the field quenchof the uniform magnetic field to understandfundamental properties of nonequilibrium dynamics of the fractional quasiparticles in the

Kitaev model. Next, 2) we consider the eﬀectof the time-dependent magnetic field pulse at
the edge of the cluster to mimic the spin injection into the Kitaev magnet.
With respect to 1), we find that two kindsof Majorana fermions are observed separately
in distinct time-scales in the time evolutionof the magnetization; the low-energy localized
Majorana fermions survive as long-lived excitations whereas the high-energy itinerant ones
are damped in the early stage [3].
With respect to 2), we find that the propagation of the spin excitations occurs even in
the absence of the magnetic field in the bulk
without spin polarization. We clarified that
this is attributed to the itinerant Majorana
fermions [2]. We also discuss the eﬀect of the
Heisenberg interaction and extension to theS = 1 Kitaev model [3] on this phenomenon.

References[1] J. Nasu, and Y. Motome: Phys. Rev. Research 1, 033007-1-16 (2019).
[2] T. Minakawa, Y. Murakami, A. Koga, and
J. Nasu: arXiv:1912.10599.
[3] A. Koga, T. Minakawa, Y. Murakami, and
J. Nasu: J. Phys. Soc. Jpn. 89, 033701-1-4(2020).

226

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

Magnetocaloric eﬀect on the Shastry–Sutherland
lattice antiferromagnet SrCu2(BO3)2
Subaru AKIMOTO and Yasuhiro H. MATSUDA
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581

We have studied the magnetocaloric eﬀect(MCE) of the Shastry–Sutherland (SS) lattice
antiferromagnet SrCu2 (BO3 )2 [1] with numerical simulations using the thermal pure quantum (TPQ) state. We used the HΦ as a numerical solver. The magnetic field dependence ofthe entropy (SM ) and specific heat (CM ) havebeen computed and compared with those experimentally obtained from the magnetic field
variation of the MCE at diﬀerent initial temperatures. Figures 1 (a) and (b) show the calculated SM and CM . The five trial curves andthe averaged result are shown for the SM and
CM , respectively.

SS lattice such as 1/4 and 1/3 plateaux [2].
The computed B–T (magnetic field – temperature) phase diagram shown in Fig. 2 is found tosemiquantitatively reproduce the experimental
results.

Figure 2: The calculated SM are plotted ina B–T plane with color graduation. Several
isentropic curves are also plotted.

Figure 1: The five trial curves and the averaged result of the calculation with TPQ state
are shown. The magnetic field is set to be 57 T
corresponding to the 1/3 plateau. (a) Calculated SM . (b) Calculated CM .

References

The distinct peak is shown in the temperature dependence of the CM at a magnetic field
and the peak positions at diﬀerent temperatures and magnetic fields would give the phaseboundaries of the quantum spin states of the[2] Y. H. Matsuda, N. Abe, S. Takeyama, H.
Kageyama, P. Corboz, A. Honecker, S. R.
Manmana, G. R. Foltin, K. P. Schmidt,
and F. Mila: Phys. Rev. Lett. 111 (2013)
137204.

[1] H. Kageyama, K. Yoshimura, R. Stern, N.
V. Mushnikov, K. Onizuka, M. Kato, K.
Kosuge, C. P. Slichter, T. Goto, and Y.
Ueda: Phys. Rev. Lett. 82 (1999) 3168.

227

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

Spin-Charge Separation in Charge Order Phase of
Organic Dirac Electron System α-(BEDT-TTF)2I3
Akito KOBAYASHI
Department of Physics, Nagoya University
Nagoya 464-8602

The organic Dirac electron system α(BEDT-TTF)2 I3 ehibits a phase transition between a massless Dirac electron phase and acharge-ordered phase under pressure as shownin Fig. 1. In the charge-ordered phase, thespin gap measured by NMR, the optical gap
in the optical conductivity, and the resistancegap determined by the DC resistance showcompletely diﬀerent values and pressure dependences. The mechanism of this novel spincharge separation has not been elucidated yet.
We have theoretically studied the mechanismof the spin-charge separation in α-(BEDTTTF)2 I3 .

tions by H. Kino and T. Tsumuraya. This factstrongly suggests that the type-I Dirac electronic state observed under hydrostatic pressure cannot be explained without the strong
correlation eﬀect. We will study the type-I
Dirac electronic state under hydrostatic pressure using the variational Monte Carlo method
in cooperation with Dr. Yoshimi of ISSP, andthen elucidate the mechanism of spin charge
separation of charge-ordered phase.

First, we have investigated the Dirac electron system with charge order in α-(BEDTTTF)2 I3 using the real-space-dependent meanfield approximation, and have shown that
there is a domain wall of the charge order thatspontaneously appears due to ferroelectricity.
We have shown that the wall-bound Dirac electrons form a one-dimensional conduction channel. Furthermore, we have shown that the linear response theory can explain the large difference between the optical gap and the resistivity gap[1].
Next, we have performed the first principles
calculation for α-(BEDT-TTF)2 I3 under hydrostatic pressure using Quantum ESPRESSO
in the Supercomputer Center, the ISSP. However, the result exhibits a type-II Diracsemimetal, which is diﬀerent from the experimental results. Similar results have been independently shown by the first-principles calcula-

Figure 1:
TTF)2 I3 .

Phase diagram of α-(BEDT-

References[1] D. Ohki, Y. Omori, and A. Kobayashi,
Phys. Rev. B 100, 075206 (2019).

228

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

Quantum Monte Carlo simulation and electronic
state calculations in correlated electron systems
Takashi YANAGISAWA
Electronics and Photonics Research InstituteNational Institute of Advanced Industrial Science and Technology (AIST)
AIST Central 2, 1-1-1 Umezono, Tsukuba 305-8568'!

!"#$%&

!"&

%"#

∆E/N

Our calculations are based on the optimized
variational Monte Carlo method[1, 2, 3, 4]. We
have investigated the ground-state phase diagram of the two-dimensional Hubbard modeland the two-dimensional d-p model. We performed parallel computation in Monte Carlocalculations. In order to reduce statistical errors, we carried out 100 ∼ 200 parallel calculations. Parallel computing is very importantto reduce Monte Carlo statistical errors.

%"&

AFI

&"#

SC
AF
SC

We used the wave function of anexp(−λK) − PG -type wave function. This
wave function is a very good many-body wave
function because the ground-state energy islowered greatly and the ground-state energy
is lower than that evaluated by any otherwave functions[2]. We can improve the wave
function systematically by multiplying by
operators PG and e−λK many times.

We show the phase diagram as a function of
the hole doping rate x in Fig. 1. The calculations were performed on a 10 × 10 lattice with
U/t = 18 and t′ = 0. There are three phases:
the antiferromagnetic insulator state (AFI),
the coexistent state of antiferromagnetism andsuperconductivity (AF+SC) and the d-wave
superonducting state (SC). Near half-filling for
approximatedly 0 ≤ x < 0.06, the ground stateis an AF insulator. The coexistent state exists
for 0.06 < x < 0.09. When the doping rate is
as large as x > 0.09, the ground state is d-wave
superconducting state. High-temperature superconductivity will occur in the SC phase inthis figure.

229

&"&
&"&&

&"&#

&"%&

&"%#

&"!&

&"!#

x

Figure 1:
Phase diagram for the twodimensional Hubbard model. The condensation energy per site as a function of the hole
doping rate x = 1 − ne on a 10 × 10 lattice.
We set t′ = 0 and U/t = 18. AFI indicatesthe antiferromagnetic insulating state and SC
denotes the d-wave SC phase. There is the coexistent state indicated as AF-S between thesestates.

References[1] T. Yanagisawa et al., J. Phys. Soc. Jpn.
67, 3867 (1998).
[2] T. Yanagisawa, J. Phys. Soc. Jpn. 85,
114707 (2016).
[3] T. Yanagisawa, J. Phys. Soc. Jpn. 88,
054702 (2019).
[4] T. Yanagisawa, Condensed Matter 4, 57(2019).

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

Dynamical mean-field theory + quantum Monte
Carlo approach to strongly correlated electron
systems with multi-orbitals
Shintaro HOSHINO
Department of Physics, Saitama University, Saitama 338-8570, Japan

We have also mapped out the phase diagramof the doped fulleride, which is also exprimentally realized [1]. To do this, we have analyzed the three-orbital Hubbard model with
antiferromagnetic Hund’s coupling by using
the dynamical mean-field theory. In this theory, the many-body problem in the lattice
is mapped onto the effective impurity problem without spatial periodicity, for which thecontinuous-time quantum Monte Carlo solver[5] has been used. Using the facilities of ISSP
supercomputer, we have calculated the susceptibilities which signal the instability of the normal state toward ordering state with spontaneously broken symmetries. The exemplary

230

0.14

orb. freez.
FO
SC
AFO

D=1, J/U=-0.25, n=1.5

0.12
0.1
T

The electrons with multiple orbital degrees
of freedom show a variety of intruguing phenomena as realized in heavy-electron materials
and also in d-electron systems such as ruthenate and iron-based superconductors. In addition, the fulleride superconductors are another example of the multiorbital electron systems [1], where the three-degenerate t1u molecular orbitals are occupied by electrons originating from alkaline doping. Due to the effectively sign-inverted Hund coupling by JahnTeller phonons, the electrons favor the low-spinstate and show superconductivity at low temperatures. Recently the experimental study
has revealed the existence of anomalous metallic state near the Mott insulator [2], and thepossibility of active orbital degrees of freedom
is pointed out theoretically [3]. The characteristic dynamical structure in the self-energy isidentified and can be a source of the anisotrpic
response function [4].

0.08
0.06
0.04
0.02
0

0

1

2

3

4
U

5

6

7

8

Figure 1: Temperature-interaction phase diagram.
The pahse boundaries for orbital ordering (FO) andsuperconductivity (SC) are shown, together withorbital freezing line.

results at the electron number n = 1.5 permolecule are shown in the figure. Here thephase diagram in the plane of temperature andCoulomb interaction is mapped out, where thehorizontal axis indicates the pressure or lattice constant dependence when one compare
the results with the experiments. We havedemonstrated that the superconductivity appears with dome-shaped transition temperatures, while the orbital ordering dominates atlarge Coulomb repulsive interaction.
References[1] For a review, see O. Gunnarsson, Rev. Mod.
Phys. 69, 575 (1997). [2] R.H. Zadik et al., Sci.
Adv. 1, e1500059 (2015). [3] S. Hoshino and P.
Werner, Phys. Rev. Lett. 118, 177002 (2017). [4]S. Hoshino, P. Werner, and R. Arita, Phys. Rev.
B 99, 235133 (2019). [5] E. Gull et al., Rev. Mod.
Phys. 83, 349 (2011).

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

Superconductivity and magnetic properties of the
Hubbard model
Atsushi Yamada
Department of Physics, Chiba University
Chiba 263-8522, Japana, Chiba 277-8581

Frustrated magnetic systems yield a rich variety of phenomena. For example, a spin liquidstate, which is a purely non-magnetic Mott insulator without spontaneously broken spatial
or spin symmetry, has attracted a lot of interest. This state is realized in geometrically frustrated systems like the charge organic transfer
salts κ-(BEDT-TTF)2 X[1] and Cs2 CuCl4 .[2] A
simple theoretical model of these compounds isthe Hubbard model on the an-isotropic triangular lattice, and spin liquid state is in fact
found in this model.[3] A spin liquid could
arise also in the intermediate coupling regionof strongly correlated systems between a semimetal and ordered state, because in this casea correlation-driven insulating gap might open
before the system becomes ordered. This possibility might be realized in the half-ﬁlled Hubbard model on the honeycomb lattice.

We have studied the magnetic and metalto-insulator transitions in the half-ﬁlled Hubbard model on the honeycomb lattice by variational cluster approximation using 10-site and16-site clusters as a reference system.[5] Thisapproach uses an exact diagonalization of the
Hubbard model deﬁned on these clusters andparts of numerical calculations were done using the computer facilities of the ISSP.Our results agree with recent large scale Quantum
Monte Carlo simulations.[4] We are currently
studying the eﬀect of the frustration on various
physically measurable quantities in the Hubbard model on the 15 -depleted square lattice.

References[1] Y. Shimizu, K. Miyagawa, K. Kanoda,
M. Maesato, and G. Saito, Phys. Rev.
Lett. 91, 107001 (2003); Y. Kurosaki, Y.
Shimizu, K. Miyagawa, K. Kanoda, and
G. Saito, Phys. Rev. Lett. 95, 177001(2005).
[2] R. Coldea, D.A. Tennant, A.M. Tsvelik,
and Z. Tylczynski, Phys. Rev. Lett. 86,
1335 (2001); R. Coldea, D.A. Tennant,and Z. Tylczynski, Phys. Rev. Lett. 68,
134424 (2003).
[3] T. Yoshioka, A. Koga, and N. Kawakami,
Phys. Rev. Lett. 103, 036401 (2009); P.
Sahebsara and D. Sénéchal, Phys. Rev.
Lett. 100, 136402 (2008); L.F. Tocchio,
H. Feldner, F. Becca, R. Valenti, and C.
Gros, Phys. Rev. B 87, 035143 (2013); A.
Yamada, Phys. Rev. B 89, 195108 (2014);
L.F. Tocchio, C. Gros, R. Valentı́, F.
Becca, Phys. Rev. B 89, 235107 (2014); A.
Yamada, Phys. Rev. B 90, 235138 (2014).
[4] S. Sorella, Y. Otsuka, and S. Yunoki, Sci.
Rep. 2, 992 (2012); F. F. Assaad andI. F. Herbut, Phys. Rev. X 3, 031010(2013); F. Parisen Toldin, M. Hohenadler,
F. F. Assaad, and I. F. Herbut, Phys. Rev.
B 91, 165108 (2015).
[5] A. Yamada, Int. J. Mod. Phys. B 30,
1650158 (2016).

231

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

Exploring low-rank tensor representation of
Matsubara Green’s functions
Hiroshi SHINAOKA
Department of Physics, Saitama University,
Saitama 338-8570

We have studied a low-rank representationof two-particle Green’s functions [1]. Twoparticle Green’s functions and the vertex functions play a critical role in theoretical frameworks for describing strongly correlated electron systems. However, numerical calculations at the two-particle level often suffer fromlarge computation time and massive memory
consumption. Example includes parquet approximation and diagrammatic non-local extensions of dynamical mean-field theory. For
example, a simulation for the single-band twodimensional Hubbard model based on the parquet approximation requires more than a TB
of memory even at a high temperature.
Thus, a compact representation of two-particle
Green’s functions is highly desired.
In Ref. [1], we introduced a sparse sampling
scheme in the Matsubara frequency domain as
well as a tensor network representation for twoparticle Green’s functions. The sparse sampling is based on the intermediate representation basis [2, 3, 4] and allows an accurate extraction of the generalized susceptibility froma reduced set of Matsubara frequencies. The
tensor network representation provides a system independent way to compress the information carried by two-particle Green’s functions.
We used System B at the ISSP supercomputer center to demonstrate efficiency of thepresent scheme for calculations of static anddynamic susceptibilities in single- and twoband Hubbard models in the framework of dynamical mean-field theory.

Figure 1 shows the dynamical lattice susceptibility computed by Bethe-Salpeter equation (BSE) based on dynamical mean-field theory for a two-band Hubbard model. In solving BSE, we compressed the local two-particle
Green’s function using the tensor-network representation. The uncompressed data takes up
700 MB even on the sparse grid, while thetensor network representation takes up only
330 kB for D = 60. Here, D is a bond dimension in the tensor network.
In future studies, we will use the new techniques to study real compounds by means of
the non-local extension of the dynamical meanfield theory and the parquet approximation.
This report is based on the collaboration
with D. Geffroy, M. Wallerberger, J. Otsuki,
K. Yoshimi, E. Gull, J. Kuneš.
Sz

1.0

y

y

0.8
0.6
0.4
0.2
0.0

X

M

X

M

X

M

Figure 1: Dynamical lattice susceptibility computed for a two-band Hubbard model. Reproduced from Ref. [1] under the CC BY 4.0 license.

References[1] H. Shinaoka, D. Geffroy, M. Wallerberger,
J. Otsuki, K. Yoshimi, E. Gull, J. Kuneš,232

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

SciPost Phys. 8, 012 (2020)[2] H. Shinaoka, J. Otsuki, M. Ohzeki and K.
Yoshimi, Phys. Rev. B 96, 035147 (2017)[3] H. Shinaoka, J. Otsuki, K. Haule, M.
Wallerberger, E. Gull, K. Yoshimi and M.
Ohzeki, Phys. Rev. B 97, 205111 (2018)[4] J. Otsuki, M. Ohzeki, H. Shinaoka andK. Yoshimi, Sparse modeling in quantum
many- body problems, J. Phys. Soc. Jpn.
89, 012001 (2020)

233

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

Photoinduced dynamics and high harmonic
generation in magnetism
Sumio Ishihara
Department of Physics, Tohoku University
Sendai 980-8578

Nonequilibrium susceptibility in Floquet
statesRecent great progress in laser light technology significantly promotes scientific researchin strongly coupled light-matter systems. Oneof the fascinating phenomena induced by thestrong light-matter coupling is the appearance of the Floquet states, in which a timeperiodic electromagnetic field behaves as photons and hybridizes with electrons. We examine the spin/charge susceptibility in photoinduced Floquet states [1]. We consider a freeelectron system coupled to a fermionic heatbath. We adopt the Hamiltonian given byH∑ ††ks εk cks cks +ν bν bν()∑+ ksν Vn u c†ks bν + b†ν cks ,=2) On the other hand, in the case of smallΩ comparable to the bandwidth, the multiple
Floquet bands cross the Fermi level of the bathand the additional peaks emerge at momenta
corresponding to nesting vectors between
the Floquet-band Fermi surfaces (Fig. 2).

∑(1)

where c†ks is a creation operator of an electronwith momentum k and spin s(=↑, ↓), and b†νis that of a fermion in a bath with a quantum
number ν.
The susceptibility is formulated on the basisof the Floquet Green s function method, andits perturbative expression is derived from aseries expansion. Through the numerical calculations, we have the following results.
1) In the case of the large electric-field frequency Ω in comparison with the electron
bandwidth, we demonstrate that the bandwidth reduction due to the dynamical localization eﬀect enhances the susceptibility anda deformation of the Fermi surface which depends on light polarization shifts the peak positions of the susceptibility (Fig. 1).

234

Figure 1: The susceptibility at q = (0, 0)
and (π, π) as a function of the light amplitude
A0 [1].
Optical response in the excitonic insulating state
The excitonic insulating (EI) state was proposed more than a half century ago in semiconductors and semimetals and has been studied
intensively on both the experimental and theoretical sides. We study the optical responses in
EI states [2]. We adopt the two-orbital Hubbard model with finite energy diﬀerence between the orbitals defined by
H = ∆

∑

(niaσ − nibσ ) +

iσ

+

∑

∑

i

α

(U

∑

tα c†iασ cjασ

⟨ij⟩ασniα↑ niα↓ + U ′ nia nib )Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Figure 2: The susceptibility, the momentum
distribution function; and the spectral functionat ω = 0, under the linearly polarized light of
A0 = 0.4 and 1.6 [1].
+ J

∑
iσσ ′

+ I

c†iaσ c†ibσ′ ciaσ′ cibσ∑ †ciα↑ c†iα↓ ciᾱ↓ ciᾱ↑ ,

(2)

iα

where ciασ is an annihilation operator of anelectron with spin σ(=↑, ↓) and the orbitalα(= a, b) at site i, and niασ = c†iασ ciασ is theparticle number operator.
The variational cluster approach method
was applied to the two-orbital Hubbard modelwith a finite energy diﬀerence between the orbitals. In the analyses of the ground state, we
considered possibilities of the low-spin state,
the high-spin antiferromagnetic state, the EI
state. The optical conductivity spectra were
formulated by the Green’s function method,where the vertex correction is taken into account. Figure 3 shows the optical conductivity
spectra. In the EI phase, a new peak aroundω/t ∼ 10 and a continuum around 3 < ω/t < 8appear. The intensity of the EI component increases with increasing J/t in the EI phase. It
is shown that the intensity of this structure at
low temperatures almost follows the EI orderparameter.
These researches have been performed withcollaboration with A. Ono, H. Li, M. Naka
and J. Ohtsuki. The numerical calculations in

235

Figure 3: Optical conductivity spectra for several values of the Hund coupling J. The
red solid lines show results including the vertex corrections, while the black dashed results show results without the vertex corrections. [2].
these works have been done using the facilitiesof the Supercomputer Center, the Institute for
Solid State Physics, the University of Tokyo.

References[1] A. Ono and S. Ishihara: Phys. Rev. B
100, 075127 (2019).
[2] H. Li, J. Otsuki, M. Naka and S. Ishihara
Phys. Rev. B 101 125117 (2020).

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

Mechanism of superconductivity and electron-holeasymmetry in molecular conductors
κ-(BEDT-TTF)2X
Hiroshi WATANABE
Cluster for Pioneering Research (CPR), RIKEN
2-1 Hirosawa, Wako-shi, Saitama 351-0198

The family of quasi two-dimensional molecular conductors κ-(BEDT-TTF)2 X has been
extensively studied as a typical example ofstrongly correlated electron system. Depending on the monovalent anion X, they show
various quantum phases such as antiferromagnetic (AF) and spin-liquid dimer-Mott insulators, and superconductivity (SC) [1]. Owing
to the similarities in the experimental phase
diagrams, κ-(BEDT-TTF)2 X system is often
compared with high-Tc cuprates, which exhibit
Mott metal-insulator transition and SC [2].
Recently, carrier doping in κ-(BEDT-TTF)2 X
has been realized using electric-double-layertransistor [3], and the direct comparison withcuprates becomes possible. The experimental
result shows the electron-hole doping asymmetry, reminiscent of the high-Tc cuprates.

electron-hole doping asymmetry in the phase
diagram where AF spin order, diﬀerent patterns of charge order, and SC compete witheach other. Hole-doping stabilizes the dimertype AF phase and promotes SC with dxy -wave
symmetry, which has similarities with high-Tc
cuprates. In contrast, in the electron-dopedside, geometrical frustration destabilizes the
AF phase and the enhanced charge correlationinduces another SC with extended-s+dx2 −y2 wave symmetry. Our results disclose the mechanism of each phase appearing in filling-control
molecular Mott systems, and elucidate how
physics of diﬀerent strongly-correlated electrons are connected, namely, molecular conductors and high-Tc cuprates.

Here, we theoretically study the doping effect in κ-(BEDT-TTF)2 X. We consider a fourband extended Hubbard model including onsite (U ) and intersite Coulomb interactions(Vij ) with κ-type geometry with the largest
hopping integral tb1 as an energy unit. The
ground state properties are studied with thevariational Monte Carlo (VMC) method. The
Gutzwiller-Jastrow type wave function is used
for the VMC trial wave function. The systemsizes for calculation are from 288 (2×12×12)
to 1152 (2×24×24).
Figure 1 shows the ground-state phase diagram of the extended Hubbard model for
κ-(BEDT-TTF)2 X [4]. We find significant

236

Figure 1: Ground-state phase diagram ofthe extended Hubbard model for κ-(BEDTTTF)2 X [4].

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

References[1] K. Kanoda: J. Phys. Soc. Jpn. 75 (2006)
051007.
[2] R. H. McKenzie: Science 278 (1997) 820.
[3] Y. Kawasugi et al., Sci. Adv. 5 (2019)eaav7282.
[4] H. Watanabe, H. Seo, and S. Yunoki: Nat.
Commun. 10 (2019) 3167.

237