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

3.3

Strongly Correlated
Quantum Systems

183

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

Machine learning and ab initio analysesof cuprate high-temperature superconductors
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 cuprate

spin liquids (QSL). It was achieved by using thesuperconductors remains one of the challenges

variational Monte Carlo tools (mVMC) [2]in physics. In this project, we have continued

combined

efforts to clarify physics and mechanism ofmachine developed by us [3,4], tensor network

high-temperature superconductivity, particularly[5,6] and the Lanczos diagonalization. In thefor the copper oxides.

This is a combined

J1J2 Heisenberg model defined on the 2D square

research activity of three different approaches:lattice with the nearest neighbor exchange(1) Parameter search of strongly correlated

interaction (J1) and the next nearest neighborelectron models by using high-accuracy solvers

exchangefor quantum many-body lattice Hamiltonians

established the existence of the QSL in the(2) Parameter-free ab initio studies of real

region 0.49<J2/J1<0.54 [7]. The analyses werematerials by using multi-scale ab initio scheme

elaborated by comparing the accuracy of ourfor correlated electrons (MACE) [1]solver with other tools in the literature. The(3) Data science approaches to expose hidden

accuracy and reliability of our solver were

physicalbyshown to be the best among existing solvers both

combining extended experimental data such as

for the ground states and the excitations. The

by angle resolved photoemission spectroscopy

reliability and robustness of the analyses was(ARPES), quasiparticle interference (QPI) by

confirmed by (i) the size insensitivity of thescanning tunnel microscope (STM) and resonant

correlation-ratio crossing points to identify theinelastic X-ray scattering (RIXS) with the help

QSL-antiferromagnetic transition and the QSLof applied mathematical and information science

totechnology such as machine-learning tools.

insensitivity about the choice of the initialquantities

andmechanisms

with

(J2),

the

we

valence-bond-solid

restrictedhaveBoltzmannfurthertransition

firmly(ii)

theIn the first approach, we have continued tovariational parameters for the wavefunction (iii)

clarify physics of electron fractionalization andbetter size extrapolation of the level crossing

formation of topological states such as quantum

point to infer the ground state phase boundary

184

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

from the level crossing points of the excitation

Accurate algorithms and tools for quantum

spectra, where larger size calculation has been

dynamics were also sought for [9]. A method to

added to enhance the accuracy of thecalculate spectral functions measurable by the

extrapolation.

ARPES representing the imaginary part of theWe have also been investigating the quantum
ground

state

of

the

Heisenberg

developed in addition to the two-particle

antiferromagnet on the pyrochlore lattice [8].

dynamical structure factor for spin and charge

After investigation of various trial wave

[10]. The application to the Hubbard model

functions, we find that the singlet ground state

revealed the d-wave gap structure for the firsthas a linear dispersion at vanishing momentum

time in the superconducting state as shown inwith a finite energy gap to triplet excitations.

Fig.1. It has shown an unrealistically large gap

Furthermore, it has turned out that the groundof the typical Hubbard model implying anstate is fully degenerate for all the point groupoversimplified nature of the Hubbard model as asymmetrymodel of the cuprate superconductors.

represented

S= 1/2electron single-particle Green’s function was

by

the

irreducible

symmetry group, which implies an emergence of

For the ab initio approaches (2), we havethe unprecedented type of quantum spin liquid.

reproduced the experimental phase diagram of
HgBa2CuO4+y [11] by solving its ab initio lowenergy effective Hamiltonian without adjustableparameters [12]. Thanks to this success, in thepresent project in 2020, more thorough andsystematic analyses for several different copperoxide compounds have started. Derivation of abinitio low-energy effective Hamiltonians for theseries of multi-layer compounds including Bi
and solving with refined mVMC tool is underway.
The third data science approach (3) has been
elaborated for the machine learning analysis of
the ARPES data. The comparison with theFig.1: Spectral function along around

literature has clarified the reliability of oursymmetry line. (f) is closest to theanalyses, which shows the emergence ofsymmetry line showing a large d-wave

prominent peak structures in the normal andgap at (π,0) [10].

anomalous parts of the self-energies and their

185

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

cancellation in the Green’s function resulting indoes not occur in the absence of thethe direct invisibility in ARPES, while thefractionalization, it can be used as the stringent

prominent peak is the origin of the hightest for the occurrence of the fractionalization.

temperature

The

This is a combined report for E project 2020-

accuracy and reliability of the present machine

Ea-0005 and 2020-Eb-0005 as well as sharedlearning method were confirmed by several

project for Fugaku project.

superconductivity[13].

This series of work has been done inbenchmark tests, which successfully reproduced

collaboration with T. Misawa, M. Hirayama, K.

the expected exact results and the established

Ido, Y. Nomura, Y. Yamaji, R. Pohle, M.

analyses in the literature. It further establishedCharlebois, J-B. Moree, F. Imoto, A. Fujimori

the robustness against unavoidable experimental

and T. Yoshida. The work is also supported bynoise and extrinsic contributions to spectral

JSPS Kakenhi 16H06345, the RIKEN Center for

functions such as background effects. The origin

Computational Science under the HPCI projectof the failure of previous studies which did not(hp200132) This project has used the software

find the prominent structure in the self-energies

HPhi, mVMC and RESPACK.

References

are exposed by faithfully following the[1] M. Imada and T. Miyake: J. Phys. Soc. Jpn. 79(2010) 112001[2] T. Misawa et al.: Comput. Phys.
Commun. 235 (2019) 447[3] Y. Nomura et al.: Phys. Rev. B.96 (2017)
205152[4] G. Carleo et al.: Nat. Commun. 9 (2018) 5322
[5] H.-H. Zhao et al.: Phys. Rev. B.96 (2017)085103[6] A. S. Darmawan et al.: Phys. Rev. B 98 (2018)
205132.
[7] Y. Nomura and M. Imada: arXiv:2005.14142.
[8] R. Pohle, Y. Yamaji and M. Imada: unpublished.
[9] M. Charlebois and M. Imada: Phys. Rev. X 10
(2020) 041023.
[10] K. Ido, M. Imada and T, Misawa: Phys. Rev.
B 101 (2020) 075124[11] T. Ohgoe et al., Phys. Rev. B 101 (2020)
045124.
[12] M. Hirayama et al.: Phys. Rev. B 99 (2019)245155.
[13] Y. Yamaji et al.: arXiv:1903.08060[14] M. Imada, arXiv:2103.10016.

assumptions by the previous studies and by
showing the error contained in the assumptions.
More thorough studies of the integrated
spectroscopy analysis have been conducted in
the combination of ARPES, and QPI or RIXS.
By using the ARPES data and their machinelearning analyses, a two-component fermion
model

is

constructedtorepresentthefractionalization of electrons supported by theprevious theoretical and experimental studies.
Then the two-component Hamiltonian was
analyzed to predict the RIXS data [14]. The
prediction shows a substantial enhancement of
the RIXS intensity in the superconducting phasein comparison to the normal phase, if theelectron fractionalization correctly describes thelowenergydynamics

ofthecupratesuperconductors. Since such an enhancement186

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

Numerical studies of quantum spin liquids byquantum mutual information
Youhei YAMAJI
GREEN, National Institute for Materials Science
Namiki, Tsukuba-shi, Ibaraki, 305-0044

Recent studies on quantum spin liquids(QSLs) revealed that a wide variety of the
QSLs is characterized by emergent fractionalized excitations and their topological nature,
instead of an absence of spontaneous symmetry breaking at zero temperature. Entanglement entropy has been proposed to distinguishtopologically trivial and non-trivial states. Although the entanglement entropy is hard tomeasure in experiments, it is applicable to any
many-body states.

the density matrix with a canonical thermal
pure quantum (cTPQ) state [1, 2] as, ρ̂pure =|Φ(β)ihΦ(β)|/hΦ(β)|Φ(β)i, where |Φ(β)i =e−β Ĥ/2 |Φ(0)i and |Φ(0)i is a normalized random vector. As emphasized in the literature,ρ̂pure is different from ρ̂: any integer power ofρ̂pure is equal to ρ̂ as E[ρ̂mpure ] = ρ̂, where E
is the average over distribution of the randomP
vector |Φ(0)i and Z(β) = n e−βEn . In contrast, the reduced ensity matrix constructedfrom ρ̂pure for the subsystem X (Y ), ρ̂X
pure
X/Y

Instead of the entanglement entropy, mutual
information (MI), a generalization of entanglement entropy, is useful at finite temperatures, where temperature scales characterizing
the fractionalization may emerge. The MI isgiven as relative entropy between the density
matrix and a product of reduced density matrices of a given system: When the system isdevided into two subsystems, X and Y , the MI 
I is given as I = tr ρ̂ ln ρ̂ − ln(ρ̂X ⊗ ρ̂Y ) ,where ρ̂ is the density matrix and ρ̂X (ρ̂Y ) is
the reduced density matrix for the subsystem
X (Y ) by taking a partial trace over the complementary subsystem Y (X).
Naively, the computational and memory
costs for I are higher than those for obtaining an eigenstate. However, typical state approaches enable us to estimate I with O(NH )computational and memory costs, where NH
is the Hilbert space dimension of the target
system, although the density matrix consisting of a typical pure state shows unphysical properties, as follows. When one define(ρ̂Ypure ), satsfies ρ̂X/Y = E[ρ̂pure ].
X/YBecause ρ̂pure and ρ̂ ln ρ̂ = −S(β), whereS(β) is entropy at inverse temperatures β, arecalculated by the cTPQ with O(NH ) computational and memory costs, ρ̂ ln(ρ̂X ⊗ ρ̂Y ) can
be calculated with the similar costs. Here, wenote that ρ̂ appears only once in ρ̂ ln(ρ̂X ⊗ ρ̂Y ).
The MI of an ab initio hamiltonian for Na2 IrO3
has been simulated and the temperature dependence will be published. The reduced density matrices are also useful to clarify thenature of correlated electrons such as hightemperature superconductivity [3].

References[1] S.Lloyd: arXiv:1307.0378; Ph.D. Thesis,
Rockefeller University (1988)[2] S. Sugiura and A. Shimizu: Phys. Rev.
Lett. 111 (2012) 010401.
[3] A. Iwano: master’s thesis, The Universityof Tokyo (2020).

187

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

Numerical studies of quantum spin liquid candidates
by highly accurate ab initio effective hamiltonians
Youhei YAMAJI
GREEN, National Institute for Materials Science
Namiki, Tsukuba-shi, Ibaraki, 305-0044

A Kitaev’s spin liquid candidate, α-RuCl3 ,has attracted considerable attention.
Although several studies on the effective hamiltonian of α-RuCl3 , there is no consensus on theeffective hamiltonian and no effective hamiltonian consistently explains experimentally observed specific heat and inelastic neutron scatterings. While these effective hamiltonians significantly deviate from the Kitaev model, theruthenium halide shows experimental signatures of Majorana excitations. To refine theeffective hamiltonian of the ruthenium halide
and to clarify whether the ruthenium halide isin the vicinity of the Kitaev model or not, wetheoretically study the series of the ruthenium
halides, α-RuH3 (H=Cl, Br, I).
Recent developments of deriving abinitio effective hamiltonians by using localized Wannier orbitals and constrained random phase or
GW approximations [1, 2] enables us to study
spin-orbit coupled Mott insulators and correlated (semi)metals. Regardless of the details ofthese methods, as the ionic radius of the halogen atom increases from Cl to I, the correlations in the ruthenium halides become weaker.
The detailed analyses on the effective abinitio
spin hamiltonians and itinerant electron hamiltonians are performed by using HΦ [3] andmVMC [4]. The spin hamiltonian seems to beinvalid for H =I and, at least, t2g hamiltonians
are required to describe the end member of theruthenium halides.
To characterize the excitation spectra of correlated electrons, we have also developed nu-

merical procedures to simulate abinitio spectral functions for the many-body Schrödinger
equation, instead of the effective lattice hamiltonians, by combining variational Monte Carlo
methods and Krylov subspace methods [5]. To
examine the applicability of the present algorithm, it is applied to the single particle spectral functions of liquid helium 3. The selfenergy is derived, which explicitly demonstrate
the strongly correlated nature of liquid helium
3 in an abinitio fashion.

References[1] M. Charlebois, J.-B. Morée, K. Nakamura, Y. Nomura, T. Tadano, Y. Yoshimoto, Y. Yamaji, T. Hasegawa, K. Matsuhira, M. Imada: arXiv:2103.09539.
[2] J.-B. Morée, S. Biermann, Y. Yamaji, and
M. Imada: unpublished.
[3] M. Kawamura, K. Yoshimi, T. Misawa,
Y. Yamaji, S. Todo, and N. Kawashima:
Compt. Phys. Communi. 217 (2017) 180.
[4] T. Misawa, S. Morita, K. Yoshimi, M.
Kawamura, Y. Motoyama, K. Ido, T. Ohgoe, M. Imada, and T. Kato: Compt.
Phys. Communi. 235 (2019) 447.
[5] N. Okada: Dynamical Variational Monte
Carlo Approaches for Excitation Spectraof Many-Body Fermions in Continuum
Space, master’s thesis, The University ofTokyo (2020).

188

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

Systematic analysis of ab initio low-energy eﬀective Hamiltonians
for Pd(dmit)2 molecular conductors
Takahiro Misawa
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581

Quantum spin liquids, which do not show
any symmetry breaking even at zero temperature, have attracted much interest since theirelementary excitations in the quantum spin liquids are expected to show exotic elementary
excitations such as Majorana particles. Because the exotic elementary excitations may
be useful for next-generation devices such as
quantum computers, a huge amount of works
on searching and identifying the quantum spin
liquids in solids has been done in a decade.

software package RESPACK [5, 6], we evaluate the transfer integrals and two-body interactions such as the Coulomb interactions. As
a result, we have found that the anisotropy of
the transfer integrals and correlations eﬀects
systematically change by changing cations.

Among several candidates of the quantum spin liquids, the molecular solids β ′ X[Pd(dmit)2 ]2 (X represents a cation) oﬀers
an ideal platform for realizing the quantum
spin liquid induced by the geometrical frustration in the magnetic interactions because
the geometrical frustration can be systematically controlled by changing cations X. In experiments, it has been proposed that the Neeltemperatures are systematically controlled bychanging cations X and the quantum spin liquid realizes in X =EtMe3 Sb [1, 2].

Moreover, we have analyzed the low-energy
eﬀective Hamiltonians using the exact diagonalization method [7, 8]. From the numerical exact analyzes of the low-energy eﬀective Hamiltonians, we have shown the significant reduction of the antiferromagnetic ordered moment occurs around X =EtMe3 Sb.
This reduction is consistent with the experimentally observed quantum spin liquid behavior in X=EtMe3 Sb. We have also shown thatthe reduction is induced by both the geometrical frustration in the hopping integrals andoﬀ-site Coulomb interactions. This result indicates that accurate evaluation of the microscopic parameters in the Hamiltonians is essential for reproducing the quantum spin liquidbehavior in the dmit-salts.

In this project, to identify the microscopic
origin of the quantum spin liquid found in
X =EtMe3 Sb, we have performed systematicand comprehensive ab initio derivations of lowenergy eﬀective Hamiltonians for available 9
compounds of β ′ -X[Pd(dmit)2 ]2 (X = Me4 Y ,EtMe3 Y , Et2 Me2 Y and , Y = As, Sb, andP) [3]. In the derivation the low-energy eﬀective Hamiltonians, we ﬁrst obtain the global
band structure for dmit-salts using Quantum
ESPRESSO [4]. Then, using the open-source

In addition to the ab initio study for thedmit-salts, we have also analyzed the longrange spin transport in the topological Diracsemimetals [9] using the real-time evolutionof the quantum systems. Furthermore, usingHΦ [7, 8], we have analyzed the magnetization process of the antiferromagnetic Heisenberg model on the kagome lattice [10]. We
have also developed an open-source library forthe shifted Krylov subspace method (Kω) [11]and an Ab initio tool for derivation of eﬀective

189

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

low-energy model for solids (RESPACK) [5, 6].

References[1] R. Kato, Chem. Rev. 104, 5319 (2004).
[2] K. Kanoda, R. Kato, Annu. Rev. Condens. Matter Phys. 2, 167 (2011).
[3] T. Misawa, K. Yoshimi, T. Tsumuraya,
Phys. Rev. Research 2, 032072 (2020).
[4] P. Giannozzi, O. Andreussi, T. Brumme,
O. Bunau, M. B. Nardelli, M. Calandra,
R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, et al., Journal of Physics: Condensed Matter 29, 465901 (2017).
[5] https://sites.google.com/view/kazuma7k6r[6] K. Nakamura, Y. Yoshimoto, Y. Nomura,T. Tadano, M. Kawamura, T. Kosugi,
K. Yoshimi, T. Misawa, Y. Motoyama,Computer Physics Communications 261,
107781 (2021).
[7] M. Kawamura, K. Yoshimi, T. Misawa,
Y. Yamaji, S. Todo, N. Kawashima, Comput. Phys. Commun. 217, 180 (2017).
[8] https://www.pasums.issp.u-tokyo.
ac.jp/HPhi/en/[9] Y. Araki, T. Misawa,
arXiv:2103.06519.

K. Nomura,[10] T. Misawa, Y. Motoyama, Y. Yamaji,
Phys. Rev. B 102, 094419 (2020).
[11] T. Hoshi, M. Kawamura, K. Yoshimi,
Y. Motoyama, T. Misawa, Y. Yamaji,
S. Todo, N. Kawashima, T. Sogabe,
Computer Physics Communications 258,
107536 (2021).

190

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

Study of phase formation, transport phenomena and effects ofexceptional points in strongly correlated quantum systems
Norio KAWAKAMI
Department of Physics,
Kyoto University, Kitashirakawa, Sakyo-ku, Kyoto 606-8502
equilibrium approaches.[1]We have studied phase formation andtransport phenomena in strongly correlated(2) Non-Hermitian phenomena

quantum systems. Especially, we have focused

Quasiparticles

on phenomena such as the formation ofwithfinitetheirlifetimesThis

non-Hermitian

provided a concise review of the non-Hermitian
topological band structures for quantum many-

We have discussed the circumstances andbody systems in equilibrium as well as theirrequirements for the emergence of non-classification.[2]Hermitian phenomena in equilibrium singleparticle properties of strongly correlated systems
relationshiplifetime.

new insights into correlated systems. We(1) Open quantum systems

theirfiniteperspective on equilibrium systems provides

becomes essential.

andGreen's

Hermitian topological phenomena because ofpoints when the non-Hermitian aspect of
bands

by

functions of equilibrium systems exhibit non-magnetic phases and the effect of exceptional
quasiparticle

described

to

Non-Hermitian phenomena offer a novel

non-Hermitian

approach to analyze and interpret spectra in thephenomena in open quantum systems (OQS)

presence of interactions. Using the density-

experiencing gain and loss. While the necessity

matrix renormalization group (DMRG), weof postselection has limited the practicality ofdemonstrated the existence of exceptional points

studying non-Hermitian phenomena in the OQS,

for the one-particle Green's function of the 1D

in closed equilibrium systems, the single-

alternating Hubbard chain with chiral symmetry,particle Green's function is described by anwith a corresponding Fermi arc at zero

effective non-Hermitian Hamiltonian without

frequency in the spectrum. They are robust andthe need of postselection. By describing

can be topologically characterized by the zeroth

quasiparticles of strongly correlated systems as

Chern number. This effect illustrates a casean OQS, we demonstrated that the non-where the temperature has a strong impact in 1D

Hermitian Hamiltonian arising in OQS andbeyond the simple broadening of spectral

Green’s functions are equivalent. We have also

features. Furthermore, we demonstrated that

demonstrated the necessity of considering theexceptional points appear even in the two-memory effect in the quantum master equation

particle Green's function (charge structure

by comparing the spectral function in thefactor), where an effective Hamiltonian is

Hubbard model computed by an OQS anddifficult to establish.[3]191

Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo(3) Two-dimensional materials
Themagnetic

phosphorene

for device applications.[5]properties

nanoribbons

of

black

have

been

(4) Many-body localization
Strong disorder in low-dimensional strongly

investigated using static and dynamical mean-

correlated systems has been discussed to induce

field theory. Besides confirming the existence offerromagnetic/antiferromagnetic

many-body

edge

the

to precisely determine the critical strength ofinteraction strength to a novel incommensurate

disorder at the localization phase transition inmagnetic phase. Furthermore, we demonstrated

one or higher spatial dimensions. We havethat the difference of the ground state energies ofconsidered the Fock-space localization of many-

the AFM and FM phase is exponentially small,body wavefunctions in an all-to-all interacting

making it possible to switch between both states

modelby a small external field.[4]offermions

anddemonstratedaquantitative agreement of the numericallyspin-dependentobtained moments of eigenstate wavefunctions

thermoelectric effect of graphene flakes with

as well as the spectral statistics to analyticalmagnetic

regime.

predictions.[6]Employing static and dynamical mean-field

References

edges

the

is

Hamiltonian. However, it has been an open issue

unit-cells find a phase transition at weak

investigatedwhich

localization of all many-body eigenstates of themagnetism, our detailed calculations using large

We

localization,in theballistic

theory, we first showed that magnetism appears[1] Y. Michishita, and R. Peters, Phys. Rev. Lett.

at the zigzag edges for a window of Coulombinteractions that increases significantly with

124, 196401 (2020).

increasing flake size. We then used the Landauer[2] T. Yoshida, R. Peters, N. Kawakami, and Y.

formalism in the framework of the non-

Hatsugai, Prog. Theor. Exp. Phys. 12A109

equilibrium

Green'sfunction

methodto(2020).

calculate the spin and charge currents in[3] R. Rausch, R. Peters, and T. Yoshida, New J.

magnetic hexagonal graphene flakes by varying
the temperature of the junction for different flake

Phys. 23, 013001 (2021).

sizes. While in non-magnetic gated graphene,[4] J. Vahedi and R. Peters, Phys. Rev. B 103,

the temperature gradient drives a charge current,

075108 (2021).

we observe a significant spin current for[5] T. T. Phung, R. Peters, A. Honecker, G. T. dehexagonal graphene flakes with magnetic zigzag
edges. Specifically, we showed that in the

Laissardiere, and J. Vahedi, Phys. Rev. B 102,

"meta" configuration of a hexagonal flake

035160 (2020).

subject to weak Coulomb interactions, a pure[6] F. Monteiro, T. Micklitz, M. Tezuka, and A.

spin current can be driven just by a temperature

Altland, Phys. Rev. Research 3, 013023 (2021).

gradient in a temperature range that is promising192

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

Theoretical study of correlated electron systems
with strong spin-orbit coupling
Yukitoshi MOTOME
Department of Applied Physics,
The University of Tokyo, Bunkyo, Tokyo 113-8656
We have theoretically studied a variety of

relevant mechanism to understand recent new

intriguing phenomena in correlated electron

generation of the topological spin crystals [9].

systems with spin-orbit coupling, ranging from(ii) Kitaev quantum spin liquids: By using theMott insulators to metals. During the last fiscal

quantum Monte Carlo simulation in theyear, we have achieved substantial progress on

Majorana

the following topics (project numbers: 2020-

investigated the effect of randomness on the

Ca-0059 and 2020-Cb-0038). We summarize

Kitaev model for the thermodynamics andthe main results below.

thermal transport [10] and the spin dynamics(i) Topological spin crystals: We clarified that[11]. We also studied the finite-temperatureinterplay between spin, charge, and orbital

phase transitions in a variety of three-degrees of freedom plays a crucial role indimensional extensions of the Kitaev model

stabilizing the magnetic hedgehog crystals[12]. In addition, by using ab initio based

recently found in B20 chiral magnets [1]. We

calculations, we performed a systematic study

also showed that the spin-charge interplay givesof Pr-based f-electron compounds as new

rise to a square-lattice type skyrmion crystals

Kitaev candidates [13]. Summarizing our recent[2]. Based on the results, we made astudies of the physics of Kitaev spin liquids andcollaboration with the experimental group to

the material design of the Kitaev magnets, weclarify the importance of the spin-charge

wrote two review papers [14,15].

coupling in GdRu2Si2 hosting the square(iii) Spin-orbit physics not requiring the spin-

skyrmion crystal [3,4]. We also studied theorbit coupling: We showed that κ-type organic

effect of anisotropy in a triangular lattice

antiferromagnets exhibit an anomalous Hall

system [5]. In addition, we investigated the

response because of an effective spin-orbit

effects of the spatial anisotropy [6], the phase

coupling generated by glide symmetry breaking

shift [7], and the twist angle [8] on theby the antiferromagnetic order [16]. We also

topological spin crystals. We wrote a review

showed that a similar mechanism predicts apaper on the itinerant frustration which is aspin current generation in perovskites with the193

fermionrepresentation,

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

C-type antiferromagnetic order [17].

(arXiv:2103.10647).

(iii) Multipole physics: We proposed channel-[10] J. Nasu and Y. Motome, Phys. Rev. B 102,

selective non-Fermi liquid behavior in the two-

054437 (2020).

channel Kondo lattice model under a magnetic[11] J. Nasu and Y. Motome, preprint

field [18]. We also clarified optical Hall

(arXiv:2103.10549).

responses in spin-orbit coupled metals with

[12] T. Eschmann, P. A. Mishchenko, K.

magnetic cluster multipole orders [19].

O'Brien, T. A. Bojesen, Y. Kato, M. Hermanns,(iv) Collaboration with experimental groups: In

Y. Motome, and S. Trebst, Phys. Rev. B 102,

addition to [3,4], we made combined studies

075125 (2020).

between our theory based on the ab initio[13] S.-H. Jang, R. Sano, Y. Kato, and Y.

calculations and experiment on topologicalMotome, Phys. Rev. Materials 4, 104420

materials [20,21].

(2020).
[14] Y. Motome, R. Sano, S.-H. Jang, Y. Sugita,

References

and Y. Kato, J. Phys: Condens. Matter 32,[1] S. Okumura, S. Hayami, Y. Kato, and Y.

404001 (2020), [Special Issue on QuantumMotome, Phys. Rev. B 101, 144416 (2020).

Spin Liquids].

[2] S. Hayami and Y. Motome, Phys. Rev.

[15] 求 幸年，那須 譲治, 固体物理 55, 297

B 103, 024439 (2021).

(2020).

[3] Y. Yasui et al., Nat. Commun. 11, 5925

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

(2020).

Yanagi, Y. Motome, and H. Seo, Phys. Rev.

[4] 東京大学大学院工学系研究科 プレスリリ

B 102, 075112 (2020).

ース, "電荷密度波で電子構造の「ねじれ」をス[17] M. Naka, Y. Motome, and H. Seo, Phys.

イッチする −トポロジカル状態の高速制御に向Rev. B 103, 125114 (2021).

けた新しい指針を開拓−".

[18] K. Inui and Y. Motome, Phys. Rev. B 102,

[5] S. Hayami and Y. Motome, Phys. Rev.

155126 (2020).

B 103, 054422 (2021).

[19] T. Sato, Y. Umimoto, Y. Sugita, Y. Kato,[6] K. Shimizu, S. Okumura, Y. Kato, and Y.

and Y. Motome, Phys. Rev. B 103, 054416

Motome, Phys. Rev. B 103, 054427 (2021).

(2021).

[7] S. Hayami, T. Okubo, and Y. Motome,[20] N. Mitsuishi et al., Nat. Commun. 11,

preprint (arXiv:2005.03168).

2466 (2020) [selected as FOCUS].

[8] K. Shimizu, S. Okumura, Y. Kato, and Y.

[21] M. Onga et al., Nano Letters 20, 4625

Motome, preprint (arXiv:2009.14569).

(2020).

[9] S. Hayami and Y. Motome, preprint

194

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

Study on newly discovered nickelate superconductivity
Yusuke NOMURA
RIKEN Center for Emergent Matter Science, 2-1 Hirosawa, Wako, 351-0198
Hubbard-type material) and cuprates (charge-

The recent discovery of superconductivity

transfer-type material).

in the doped nickelate Nd0.8Sr0.2NiO2 has
attracted much attention because the nickelates(B) Then, it is of great interest to investigate

may serve as a cuprate analog superconductor

the strength of magnetic exchange coupling J in[1]. Indeed, the electronic structure of thethe d9 nickelates. Note that a large value of J of

nickelate is similar to that of the cuprates in theabout 130 meV is a characteristic feature of thesense that the strongly correlated x2−y2 orbital

cuprates. We study the J value in theoretically-is present. However, as a crucial difference

designed RbCa2NiO3 and A2NiO2Br2 (A: afrom the cuprates, which shows Mott insulating

cation with the valence of +2.5) (Fig. 1). We

behavior, RNiO2 (R=Nd, Pr) is not a Mott

show that these nickelates have a sizeable

insulator due to the carrier doping from themagnetic exchange coupling as large as about

rare-earth layer (self-doping) [2].

80-100 meV, which is not far smaller than thatof the cuprates [4].

In the present study, we investigate thefollowing fundamental questions: (A) Is it
possible to design nickelates whose electronic
structure is more similar to that of the cuprates?
(B) Can we expect large magnetic exchange
coupling J in the cuprate-analog nickelates?
For (A), using a concept of “block layers”,
we perform a systematic materials design oflayered nickelates. Then, we find several
dynamically stable cuprate-analog nickelates[3]. In such nickelates, the self-doping is absent,
and hence the x2−y2 orbital becomes half-filled.
Therefore, the correlation effect will induce

Fig. 1: Crystal structure of RbCa2NiO3 andMott insulating behavior in such d9 nickelates.

A2NiO2Br2 (A: a cation with the valence of +2.5).

Using the d9 nickelates, we can perform a faircomparison between the nickelates (Mott-

195

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

References[1] D. Li et al., Nature 572, 624 (2019)[2] Y. Nomura et al., Phys. Rev. B 100, 205138(2019) [Editor's suggestion][3] M. Hirayama and T. Tadano et al., Phys.
Rev. B 101, 075107 (2020)[4] Y. Nomura et al., Phys. Rev. Research 2,
043144 (2020).

196

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

Studies on superconductivity in multiorbital systems with an incipient bandand designing of new nickelate superconductors
KAZUHIKO KUROKI
Department of Physics, Osaka University
1-1 Machikaneyama, Toyonaka, Osaka, 560-0043, Japan

SUPERCONDUCTING MECHANISM OF A
NEW CUPRATE SUPERCONDUCTOR
BA2 CUO3+δ

Recently, Li et al.[1] reported high-Tc superconductivity in a new cuprate Ba2 CuO3+δ witha K2 NiF4 -like layered structure. There, they ﬁnd
several unique features which strongly suggest thatthe material is a diﬀerent type of cuprate superconductor that opens up a new paradigm. Namely,a large amount of oxygen vacancies are present
within the CuO2 planes, and a great amount ofholes are doped which should cause a large deviation of the Cu valence from 2+. Also, the distance
between Cu and the apical O is shorter than thein-plane Cu-O distance, so that the octahedron iscompressed along the c axis. This should result ina crystal ﬁeld very diﬀerent from that of the ordinary cuprates, namely, the 3d3z2 −r2 orbital lifted inenergy above the 3dx2 −y2 orbital. These ﬁndings
suggest that the mechanism of superconductivity
in the Ba2 CuO3+δ may be considerably diﬀerent
from that for the conventional cuprates.
In 2020 ﬁscal year, we have studied possible superconducting mechanisms of Ba2 CuO3+δ [2]. In
order to study the superconducting mechanism of
the new superconductor, a realistic model Hamiltonian is required, which has to start with a determination of the crystal structure. For the 2-1-3 composition in particular, the chain structure is known
to be stable in an actual material Sr2 CuO3 . As another possibility within the 2-1-3 composition, we
considered a Lieb lattice type structure. First principles total-energy calculation performed with the
VASP code shows that the chain and the Lieb lattice structures are close in the total energy, so thatthe latter may also be considered as a candidate.
Focusing on these two structures, we have calculated their band structures, and extracted theWannier orbitals to construct multiorbital Hubbard models, i.e, a two-orbital model for the chain
structure and a six-orbital model for the Lieb typestructure. We have applied ﬂuctuation exchange
approximation to the eﬀective models, and discuss
the possibility of superconductivity by solving thelinearized Eliashberg equation. We have shownthat s-wave and d-wave pairings closely competewith each other and, more interestingly, that a coexistence of intra- and interorbital pairings arises.

197

We also reveal an intriguing relation of the Liebmodel with the two-orbital model for the usual
K2 NiF4 -type cuprate where a close competition between s- and d-wave pairings is known to occur.
We have further shown that s±-wave superconductivity is strongly enhanced when the d3z2 −r2band is raised in energy so that it become nearly
incipient with the lower edge of the band close
to the Fermi level within a realistic band ﬁlling
regime. The enhanced superconductivity in thepresent model is in fact shown to be related to an
enhancement found previously in the bilayer Hubbard model with an incipient band. Namely, we
can show that, by orbital basis transformation[3],the two-orbital model can be transformed to the bilayer Hubbard model when all the intra and interorbital interactions have the same magnitude. In
reality, the interactions are diﬀerent in magnitude,
but we ﬁnd that the similarity between the twomodels holds to some extent. The strong enhancement of s±-wave superconductivity in the bilayer
model with an incipient band has been shown by
a number of previous studies, including our numerical study[4] based on multi-variable variational
Monte Carlo method[5, 6].

DESIGNING NICKELATESUPERCONDUCTORS

Inspired by the proposed superconducting mechanism for Ba2 CuO3+δ , in 2020 ﬁscal year,
we have also theoretically designed unconventional nickelate superconductors with d8 electron
conﬁguration[7]. Although the materials that wedesign apparently have nothing to do with bilayer
systems at ﬁrst glance, electronic structure of thesematerials resemble that of the bilayer Hubbard
model, thereby strongly enhancing superconductivity. It is known that bilayer Hubbard model withappropriately large interlayer electron hopping exhibits s±-wave superconductivity with extremely
high-Tc , but it is generally diﬃcult to realize sucha situation in actual bilayer-type materials due tovarious restrictions regarding the atomic orbitals.
We have adopted a completely diﬀerent strategy torealize the desired situation, by which we design amixed-anion nickelates as candidates for new highTc superconductors.

Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo2
mately equivalent to large interlayer hopping in abilayer system, as schematically depicted in Fig.2.
Applying the ﬂuctuation exchange approximation
to the ﬁve orbital model of these materials, we haveshown that the maximum superconducting transition temperature of some of them may be evenhigher than that of the high-Tc cuprates.
FIG. 1.
Band structures of (a) La2 NiO4 , (b)Ca2 NiO2 Cl2 and (c) (b) Ca2 NiO2 H2 (taken from
Ref.[7]).
two-orbital

bilayerorbital-1

antibonding

EF

DE

EF

orbital-22t⊥bonding

FIG. 2. Schematic figure of the correspondence between
the two-orbital model and the bilayer model (takenfrom Ref.[7]).

8More speciﬁcally, we have considered d nickelates with K2 NiF4 structure, where the apical oxygens are replaced by halogens or hydrogens. The
lattice parameters are determined by structural optimization using VASP. The The band structures of
two of the proposed candidates, Ca2 NiO2 Cl2 and
Ca2 NiO2 H2 , are shown in Fig.1 together with thatof a reference material La2 NiO4 . The key point
here is that the large crystal ﬁeld splitting of theorbital energy levels, induced by replacing the apical oxygens by halogens or hydrogens, is approxi-

198

We have further found that our theory may also
be related to another newly discovered superconductor (Nd,Sr)NiO2 [8]. In this superconductor, thenickel electron conﬁguration is expected to be close
to d9 , but this might be aﬀected, e.g., by the reduction process during the synthesis of the material. Since the band structure of (Nd,Sr)NiO2 itself, due to the absence of apical oxygens, resembles that of the proposed mixed anion-nickelates,(unintentional) reduction of the electron number in(Nd,Sr)NiO2 may also result in a similar s±-wave
superconductivity.

[1] W.M. Li et al., Proc. Natl. Acad. Sci. USA 116,
12156 (2019).
[2] K. Yamazaki et al. Phys. Rev. Res. 2, 033356(2020).
[3] H. Shinaoka et al., Phys. Rev. B 92, 195126 (2015).
[4] D. Kato and K. Kuroki, Phys. Rev. Res. 2, 023156(2020).
[5] D. Tahara and M. Imada, J. Phys. Soc. Jpn. 77,
114701 (2008).
[6] T. Misawa et al., Comp. Phys. Comm. 235, 447(2019).
[7] N. Kitamine, M. Ochi and K. Kuroki, Phys. Rev.
Res. 2, 042032(R).
[8] D. Li et al., Nature (London) 572, 624 (2019).

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

Bulk-edge correspondence for non-Hermitian
topological systems
Tsuneya Yoshida
Department of Physics,
1-1-1 Tennodai, Tsukuba, Ibaraki, JAPAN 305-8571
In this year, we have studied topologicalphases in non-Hermitian systems [1-6].
In particular, by employing numerical
diagonalizationofmatrix,

have

we

a

non-Hermitianelucidatedthetopological properties of non-Hermitian
fractional quantum Hall state [1].

Fig.1 Spectrum of the Liouvillian matrix

Specifically, we have analyzed a two-describing fractional non-Hermitian quantum

dimensional open quantum system with- Hall states.

two-body loss. Such a system is described

dimensions and symmetry, e.g., one

by Lindblad equation. Vectorizing thedimensional open quantum systems withdensity matrix, we have mapped theinversion symmetry.

Liouvillian to the Liouvillian matrix.

We have also analyzed classical systems

Diagonalizing
analyzed

thismatrix,

topological

we

have

in term of topological band theory[7-9].

properties.

Our

Our

analysishave

elucidated

that

analysis elucidate that for a system with

topological phenomena are observed a

two-body

variety of systems beyond quantum

loss,

the

Liouvilliangapremains open (see Fig. 1). In addition,
introducing
number,we

topologyofthepseudo-spin

haveelucidated

systems.

Chern

References

thatfractional[1] T. Yoshida, K. Kudo, H. Katsura, andquantum Hall states is maintained for

Y. Hatsugai, Phys. Rev. Research 2

this open quantum systems.

033428 (2020)

Our approach of the characterization[2] T. Yoshida, Phys. Rev. B 103 125145

can be also generalized to other cases of(2021).

non-Hermitian

199

Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[3] H. Araki, T. Yoshida, and Y. Hatsugai,(2020).

J. Phys. Soc. Jpn. 90 053703 (2021).

[7] T. Yoshida, T. Mizoguchi, and Y.

[4] P. Delplace, T. Yoshida, Y. Hatsugai,

Hatsugai arXiv:2012.05562.

arXiv: 2103.08232.

[8] T. Yoshida, T. Mizoguchi, Y. Kuno and[5] R. Rausch, R. Peters, and T. Yoshida,

Y. Hatsugai arXiv: 2103.11305

New J. Phys. 23 013011 (2021).

[9] T. Yoshida and Y. Hatsugai, Sci. Rep.

[6] T. Yoshida, T. Mizoguchi, and Y.

11 888 (2021).

Hatsugai Phys. Rev. Research 2 022062.

200

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

Numerical analyses on quantum spin liquids in strongly
correlated electron systems
Kota IDO
Institute for Solid State Physics,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
β’- Pd(dmit)2 salts have attracted attention

find that the energy of quantum spin liquids

as candidate materials where a quantum spinusing our proposed wavefunction is much

liquid emerges [1]. This system consists of thelower than those using the conventional pairing

dimers of two Pd(dmit)2 molecules and cations.

wavefunction. In addition, we performed

Since the dimers form an anisotropic triangularvariance extrapolations of energies by using thelattice, this system is considered as a two-

power Lanczos method [4]. As the results, we

dimensional correlated electron system withsucceed in reproducing the stability of thestrongAlthoughquantum spin liquid obtained by using thenumerical studies play an important role in

fixed-node method [5]. Our results suggest thatinterpretingthethis approach is a new efficient way to analyze

candidate materials for quantum spin liquids, it

the ground state in strongly correlated electron

is difficult to accurately analyze large two-systems with the geometrical frustration and todimensional systems with strong frustration andclarify the nature of the spin liquids.

geometricalfrustration.

experimental

resultsinmany-body correlations.

In addition, we have analyzed the groundIn this study, to perform the ground state

states of an extended Kitaev model using theanalysis on strongly correlated electron systems

VMC method. It is a future issue to study notwith geometrical frustration, we improved theonly the ground states but also nonequilibrium

accuracy of the trial wavefunction of thephenomena in quantum spin liquids.

variational Monte Carlo (VMC) method. To
achieve this goal, we proposed a two-

Referencescomponent[1] K. Kanoda, R. Kato, Annu. Rev. Condens.

pairing

wavefunction

as

the

fermionic wavefunction. The difference from

Matter Phys. 2, 167 (2011).

the conventional pairing wavefunction is that[2] G. Carleo and M. Troyer, Science 355, 602

the strength of the pairing is dependent on the(2017).

number of the local density. By using this[3] Y. Nomura, A. S. Darmawan, Y. Yamaji, andwavefunction

recently

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

developed neural network correlators [2, 3], we[4] E. S. Heeb, T. M. Rice, Z. Phys. B 90, 73

performed benchmarks for the t-t’ Hubbard(1993).

model on the square lattice and the triangular[5] L. F. Tocchio, A. Parola, C. Gros, and F.

lattice with one-dimensional anisotropy. We

Becca, Phys. Rev. B 80, 064419 (2009).

combined

with

201

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

Search for multiple-Q magnetic orders in systems
with bond-dependent anisotropic interactions
Satoru Hayami
Department of Applied Physics, The University of Tokyo, Tokyo 113-8656

Noncoplanar spin textures have attracted
considerable attention for potential applications to next-generation electronic and spintronic devices. In the project with numbers:
2020-Ca-0002 and 2020-Cb-0005, we have theoretically investigated a plethora of noncoplanar spin textures in correlated electron systems. We have presented the main results thisyear below.
(i) Magnetic hedgehog lattice: A magnetic
hedgehog lattice (HL), which is characterized
by a periodic array of magnetic monopoles
and antimonopoles, is one of the noncoplanar multiple-Q states in the three-dimensionallattice system. The emergence of the HLs
has been suggested in the noncentrosymmetricmetal MnSi1−x Gex [1, 2] by experiments, but
their microscopic origins have not been fully
clarified owing to the unconventional shortmagnetic periods of the spin textures. We haveinvestigated the stability of the HLs by analyzing an eﬀective spin model derived from anitinerant electron model [3]. By using the variational calculations and simulated annealing,
we obtained two types of HLs in the wide rangeof model parameters in the ground-state phase
diagram. The results indicate the importanceof the interplay between the DzyaloshinskiiMoriya(DM)-type interaction by the spin-orbit
coupling and the multiple-spin interactions by
the spin-charge coupling.
(ii) Square skyrmion crystal: We have studied an instability toward a square skyrmion
crystal (SkX), which consists of two helices with equal weight, in a centrosym-

202

metric tetragonal lattice system [4].
We
considered an eﬀective spin model from anitinerant electron model on a square lattice. Reflecting the itinerant nature of electrons, the spin model exhibits the bilinear and biquadratic interactions in momentum space [5]. We also introduced fourfoldsymmetric bond-dependent anisotropic andeasy-axis anisotropic interactions to stabilize
the multiple-Q states. By performing simulated annealing, we have shown that the square
SkX appears by the synergy eﬀect among thebiquadratic, bond-dependent anisotropic, andeasy-axis anisotropic interactions in an external magnetic field. We have also shown thatthe magnetic phase diagram obtained by simulated annealing well reproduces the experimental results for the SkX hosting material
GdRu2 Si2 [6, 7].
In addition, we have investigated the vorticity and helicity of the square SkX in centrosymmetric itinerant magnets. Owing to theabsence of the DM interaction in centrosymmetric magnets, the vorticity and helicity of
the SkX are arbitrary. In other words, theNéel SkX, Bloch SkX, and two anti-SkXs aredegenerate. In such a situation, we foundthat the anisotropic Ruderman-Kittel-KasuyaYosida (RKKY) interaction depending on the
wave vectors lift the degeneracy of the SkXs
with the diﬀerent vorticity and helicity [8].
Our result opens a possibility of controlling thetypes of the SkXs through the electronic band
structures.
(iii) Triangular skyrmion and meron crys-

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

tals: We have studied the SkXs and meron
crystals in itinerant hexagonal magnets [9].
On the basis of the eﬀective spin model withthe bilinear and biquadratic interactions in
momentum space [5], we examined the eﬀectof two magnetic anisotropies on the stabilityof the multiple-Q states: The one is the local single-ion anisotropy and the other is the
sixfold-symmetric bond-dependent anisotropy.
In the case the single-ion anisotropy, we foundthat both the SkXs with the skyrmion number of one and two, which are denoted as thensk = 1 SkX and nsk = 2 SkX, appear under the single-ion anisotropy, although theirstability against the single-ion anisotropy is
diﬀerent with each other. For the nsk = 2
SkX, it is found that the stable region underthe easy-axis anisotropy is wider than that for
the easy-plane anisotropy at zero field. Meanwhile, the critical magnetic field to destabilize the nsk = 2 SkX is larger for the easyplane anisotropy than that for the easy-axis
anisotropy. For the nsk = 1 SkX, there is a
drastic eﬀect of the single-ion anisotropy on its
stability; the nsk = 1 SkX is very weak (strong)against the easy-plane(-axis) anisotropy.
Under the bond-dependent anisotropy, we
also obtained both the nsk = 1 and nsk = 2
SkXs in the wide range of model parameters.
Notably, we found a variety of chiral magnetic
states with nonzero scalar chirality, which isdistinct from the SkXs. Especially, we obtained two types of the meron crystals: One
is the nsk = 1 meron crystal, which consists ofone meron-like and three antimerion-like spin
textures in the magnetic unit cell and the otheris the nsk = 2 meron crystal consisting of fourmeron-like spin textures in the magnetic unitcell.
The eﬀective spin model in itinerant magnets will provide a deep understanding of themechanism of the multiple-Q states, whichhave been recently observed in experiments, inan eﬃcient way. Indeed, we reproduced themagnetic phase diagrams in Gd3 Ru4 Al12 [10]203

and CeAuSb2 [11], which host the multiple-Q
states, by analyzing the extended eﬀective spinmodel.

References[1] T. Tanigaki et al., Nano Lett. 15 (2015)5438.
[2] Y. Fujishiro et al., Nat. Commun. 10(2019) 1059.
[3] S. Okumura, S. Hayami, Y. Kato, and
Y. Motome, Phys. Rev. B 101 (2020)
144416.
[4] S. Hayami and Y. Motome, Phys. Rev. B
103 (2021) 024439.
[5] S. Hayami, R. Ozawa, and Y. Motome,
Phys. Rev. B 95 (2017) 224424.
[6] N. D. Khanh et al., Nat. Nanotech. 15(2020) 444.
[7] Y. Yasui et al., Nat. Commun. 11 (2020)
5925.
[8] S. Hayami and R. Yambe, J. Phys. Soc.
Jpn. 89 (2020) 103702.
[9] S. Hayami and Y. Motome, Phys. Rev. B
103 (2021) 054422.
[10] M. Hirschberger, S. Hayami, and Y.
Tokura, New J. Phys. 23 (2021) 023039.
[11] S. Seo et al., Commun. Phys. 4 (2021) 58.

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

Time-Dependent DMRG Study of Spectral Shape in
the Optical Conductivity of Two-Dimensional
Hubbard Model
Takami TOHYAMA
Department of Applied Physics, Tokyo University of Science, Tokyo 125-8585

The parent compound of cuprate superconductors is a Mott insulator in two dimensions.
The optical conductivity in the compound has
been measured and its spectral features havebeen discussed experimentally. However, theoretical understanding of its spectral shape isstill far from complete. This is due to theoretical diﬃculty caused by strong coupling of
localized spins and holon/doublon created byphotoexcitations. For example, it is unclear
whether there is a magnetically induced excitonic peak at the absorption edge or not [1, 2].
In order to clarify spectral properties of the
optical conductivity in two-dimensional Mott
insulator, we use time-dependent densitymatrix renormalization group (tDMRG) technique for a half-ﬁeld Hubbard model withsecond-neighbor hopping (t-t′ -U model) andcalculate time-dependent current induced byelectric ﬁeld applied along the x direction, fromwhich we can obtain the optical conductivity.
In our tDMRG the time-evolution operator
necessary for the calculation of time-dependent
wave function is evaluated by using a kernel
polynomial method [3], which is an approach
good for lattices more than one dimension.
We use a 6×6 square lattice of the t-t′ -Umodel with U/t = 10 under open boundary
condition [4]. The number of states kept in
tDMRG procedure is 4000. We note that we
also use an 8×8 lattice to check size dependence for a t-U model and ﬁnd that the difference of the optical conductivity in the twolattices is small. For the 6×6 lattice withoutt′ , the optical conductivity shows a prominent
peak at the absorption edge followed by continuous spectral weight distribution above thepeak position. The peak intensity depends onthe value of t′ , exhibiting a maximum of thepeak intensity at t′ = 0 as a function of t′ .
Since t′ introduces a diagonal antiferromagnetic exchange interaction, the eﬀect of frustration competing with the nearest-neighborexchange interaction appears. The frustrationeﬀect should not be dependent on the sign oft′ . This is the case of the present result. We
can say that magnetic interactions contributeto the peak at the absorption edge, indicating
magnetic origin of the peak, i.e., magnetic exciton. However, we note that there is no gap
between the peak and continuum mentioned
above. This is diﬀerent from a standard exciton where the exciton is deﬁned as a boundstate with well-deﬁned binding energy.

References[1] T. Tohyama, Y. Inoue, K. Tsutsui, and
S. Maekawa: Phys. Rev. B 72 (2005)
045113.
[2] X.-J. Hana et al.: New J. Phys. 18 (2016)
103004.
[3] K. Shinjyo and T. Tohyama: Phys. Rev.
B 103 (2021) 035141.
[4] K. Shinjyo, Y. Tamaki, S. Sota, and T.
Tohyama, in preparation.

204

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

Development and application of DFT+DMFT
software DCore
Hiroshi SHINAOKA
Department of Physics, Saitama University,
Saitama 338-8570

To make DFT+DMFT available to moreusers and to promote the development of acommunity of users and developers, we needan open-source program package with a userfriendly interface.
In this project, we have developed an opensource program, DCore v3.0.0 [1], that implements DMFT. This program features an interface based on text and HDF5 files, allowing
DMFT calculations of tight-binding models tobe performed on predefined lattices as well as
first-principles models constructed by external
DFT codes through the Wannier90 package.
Furthermore, DCore provides interfaces to avariety of quantum impurity solvers such as
quantum Monte Carlo codes developed in the
ALPS projects. DCore is implemented on the
top of the TRIQS Python library (Toolbox for
Research on Interacting Quantum Systems).
DCore v3.0.0 supports Python 3.x and TRIQS
3.0.

We used DCore to generate hybridization
functions for a five-orbital 2 × 2 cluster impurity model for LaAsFeO [2]. The impurity
model involves 40 spin orbitals. The hybridization functions computed by DCore are shownin Fig. 1. The hybridization functions havelarge off-diagonal elements and decay exponentially in the intermediate representation (IR)basis. We found that the exact hybridization
functions with large off-diagonal elements canbe fitted accurately with 332 bath sites using
sparse modeling techniques.
This report is based on the collaboration
with Y. Nagai, J. Otsuki, M. Kawamura, N.
Takemori, K. Yoshimi.
La( c ) W it h o u t r a n d o mO

p o t e n t i a Asl

Fe2
Fe1

1 0 1
∆i j ( l )Dynamical mean-field theory (DMFT) has
become a standard theoretical tool for studying strongly correlated electronic systems. In
a DMFT calculation, an original lattice model
is mapped to an effective Anderson impurity problem whose bath degrees of freedom
are self-consistently determined. Although
the DMFT was originally proposed for solving models such as Hubbard models, it canbe combined with density functional theory(DFT) based on ab initio calculations for describing the electronic properties of strongly
correlated materials. This composite framework, called DFT+DMFT, has been applied
to various types of materials.

1 0∆1 1 ( l )∆1 6 ( l ) ×−4∆1 7 ( l ) ×

1 0 −2
1 0 −4

1 0 −91 0 −1 4
0

1 0

2 0

3 0
l

4 0

5 0

6 0

Figure 1: Illustration of the DMFT selfconsistency cycle implemented in DCore

References[1] H. Shinaoka, J. Otsuki, M. Kawamura, N. Takemori, K. Yoshimi,
arXiv:2007.00901v3[2] H. Shinaoka and Y. Nagai, PRB 103,
045120 (2021)

205

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

Classical Monte Carlo study of J1-J2 Heisenberg
antiferromagnet on the kagome lattice
Hiroshi SHINAOKA
Department of Physics, Saitama University,
Saitama 338-8570

Frustrated XY antiferromagnets have been
investigated extensively in the past decades
due to exotic phenomena induced by geometrical frustration. The classical nearest-neighbor
XY antiferromagnet on the kagome lattice isa typical two-dimensional frustrated magnet.
This model hosts a macroscopically degenerate
ground-state manifold. Small perturbations
such as thermal fluctuations lift the macroscopic degeneracy. This gives rise to various
exotic phenomena. An interesting question is
how next-nearest-neighbor interaction affects
the competition of degenerated states at finite
temperature.

Figure 1: q = 0 and
tions

√

3×

√

3 spin configura-based on the collaboration with F. Kakizawa
and T. Misawa.

In this project, we investigated the finitetemperature phase diagram of classical J1 J2 XY antiferromagnets on the kagome lattice through extensive Monte Carlo simulations [1, 2] using loop updates. We implemented the code in Julia, which is a modern
programming language for scientific computing. The code is parallelized using MPI.jl.
We found that a weak antiferromagnetic J2
induces an unconventional first-order transition at low temperatures. Furthermore, in thevicinity of the first-order transition, we find an
octupole ordered phase between two competing quasi-long-range orders with different spin
configurations (See Fig. 1). At antiferromatic
J2 = −0.03, the ground state is the q = 0
state. As shown in Fig. 2, there is a firstorder transition between the high-temperatureparamagnetic phase and the q = 0 quasi-longrange-ordered phase. At J2 = −0.005, theBerezinskii-Kosterlitz-Thouless transition appears at a higher temperature. This report is

Figure 2: Temperature (T ) dependence of thespecific heat C computed for J2 = −0.03 andJ2 = −0.005

References[1] F. Kakizawa, T. Misawa, H. Shinaoka,
APS March Meeting 2021, Y39.00012.
[2] 柿澤文哉, 三澤貴宏, 品岡寛, J1 –J2 古典
反 強 磁 性XY カ ゴ メ 模 型 に お け る 新 奇相転移、 日本物理学会第76回年次大会（2021年), 12pC1-8.

206

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

Research of Quantum Critical Points
Emerging between Two-Channel Kondo andFermi-Liquid States
Takashi HOTTA
Department of Physics, Tokyo Metropolitan University
1-1 Minami-Osawa, Hachioji, Tokyo 192-0397

In this research, we have investigated theemergence of quantum critical points (QCP)
near a two-channel Kondo phase by evaluating
an f -electron entropy of a seven-orbital impurity Anderson model hybridized with three (Γ7and Γ8 ) conduction bands with the use of a numerical renormalization group (NRG) method.
First we consider the local f -electron Hamiltonian Hloc . For the purpose, we introduce one
f -electron state, which is deﬁned by the eigenstate of spin-orbit and crystalline electric ﬁeld
(CEF) potential terms. Under the cubic CEF
potential, we obtain Γ7 doublet and Γ8 quartet
from j = 5/2 sextet, whereas we ﬁnd Γ6 doublet, Γ7 doublet, and Γ8 quartet from j = 7/2
octet. By using those one-electron states asbases, we describe Hloc 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 λ isthe spin-orbit coupling of f electron. In thisresearch, we set λ = 0.1 and 0.11 for Pr and Ndions, respectively. 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 forthe angular momentum ℓ = 3. Note that thesixth-order CEF potential term B60 does notappear for j = 5/2, since the maximum size ofthe change of the total angular momentum isless than 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.
Now we consider the matrix element I of
the Coulomb interaction. Here we skip the details of the derivation of I, but they are expressed with the use of four Slater-Condon 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 ruleinteraction among f orbitals. it is reasonable
to set U as 1 eV in this research.
In Fig. 1, we show the local CEF groundstate phase diagram for n = 3, obtained from207

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

1.5

0.015

G6 doublet
W=10-3
x=0.0

0.005

1

W=10
x=-1.0-3

0.0-0.005-0.01-0.015-0.2

Simp

B60 [meV]0.01

G8 quartet

n=3
U=1
l=0.11

0.5 log 2-0.1

0.0
B40 [meV]

0.1

0.2

the diagonalization of Hloc . For n = 3, theground-state multiplet for B40 = B60 = 0
is characterized by total angular momentum
J = 9/2. Under the cubic CEF potentials, thedectet of J = 9/2 is split into three groups asone Γ6 doublet and two Γ8 quartets.
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 by∑εkc†kµτckµτ +k,µ,τ

+ Hloc ,

log f

0.5

Figure 1: Local CEF ground-state phase diagram on the plane of (B40 , B60 ) for n = 3.
Red line denote the trajectory of B40 = W x/15
and B60 = W (1 − |x|)/180 in the range of−1 ≤ x ≤ 0 for W = 10−3 .

H =

L=8.0
x=0.0
M=5000x=-0.3x=-0.372
x=-0.4x=-0.5n=3
U=1
l=0.11
W=10-3
V=0.7∑

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
states. Here we use Λ = 8 and M = 5, 000.
Note that the temperature T is deﬁned as T =DΛ−(N −1)/2 in the NRG calculation, where N

0 -9
10

-8

10

-7

-6

-5

-4

10 10 10 10
temperature T [eV]

-3

10

-2

10

Figure 2: Entropies Simp vs. temperature forn = 3 in the range of −1 ≤ x ≤ 0 with W =√10−3 . Note that ϕ is given by ϕ = ( 5 + 1)/2.
is the number of the renormalization step andD is a half of conduction band width. Here we
set D = 1 eV.
Here we brieﬂy discuss the typical results forn = 3 [2]. In Fig. 2, we show entropies Simp vs.
temperature, when we change the CEF parameters from x = 0 (Γ6 doublet) to x = −1 (Γ8quartet) for n = 3, W = 0.001, and V = 0.7.
For x = 0 and −0.3, we obtain the magnetic
two-channel Kondo phase, while for x = −0.4
and −0.5, the screened Kondo phase appears,
since the Γ8 quartet eﬀectively expressed by
S = 3/2 spin is screened by three conduction
electrons. At x = −0.372 between two-channel
Kondo and Fermi-liquid phases, we ﬁnd an entropy plateau of log ϕ. In this case, we could
not observe the residual entropy of log ϕ at low
enough temperatures, but the signal of QCP is
considered to be obtained.

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

208

-1

10

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

Disorder eﬀect on Kitaev quantum spin liquids
Joji Nasu
Department of Physics, Yokohama National University
Hodogaya, Yokohama 240-8501, Japan

Quantum spin liquid has been actively studied in strongly correlated electron systems. It
is a state that does not show magnetic order down to the lowest temperature, whichappears in insulating magnets and thought
to appear owing to the cooperation of strong
quantum ﬂuctuations and many-body eﬀects.
Thus far, one considers geometrical frustration, which is essential to realize this state in
the Heisenberg quantum spin model in the lattice structure. However, in this situation, one
usually suﬀers from the negative sign problem
in the quantum Monte Carlo method. Hence,the properties of quantum spin liquid have notbeen well understood theoretically.
On the other hand, attempts to realize quantum spin liquids other than the Heisenberg
model have been made. Among them, the Kitaev model has recently attracted considerable
attention because its ground state is a quantum spin liquid as an exact solution. Experimentally, many physical quantities such as speciﬁc heat, entropy, magnetic excitation, andthermal transport properties have been measured in ruthenium and iridium compounds,
which are believed to be described by the Kitaev model, and have been compared with theoretical studies. In the real systems, eﬀects
of the disorder are inevitably present, which
complicates extracting the nature of the Kitaev quantum spin liquid. Thus, calculations
incorporating such eﬀects are necessary to separate the properties intrinsic to Kitaev quantum spin liquids from others.

and bond randomness using quantum Monte
Carlo simulations. The spin Hamiltonian isrewritten into the free Majorana fermion system coupled to Z2 variables. The Majorana
fermion system is diagonalized using the LAPACK library, and the conﬁgurations of theZ2 variables are updated by the Markov chain
Monte Carlo method. To avoid the freezing of
the conﬁgurations, we employ the parallel tempering technique using the MPI library. Here,
we calculate thermodynamic quantities such as
the speciﬁc heat and thermal transport [1]. We
ﬁnd that the speciﬁc heat remains largely intact even in the presence of disorder at higher
temperatures compared to the spin-exchange
energy, while the low-temperature behavior issensitive to the disorder. For thermal transport, the longitudinal component of the thermal conductivity is suppressed by the disorder. This behavior is common to site dilution
and bond randomness. On the other hand, thethermal Hall conductivity is aﬀected by these
disorders diﬀerently. We ﬁnd that the quantization of the thermal Hall conductivity is robust against the bond randomness but is fragile for the site dilution.

References[1] J. Nasu and Y. Motome: Phys. Rev. B
102 (2020) 054437.

In this study, we perform numerical calculations for the Kitaev model with site dilution

209

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

Effect of quenched disorder on a triangular lattice fermionicmodel
Chisa HOTTA
Department of Basic Science,
The University of Tokyo, Meguro-ku Komaba, Tokyo 153-8902
We studied the fermionic model on amoving the location of pins, and eachtriangular lattice with quenched disorder[1]. We

timediagonalizing the ball-fermions under theare targeting the so-called tV model on ainteraction from pins, the system is safely

triangular lattice consisting of transfer integrals

thermalized to the pinball-liquid phase. Wet and nearest neighbor Coulomb repulsions V,

made full use of the ISSP supercomputing

whose ground state is known to host a Pinballsystem to perform the calculation up to theliquid phase [2,3] at large V/t, which is a statesystem size of 24x24; the finite temperature

where part of the fermions localize on one ofphase transition is observed as a divergence inthe three-sublattices, and form a symmetry

the specific heat. Also the low energy profile ofbreaking long range order. The rest of thethe pinball liquid state is well reproduced. We

fermions remain metallic and propagate along

also found that at extremely low temperature,

2/3 of the lattice sites forming a honeycomb

the introduction of the quenched disorder willband. To be precise, the long range ordered

drive the system to some sort of a glassy phase,fermions (pin) and metallic fermions(ball) are

where there starts to appear a domain wall

not fully separate, but exchange by the quantum

structure of pinball phase. This domain

fluctuation. In this study, we focus on the finite

becomes distinct as the system size becomes

temperature properties of this phase, by adding

larger and the degree of disorder is increased,

a quenched disorder to the transfer integral.

which is detected in the broad structure of the

However, it is quite difficult to fully solve

Replica overlap distribution function.

the standard tV model at finite temperature in a

References

reasonably large system size. Therefore, we[1] Keita Urakami and Chisa Hotta, in preparation.

introduce[2] Chisa Hotta and Nobuo Furukawa , Phys. Rev.

the

Falikov-Kimball

type

of

approximation and set the pins to be classical,

B 74, (2006)193107.

while keep the balls to be quantum, and[3] Chisa Hotta, Nobuo Furukawa, Akihikoperform the Monte Carlo calculation. By

Nakagawa, and Kenn Kubo , J. Phys. Soc. Jpn. 75,

changing the ratio of pin and ball and by(2006) 123704.

210

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

Ordered states in strongly correlated Dirac electron
systems of organic conductors
Akito KOBAYASHI
Department of Physics, Nagoya University
Nagoya 464-8602

In the low temperature phase (T < Tco =135K) of the organic Dirac electron system α(BEDT-TTF)2 I3 , it is believed that the electron correlation breaks the spatial inversion
symmetry and transitions to a stripe chargeordered insulator. The electron correlationinduces the reshaping of Dirac cones andthe anomalous spin fluctuation. We investigated the electron correlation eﬀects on theanomalous thermoelectric eﬀects and the spin
fluctuations[2, 4, 6]. In contrast, the elementsubstituted α-(BETS)2 I3 shows insulating behavior at low temperatures below 50 K, but nocrystal symmetry breaking or charge densitychange is observed in X-rays, and no sign ofmagnetic transition is obtained from NMR. In
this study, we constructed an extended Hubbard model based on first-principles calculations and analyzed it using mean-field approximation. It was shown that a magnetic transition occurs due to a short-range Coulomb interaction as shown in Fig. 1[3]. However, sinceno magnetic transition was observed in the experiment, it is necessary to investigate otherpossibilities. In addition, we showed possibleedge magnetism in Dirac nodal line systems of
single component molecular conductors [1, 5].

References[1] T. Kawamura, B. Zhou, Akiko Kobayashi,
Akito Kobayashi, accepted to J. Phys.
Soc. Jpn. (2021/03/30)

Figure 1: A candidate of ordered states in α(BETS)2 I3 .
[2] M. Hirata, A. Kobayashi, C. Berthier and
K. Kanoda, Rep. Prog. Phys. 84 036502(2021)[3] D. Ohki, K. Yoshimi, and A. Kobayashi,
Phys. Rev. B 102, 235116 (2020)[4] D. Ohki, M. Hirata, T. Tani, Kazushi K.,
and A. Kobayashi, Phys. Rev. Research 2,033479 (2020)[5] T. Kawamura, D. Ohki, B. Zhou, Akiko
Kobayashi, Akito Kobayashi, J. Phys.
Soc. Jpn. 89, 074704 (2020)[6] D. Ohki, Y. Omori, A. Kobayashi, Phys.
Rev. B 101, 245201 (2020)

211

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

Development of electronic structure calculation forsolids with post-Hartree-Fock calculations
Kazutaka NISHIGUCHI
Department of Physics, Osaka University
Toyonaka, Osaka 560-0043, Japan– Electronic structure calculation for solids
with post-Hartree-Fock calculations –Continuing from our preceding-term study, we
have theoretically studied a new method ofelectronic structure calculation for solids withquantum chemistry calculations. Our proposed procedure is as follows. We first startwith Hartree Fock (HF) calculation for solids
and construct maximally localized Wannierfunctions as a localized basis set. Then we
consider a real-space cluster in the Wannier
representation and solve it numerically withequation-of-motion coupled cluster (EOM-CC)
theory. Here we have used VASP (The Vienna
Ab initio Simulation Package) [1, 2, 3, 4, 5]
in the HF calculations, Wannier90 [6, 7, 8, 9]in the wannierization method, and GELLAN
program [10] in the EOM-CC calculations.
We apply this method to typical semiconductors, namely, bulk Si and SrTiO3 , and evaluate these energy band gaps as a benchmark.
Developing the preceding-term study, we use
symmetry-adapted Wannier functions [9] inconstructing the localized basis set with Wannier90. We then can observe that the localizedbasis set can be obtained stably and the band
gap can be improved compared to that with
the HF calculations. However, it is still insuﬃcient because the number of the localizedbases is not suﬃcient to take into account theelectron correlations. To increase the numberof the localized bases will be an interesting future work.

CaZn2 X2 (X = Sb, As, P) –Continuing from our preceding-term study,
we have numerically evaluated a formationenergy of n-type doping in a strong candidate of high performance thermoelectric materials CaZn2 X2 (X = Sb, As, P) by the
first-principles calculations based on the density functional theory (DFT). We consider thechemical doping into the interstitial site andthe element substitution of Ca, where thedopant is assumed to be the alkaline earth
metals (Mg, Ca, Sr, Ba) and group 3 elements (Sc, Y, La). The ionization energy ofthese elements is relatively small as cations
so that the smaller formation energy is expected. To take into account the chemical doping, the 3 × 3 × 2 supercell calculations are
performed (Figure 1) with the Perdew Burke
Ernzerhof (PBE) functional, namely, the generalized gradient approximation (GGA). Here
we have used VASP [1, 2, 3, 4, 5] in the firstprinciples calculations.

Figure 1: 3 × 3 × 2 supercell of CaZn2 X2 (X =Sb, As, P). (Left) Doping into the interstitialsite. (Right) Element substitution of Ca.

– Formation energy of n-type doping in

For the interstitial doping one can observe212

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

in Figure 2 that the chemical doping of Ca
and group 3 elements is favorable according
to the formation energy. In addition, for theelement substitution of Zn one can also see in
Figure 3 that the element substitution for Sc
is stable. This tendency of the formation energy is considered to depend on the ionizationenergy of these doping elements and the ionicradius against the interstitial site and Ca vacancy site. The improvement of the exchange
correlation functional such as the hybrid functional and the inclusion of the finite-size supercell eﬀect will be important future works.

Figure 3: Formation energy of the element substitution of Ca.
[3] G. Kresse and J. Furthmüller, Comput.
Mater. Sci. 6 15 (1996).
[4] G. Kresse and J. Furthmüller, Phys Rev.
B 54 11169 (1996).
[5] VASP (The Vienna Ab initio Simulation
Package): https://www.vasp.at[6] N. Marzani and D. Vanderbilt, Phys. Rev.
B 56 12847 (1997).
[7] I. Souza, N. Marzani, and D. Vanderbilt,
Phys. Rev. B 65 035109 (2001).
Figure 2: Formation energy of the doping into
the interstitial site.

[8] Wannier90: http://www.wannier.org[9] R. Sakuma, Phys. Rev. B 87 235109(2013).

References[1] G. Kresse and J. Hafner, Phys. Rev. B 47
558(R) (1993).

[10] GELLAN: A hierarchical massivelyparallel quantum chemistry program,
Kobe University[2] G. Kresse and J. Hafner, Phys. Rev. B 49
14251 (1994).

213

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

Development of electronic structure calculation forsolids with quantum chemistry calculations
Kazutaka NISHIGUCHI
Department of Physics, Osaka University
Toyonaka, Osaka 560-0043, Japan– Electronic structure calculation for solids
with quantum chemistry calculations –A new method of electronic structure calculation for solids with quantum chemistry calculations have been studied theoretically. We firststart with Hartree Fock (HF) calculation forsolids and construct maximally localized Wannier functions as a localized basis set. Then
we consider a real-space cluster in the Wannier
representation and solve it numerically withequation-of-motion coupled cluster (EOM-CC)
theory. Here we have used VASP (The Vienna
Ab initio Simulation Package) [1, 2, 3, 4, 5]
in the HF calculations, Wannier90 [6, 7, 8] inthe wannierization method, and GELLAN program [9] in the EOM-CC calculations.
We apply this method to typical semiconductors, namely, bulk Si and SrTiO3 , and evaluate these energy band gaps as a benchmark.
We can observe that the band gap can be improved compared to that with the HF calculations, but it is not suﬃcient because the localized basis set is diﬃcult to converge stably.
To achieve stable construction of the localizedbasis set will be an interesting future work.
– Formation energy of n-type doping in
CaZn2 Sb2 –A formation energy of n-type doping in astrong candidate of high performance thermoelectric materials CaZn2 Sb2 have been numerically evaluated by the first-principles calculations based on the density functional theory(DFT). We consider the chemical doping into
the interstitial site and the element substitu-

tion of Zn, where the dopant is assumed to be
the alkaline earth metals (Mg, Ca, Sr, Ba) andgroup 3 elements (Sc, Y, La). The ionizationenergy of these elements is relatively small as
cations, so that the smaller formation energy
is expected. To take into account the chemical
doping, the 3 × 3 × 2 supercell calculations are
performed (Figure 1) with the Perdew Burke
Ernzerhof (PBE) functional, namely, the generalized gradient approximation (GGA). Here
we have used VASP [1, 2, 3, 4, 5] in the firstprinciples calculations.

Figure 1: 3 × 3 × 2 supercell of CaZn2 Sb2 .
(Left) Doping into the interstitial site. (Right)
Element substitution of Zn.
For the interstitial doping one can observe in
Figure 2 that the chemical doping of Ca andgroup 3 elements is favorable according to theformation energy. In addition, for the element
substitution of Zn one can also see in Figure
3 that the element substitution for Mg, Ca,Sr, and group 3 elements is stable. This tendency of the formation energy is considered todepend on the ionization energy of these doping elements and the ionic radius against theinterstitial site and Zn vacancy site. On theother hand, the calculated formation energy is214

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

now understood to be qualitative because thePBE functional is not suﬃcient to evaluate theaccurate energy band gap for semiconductors
and the finite-size supercell eﬀect is not included in these calculations, but the tendencyof the formation energy is important for material synthesis of n-type semiconductors. The
improvement of the exchange correlation functional such as the hybrid functional and theinclusion of the finite-size supercell eﬀect will
be important future works.

[4] G. Kresse and J. Furthmüller, Phys Rev.
B 54 11169 (1996).
[5] VASP (The Vienna Ab initio Simulation
Package): https://www.vasp.at[6] N. Marzani and D. Vanderbilt, Phys. Rev.
B 56 12847 (1997).
[7] I. Souza, N. Marzani, and D. Vanderbilt,
Phys. Rev. B 65 035109 (2001).
[8] Wannier90: http://www.wannier.org[9] GELLAN: A hierarchical massivelyparallel quantum chemistry program,
Kobe University

Figure 2: Formation energy of the doping into
the interstitial site.

Figure 3: Formation energy of the element substitution of Zn.

References[1] G. Kresse and J. Hafner, Phys. Rev. B 47
558(R) (1993).
[2] G. Kresse and J. Hafner, Phys. Rev. B 49
14251 (1994).
[3] G. Kresse and J. Furthmüller, Comput.
Mater. Sci. 6 15 (1996).

215

Activity Report 2020 / 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

the band parameters.
0.006

PS

∆E /N

We carried out numerical calculations on
the basis of the optimized variational Monte
Carlo method[1, 2, 3, 4]. We have investigated the ground-state phase diagram of thetwo-dimensional Hubbard model and the twodimensional d-p model. We performed parallel computations by using a Monte Carlo algorithm. In order to reduce statistical errors,
we carried out 200 ∼ 500 parallel calculations.
Parallel computing is very necessary to reduce
Monte Carlo statistical errors.

0.004

AF
0.002

AFI
SC
0

We employed the improved wave function ofan exp(−λK) − PG -type[1]. This wave function is a very good many-body wave function because the ground-state energy is lowered
greatly and the ground-state energy is lower
than those that are evaluated by any otherwave functions[2]. We can improve the wave
function systematically by multiplying by operators PG and e−λK several times.
We exhibit the condensation energy as afunction of the doping rate x for the three-band
d-p model in Fig. 1 where Ud = 10tdp and thelevel difference is taken as ∆dp ≡ ǫp − ǫd =tdp [6]. There is the AF region when x is small
and the SC region exists near the optimum region x ∼ 0.2. We should mention that there
is a phase-separated region in the low-doping
region where x < 0.07. This indicates the existence of AF insulator phase in the low doping
region. This is similar to the phase diagramof the 2D Hubbard model. The phase separation (PS) is, however, dependent on the level
difference ∆dp . As ∆dp decreases, the phaseseparated region decreases and vanishes when
∆dp approaches zero. Thus the area of phaseseparated region can be controlled by changing

216

0

0.05

0.1

0.15

0.2

0.25

0.3

x

Figure 1: AF and SC condensation energies asa function of the hole doping rate x for ∆dp = 1
on a 8 × 8 lattice with tpp = 0.4 and Ud = 10 inunits of tdp . There is a phase-separated region
when x is small.

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).
[5] T. Yanagisawa, Phys. Lett. A403, 127382(2021).
[6] T. Yanagisawa et al., EPL (2021).

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

Calculation of the electronic state in a large model of theorganic charge ordering system using mVMC
Masatoshi SAKAI, Riku TAKEDA, Ryoma ISHII
Department of Electrical and Electronic Engineering,
Chiba University, 1-33 Yayoi-cho Inage, Chiba 263-8522
We have studied the electronic state of thethe system.

large-size model of the organic charge order

The other work in progress is about -(BEDT-

system by using the supercomputer. We made aTTF)2PF6, of which metal-insulator phasebase model of -(BEDT-TTF)2I3 by citingtransition is 297 K [2]. In this work, we could

reference [1] and carried out a 28 x 28 site

not find the physical parameters in previous

calculation to reveal the electron number

works, so we began to calculate from the band

dependence of the charge order structure of thestructure before carrying out model calculation

correlated electron system. The simulationof correlated electronic system. We usedsoftware mainly used in this work was mVMC,

Quantum Espresso and Respack to obtainas preinstalled in the ENAGA and OHTAKA

physical parameters to build the calculation

supercomputer systems. The electron number

model used in mVMC or HPhi. We learned

dependence of the charge order structure is

Quantum Espresso and Respack, and now wesignificant for us to support our experimental

have just achieved the end of the calculation ofwork. Our experimental work was the metal-

the Respack. Therefore, we will verify theinsulator phase transition transistor which was

results and apply the next step by using mVMC

observed in organic charge order material of -or HPhi to reveal the physical principle of the(BEDT-TTF)2I3 near the metal-insulator phase

phase transition field-effect transistor.

transition at 137 K. We expect that thecalculation of electron number dependence ofReferencescharge order structure and energy will give a[1] Daigo Ohki, Genki Matsuno, Yukiko Omori

theoretical backbone to our experimental result.

and Akito Kobayashi: Crystals 8, (2018) 137.

Our group made an effort to learn about[2] Hayao Kobayashi, Takehiko Mori, Reizo

supercomputer and mVMC, and now we are

Kato, Akiko Kobayashi, Yukiyoshi Sasaki,

finally in routine calculation by varying

Gunzi Saito and Hiroo Inokuchi: Chem. Lett.

electron number from 376 to 404 electrons in(1983) 581.

217

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

Study for high-Tc cuprates using four-band d-p model
Hiroshi WATANABE
Research Organization of Science and Technology, Ritsumeikan University
1-1-1 Noji-Higashi, Kusatsu-shi, Shiga 525-8577

The discovery of superconductivity incuprates [1] has brought about the significant
progress in strongly-correlated electron systems. However, we have not yet reached a unified understanding of their properties, including their material dependence of the superconducting transition temperature Tc . Recent development of theories and experiments sheds
light on the importance of the orbital degreeof freedom in cuprates [2, 3]. To elucidate thematerial dependence of Tc , we study the fourband d-p model for La2 CuO4 and HgBa2 CuO4
systems that properly includes the orbital degree of freedom. The ground state properties
are studied with the variational Monte Carlo(VMC) method. The Gutzwiller-Jastrow typewave function is used for the VMC trial wave
function. The system size for the calculation is N =24×24=576 unit cells (and thus
576×4=2304 orbitals in total), which is large
enough to avoid finite size eﬀects.
The superconducting correlation function
vs hole doping rate δ for the La2 CuO4
and HgBa2 CuO4 systems are shown in Fig. 1.
The dome-shaped behavior and the material
dependence are consistent with experiments.
We show that the Cu dz 2 orbital contribution
around the Fermi energy is destructive for dwave superconductivity, and thereby becomesa key factor to determine Tc . The energy diﬀerence ∆dp between Cu dx2 −y2 and O p orbitals
is also shown to be a key factor. Our result
accounts for the empirical correlation between
Tc and model parameters and gives the unified
description beyond the usual one-band Hubbard, t-J, and even three-band d-p models.
P dd

Figure 1: Superconducting correlation function P dd vs hole doping rate δ for the La2 CuO4
and HgBa2 CuO4 systems.

References[1] J. G. Bednorz and K. A. Muller, Z. Phys.
B 64, 189 (1986).
[2] H. Sakakibara et al., Phys. Rev. Lett. 105,
057003 (2010).
[3] C. E. Matt et al., Nat. Commun. 9, 972(2018).

218

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

Studies of the superconductivity and magnetic statesin the strongly correlated electron systems using
Hubbard models.
Atsushi Yamada
Department of Physics, Chiba University
Chiba 263-8522, Japana, Chiba 277-8581

A spin liquid state has attracted a lotof interest. This state is realized in geometrically frustrated systems like the chargeorganic transfer salts κ-(BEDT-TTF)2 X[1]and Cs2 CuCl4 .[2] Hubbard model on the anisotropic triangular lattice is a simple theoretical model of these compounds, and spin liquid state is found in this model.[3] A spin liquid could arise also in the intermediate coupling region of strongly correlated systems between a semi-metal and ordered state, because
in this case a 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, where a semi-metal is realized at
U = 0.
We have studied the magnetic and metal-toinsulator transitions by variational cluster approximation using 10-site and 16-site clusters
as a reference system. Parts of numerical calculations were done using the computer facilities of the ISSP. We found that UAF = 2.7 and
UMI = 3.0 for 10-site cluster, and UAF = 2.7
and UMI = 3.2 for 16-site cluster.[5] This result
also rules out the existence of the spin liquid inthis model. Our results agree with recent large
scale Quantum Monte Carlo simulations.[4]We are currently improving our program using MPI technique so that we will be able tostudy a larger cluster size system.

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).

219