3.4

Cooperative Phenomena

in Complex Macroscopic Systems

238

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

Tensor-network-based methods and their application todiscrete lattice models
Naoki KAWASHIMA
Institute for Solid State Physics,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
Tensor network (TN) is a useful concept in

the local-minima-trapping due to the formation

various fields, not only in the condensed matteroftheory, but also in data science and machine

benchmark showed that the new method

learning. Even in the condensed matter theory

efficiently removes such loops when the target

alone, the application of the TN concept is

tensor has, at least approximately, the structure

multi-fold. In variational principles calculation,of the corner double line tensor, which is

the TN is used as a new type of variational

typical in TN-based RG calculations.

wave

function.

Inthereal-space

redundantIn[2,3],

entanglement

wedeveloped

loops.

theOur

tensorrenormalization group (RG) theory, it serves asrenormalization group (TRG) algorithm fora natural representation of a renormalized

statistical systems with open boundaries. Thissystem. In the pure mathematical physics, it isallows us to investigate properties of boundary

used as a new language for discussing theoperators. We applied the method to several

topological nature of the quantum states. In

two-dimensional classical spin models with

SY2019, we developed new algorithms foropen boundaries. The boundary conformal

generalthespectrumwasobtainedboundary conformal-field theory [2,3], theconsistent

with

the

classical

conformal field theories (BCFTs) for the Ising,TN-based
dynamics

computation
[4],

and

[1],

variational

optimization of topological quantum states.

andturnedunderlyingoutboundary

tricritical Ising, and three-state Potts models.

[5,6]

In [4], we developed a method for tracing

In [1], we proposed a new method forthetime-evolution

of

the

probability

obtaining an efficient representation of a many-

distribution function of a given dynamicalranked tensor in terms of a ring network ofsystems, which is otherwise difficult. We

smaller-ranked tensors. This kind of operations

studied (1+1)-dimensional directed percolation,is vital both in TN application to the real spacea canonical model of a nonequilibrium

renormalization group and in data science. The

continuous transition. We find a universal

central technical problem here is, how to avoid

relaxation of Rényi entropy at the absorbing

239

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

transition point. We also detected a singularity

phase appears between the low field zig-zag

in the "entanglement entropy" in the active (i.e.,magnetic order and the high field polarized

non-absorbing) phase. Although the absorbing

state. We found that the chiral Kitaev spinstates is exponentially rare in this phase, its rareliquid occupies an area in the phase diagram,appearance causes this singularity.

smaller than predicted by other groups.

In [5], we studied the Kitaev model by the
TN representation. We proposed a new wave

References

function that is extremely simple in the TN[1] Hyun-Yong Lee and Naoki Kawashima: J.

representation but still shares many essential

Phys. Soc. Jpn. 89, 054003 (2020)

properties with the ground state of the Kitaev[2] Shumpei Iino, Satoshi Morita, and Naoki

model, such as the quasi-long-range correlation,

Kawashima: Phys. Rev. B 101, 155418 (2020)

the short-range magnetic correlation, the flux-[3] Shumpei Iino, Satoshi Morita, and Naoki

free property, the Z2 gauge structure, and the

Kawashima: Phys. Rev. B 100, 035449 (2019)2D Ising universality. We call this state the[4] Kenji Harada and Naoki Kawashima: Phys.

"loop gas state" because of its mathematical

Rev. Lett. 123, 090601 (2019)equivalence to the loop gas model, one of the[5] Hyun-Yong Lee, Ryui Kaneko, Tsuyoshi

standard and solvable statistical physicalOkubo, and Naoki Kawashima: Phys. Rev. Lett.

models. We showed that the loop gas state is123, 087203 (2019)adiabatically connected to the ground state of[6] Hyun-Yong Lee, Naoki Kawashima, andthe Kitaev model, indicating a similar role to

Yong Baek Kim:arXiv:1911.07714

the AKLT state for the S=1 AFH chain.

[7] Hyun-Yong Lee, Ryui Kaneko, Tsuyoshi
Okubo, Naoki Kawashima: Phys. Rev. B 101,

Later we extended our method to several
other models: S=1 Kitaev model [6], and035140 (2020)

Kitaev-Gamma model [7,8]. In particular, in[8] Hyun-Yong Lee, Ryui Kaneko, Li Ern[8], we studied field-induced quantum phases

Chern, Tsuyoshi Okubo, Youhei Yamaji, Naoki

in theoretical models of the Kitaev magnets,

Kawashima and Yong Baek Kim: Naturesuch as α-RuCl3 for which the half quantizedCommunications 11, 1639 (2020)

thermal Hall conductivity was discovered. Our
TN calculation suggested that a field-induced240

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

Kármán-Vortex Cavitation Around Circular-Cylinder Arrays
Yuta ASANO1, Hiroshi WATANABE1, 2, and Hiroshi NOGUCHI1
1

Institute for Solid State Physics, University of Tokyo,
Kashiwa-no-ha, Kashiwa, Chiba 277-8581, Japan2

Department of Applied Physics and Physico-Informatics, Keio University,
Yokohama, Kanagawa 223-8522, Japan

Cavitation

is

a

flow

phenomenon

bubbles and their effects on the surrounding

accompanied by the formation, growth, andflow are investigated by MD calculations [2].

collapse of bubbles caused by local pressure

3-dimensional MD simulation is carried out tochanges in a liquid with the high flow velocity.

analyze the flow around the periodicallySince cavitation has adverse effects on fluidaligned cylinders. The Lennard-Jones (LJ)machinery such as performance degradation,interaction is used for the interaction between

vibration and noise, and erosion, the elucidation

the fluid particles. The cylinder is modeled byof its mechanism is of great engineering

fixing LJ particles on its surface. Periodic

importance. The cavitation inception is thoughtboundary conditions are used for all directions

to be caused by the growth of small bubble

in the system, and two cylinders are arranged

nuclei in a low-pressure region for enough time

perpendicular to the flow to represent an

to grow [1]. However, since it is extremely

infinite row of cylinders. In this way,difficult to discuss the dynamics of microscopic

interference of the vortices behind the upper

bubble nuclei in a macroscopic flow field, theand lower cylinders can be taken into account.

mechanism of cavitation generation is not well

The temperature and velocity are controlled byunderstood. In addition, because the phase

the Langevin heat bath at the downstream sidetransition has to be directly discussed in theof the cylinder. Fluids in the supercriticalflow field, it is extremely difficult to analyze

region are inflowed. Cavitation is generated by

the cavitation inception of the liquid in whicha decrease in the local temperature of the

no bubble nuclei are contained as in the case of

inflow fluid.

flow in hydraulic machines.

Figure 1 shows a snapshot of the density (left

In the present study, molecular dynamics

panels) and vorticity (right panels) fields. At(MD) calculations are performed to analyze

high temperature (temperature T=2), there arecavitation around periodically aligned cylinders,no bubbles and therefore non-cavitating flow isand the formation and growth processes ofobtained. An anti-phase synchronized Kármán

241

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

vortex is generated behind the cylinder. It ischanges occur due to temperature changes.

also found that this synchronization amplifies

The change in the vortex structure should

the vibrations acting on the cylinder. At thechange the characteristics of the vibrations

lower

are

excited by the vortex. It is found that thegenerated in conjunction with the shedding

bubbles generated near the cylinder inhibit thecycle of the Kármán vortex. The Kármán

propagation of the vibration associated with thevortex remains synchronized in the anti-phase

ejection of the vortex, and the vibration

as in the non-cavitating flow. At T=1.25, a gas-

amplitude decreases and eventually disappears.

phase region is formed behind the cylinder andIn summary, the analysis from the molecular

the distance to the formation of the Kármán

scale reveals that bubbles generated near thevortices increases. The phase of the Kármán

cylinder significantly change the properties ofvorticesthe lift and flow fields.

temperature,isslightly

T=1.3,

bubbles

disturbed.

As

the

temperature is further reduced to T=1.2, thegas-phase region behind the cylinder is further

References

expanded and the upper and lower Kármán

[1] C. E. Brennen, Cavitation and Bubble

vortices become asymmetric. It is also found

Dynamics (Oxford University Press, 1995).

that this asymmetric structure is switched by a[2] Y. Asano, H. Watanabe, and H. Noguchi: J.

long period of time. In Newtonian fluids in

Chem. Phys. 152 (2020) 034501.

which no phase transitions occur, no such

Figure 1: Snapshots of density (left panels) and vorticity (right panels)
fields of LJ fluid. The black circles in the figure represent cylinders.

242

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

Detachment of fluid membrane from substrate byprotein adhesion
Hiroshi NOGUCHI
Institute for Solid State Physics, University of Tokyo, Kashiwa, Chiba 277-8581

Lipid membranes supported on a solid substrate are considered as a model system for biological membranes and are extensively used to
study protein functions and membrane properties. Recently, the membrane detachmentfrom the substrate was reported by the adhesion of Ca2+ -dependent membrane-binding
proteins annexin. The membrane rolling andblebbing have been observed. To understandthe mechanism to form these membrane structures, we have perfromed meshless membrane
simulations with the isotropic spontaneous
curvatures. We focus on the generation of the
spontaneous curvature by the protein adhesion
here.
Before considering the detachment from thesubstrate, we investigated the closing transition from a membrane patch to a vesicle at a
finite spontaneous curvature [1]. We have determined the theoretical condition of the closure and proposed the estimation method of
the membrane Gaussian modulus from the curvature of the membrane patch.
Next, we added the attraction to a solid substrate by the Lennard-Jones potential. When ahigh curvature C0 is imposed, a membrane disk
is detached from the substrate, and vesicles areformed [see Fig. 1] [2]. The membrane edge
is rolled up, and the edge undulation expands
leading to the vesicle formation. The averagevesicle size is several times larger than the minimum vesicle size formed by the membrane closure of a cup-shaped membrane patch.
With strong adhesion, the competition between the bending and adhesion energies determines the minimum value of the spontaneous curvature for the detachment. In contrast, with weak adhesion, a detachment occurs at smaller spontaneous curvatures due

243

Figure 1: Sequential snapshots of the detachment and vesiculation of the disk-shaped membrane from a solid substrate.
to the membrane thermal undulation. Whenparts of the membrane are pinned on the substrate, the detachment becomes slower, and aremained membrane patch forms straight orconcave membrane edges. The edge undulation induces vesiculation of long strips anddisk-shaped patches. Therefore, membrane
rolling is obtained only for membrane strips
shorter than the wavelength for deformation
into unduloid. This suggests that the rollingobserved for the annexins A3, A4, A5, andA13, results from by the anisotropic spontaneous curvature induced by the proteins.

References[1] H. Noguchi: J. Chem. Phys. 51, 094903(2019).
[2] H. Noguchi: Soft Matter 15, 8741 (2019).

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

Brittle structure of semicrystalline polymers bylarge-scale coarse-grained molecular dynamics
simulation
Yuji HIGUCHI
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581

We have studied mechanical properties ofsemicrystalline polymers [1, 2, 3] and doublenetwork gels [4] via coarse-grained molecular dynamics simulations. Understanding thefracture process of the semicrystalline polymers at the molecular level is essential to increase its toughness against the stretching.
This year, void formation and growth processes in the fracture process of semicrystalline
polymers were studied on System B.
Figure 1 shows the fracture process of
the lamellar structure consisting of 1.5 × 107
monomers. The lamellar structure consists ofamorphous and crystalline layers. Voids are
generated in amorphous layers, and then theygrow parallel to the stretching direction cylindrically. This process is diﬀerent from thosein elastomers and glassy polymers. The void
structures during their growth process at themolecular level are still debatable in experiments. The maximum length of voids along
x, y, and z-directions are monitored during
the stretching. At ﬁrst, the maximum lengthof voids in x and y-directions, which are perpendicular to the stretching direction, is larger
than that in z-direction, which is parallel tothe stretching direction. Then, the length inz-direction becomes larger than those in x andy-directions. This result is consistent withthat observed in the experiment [5]. Thus,
the void generation and growth processes in
the lamellar structure of crystalline polymers
at the molecular level are successfully revealed.

Snapshot
Amorphous
CrystalStrain = 0.00
Strain = 0.75

Void growth

Strain = 1.00

Strain = 1.50

Void

Strain = 0.00
Strain = 0.75

Strain = 1.00 Strain = 1.50

Figure 1: The deformation and fracture process of the lamellar structure. (a) Snapshotand (b) void generation and growth processes.

References[1] Y. Higuchi and M. Kubo: Macromolecules
50 (2017) 3690.
[2] Y. Higuchi: Polym. J. 50 (2018) 579.
[3] Y. Higuchi: Macromolecules 52 (2019)6201.
[4] Y. Higuchi, K. Saito, T. Sakai, J. P. Gong,
and M. Kubo: Macromolecules 51 (2018)3075.
[5] S. Humbert, O. Lame, J. M. Chenal, C.
Rochas, and G. Vigier: Macromolecules
43 (2010) 7212.

244

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

Topological Magnetization Plateau of Quantum
Antiferromagnetic Spin Chain with Anisotropies
Tôru SAKAI1,2 , Hiroki NAKANO1,2 , and Tokuro SHIMOKAWA3 ,1
Graduate School of Material Science, University of Hyogo,
Kouto, Kamigori, Hyogo 678-1297, Japan
2
Synchrotron Radiation Research Center, Kansai Photon Science Institute,
Quantum Beam Science Research Directorate,
National Institute for Quantum and Radiological Science and Technology (QST)
SPring-8, Kouto, Sayo, Hyogo 679-5148, Japan
3
Okinawa Institute of Science and Technology Graduate University, Onna,
Okinawa 904-0945, JapanThe symmetry protected topologicalphase[1, 2] is one of interesting topics inthe strongly correlated electron systems.
The recent study[3, 4] using the numerical
diagolaization of finite-size clusters and thelevel spectroscopy analysis revealed that thesymmetry protected topological phase appears
in the S = 2 antiferromagnetic chain withthe single-ion anisotropy D and the coupling
anisotropy ∆. However, the region of thisphase is quite small in the ground state phase
diagram. Thus it would be diﬃcult to discoverin some real experiments. In the present
study we focused on the symmetry protected topological phase in the magnetization
plateau of the S = 2 antiferromagnetic chain,using the same numerical diagonalizationand level spectroscopy analysis. As a result,
it is successfuly found that the symmetry
protected topological phase appears in muchlarger region of the phase diagram of the
1/2 magnetization plateau. In addition thebiquadratic interaction is revealed to enhanse
this phase[5].
The 1/3 magnetization plateau of the distorted diamond spin chain was also investigated using the nimerical diagonalization andthe finite-size scaling analysis. As a result it

245

is found that the magnetization plateau vanishes for suﬃciently large XY -like coupling
anisotropy. The phase diagram of the 1/3 magnetization plateau was also obtained[6, 7].
As related topics to the magnetization
plateau, we invesitgated the spin gap of
the kagome-lattice antiferromagnet. We performed the large-scale numerical diagonalization of the spin gap and the magnetization
curve of the system and tried to compare themwith the recent experimental results[8, 9]References[1] F. Pollmann, A. M. Turner, E. Berg and
M. Oshikawa, Phys. Rev. B 81, 064439(2010).
[2] F. Pollmann, E. Berg, A. M. Turner and
M. Oshikawa, Phys. Rev. B 85, 075125(2012).
[3] T. Tonegawa, K. Okamoto, H. Nakano, T.
Sakai, K. Nomura and M. Kaburagi, J.
Phys. Soc. Jpn. 80, 043001 (2011).
[4] K. Okamoto, T. Tonegawa and T. Sakai,
J. Phys. Soc. Jpn. 85, 063704 (2016).

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[5] T. Sakai, K. Okamoto and T. Tonegawa,Phys. Rev. B 100, 054407 (2019).
[6] Y. Ueno, T. Zenda, Y. Tachibana, K.
Okamoto and T. Sakai, JPS Conference
Proceedings 30, 011085 (2020).
[7] T. Zenda, Y. Tachibana, Y. Ueno, K.
Okamoto and T. Sakai, JPS Conference
Proceedings 30, 011084 (2020).
[8] H. O. Jeschke, H. Nakano and T. Sakai,Phys. Rev. B 99, 140410(R) (2019).
[9] R. Shirakami, H. Ueda, H. O. Jeschke,
H. Nakano, S. Kobayashi, A. Matsuo, T.
Sakai, N. Katayama, H. Sawa, K. Kindo,
C. Michioka and K. Yoshimura, Phys.
Rev. B 100, 8722674 (2019).

246

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

Numerical Diagonalization Study on Spin Nematic
Phase of Low-Dimensional SystemTôru SAKAI1,2 , Hiroki NAKANO1,2 , and Tokuro SHIMOKAWA3 ,1
Graduate School of Material Science, University of Hyogo,
Kouto, Kamigori, Hyogo 678-1297, Japan
2
Synchrotron Radiation Research Center, Kansai Photon Science Institute,
Quantum Beam Science Research Directorate,
National Institute for Quantum and Radiological Science and Technology (QST)
SPring-8, Kouto, Sayo, Hyogo 679-5148, Japan
3
Okinawa Institute of Science and Technology Graduate University, Onna,
Okinawa 904-0945, JapanThe spin nematic phase has attracted a lotof interests in the field of the quantum spin systems. It was theoretically predicted to appear
around the saturation magnetization in themagnetization process of the frustrated spin
chain[1]. Our present study using the numerical diagonalization and the finite-size scaling
analysis revaled that the spin nematic liquid
phase apears in the S = 1 spin ladder systemwith the easy-axis single-ion anisotropy D evenfor a small magnetic field H. The phase diagram in the DH plane and the DH plane wereobtained[2, 3].
The research of the three-leg spin tube withthe ring-exchange interaction J4 using the unmerical diagonalization and the finite-size scaling analysis[4] indicated that the spin nematic
liquid pahse appears for suﬃciently large J4 .
In addition the spin nematic liquid phase wasalso revealed to appear in the spin ladder system with the antiferromagnetic leg and the ferromagnetic rung exchange interaction in thepresent of the XXZ anisotropy[5].

References[1] T. Hikihara, L. Kecke, T. Momoi and A.
Furusaki, Phys. Rev. B 78, 144404 (2008).

247[2] R. Yamanokuchi, K. Okamoto and T.
Sakai, Proceedings of 2018 16th International Conference on Megagauss Magnetic Field Generation and Related Topics (MEGAGAUSS) (2019) 8722674.
[3] T. Sakai and K. Okamoto, JPS Conference Proceedings 30, 011083 (2020).
[4] T. Sakai, JPS Conference Proceedings 29,
014004 (2020).
[5] K. Okamoto, T. Tonegawa, M. Kaburagi
and T. Sakai, JPS Conference Proceedings 30, 011024 (2020).

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

Aggregation process of full-length amyloid-beta peptides
studied by molecular dynamics simulations
Hisashi Okumura
Exploratory Research Center on Life and Living Systems,
Institute for Molecular Science, Okazaki, Aichi 444-8585
Alzheimer's disease is thought to be caused

number of combinations between replicas andby aggregates of amyloid-beta peptides, such asparameters without the loss of sampling

oligomers and amyloid fibrils. Amyloid-beta

efficiency. For comparison, we applied the

peptide is a peptide consisting of 39-43 amino

replica sub-permutation, replica-permutation,acid residues. The oligomer structure of theand replica-exchange methods to a β-hairpin

amyloid-betaWemini protein, chignolin, in explicit water. We

investigated the aggregation process of full-calculated the transition ratio and number oflength amyloid-beta peptides by moleculartunneling events in the parameter space, thedynamics simulations. Several amyloid-beta

number of folding–unfolding events, and the

peptides were placed in explicit water, and theirautocorrelation function to compare sampling

dynamics and structural changes were observed.

efficiency. The results indicate that among the

As a result, aggregates with intermolecular

three methods, the RSPM is the most efficientbeta-sheet structure was generated. During thisin both parameter and conformational spaces.

peptidesis

notclear.

oligomerization process, Aβ peptides are first
aggregated by electrostatic force. Then, they

References

are aggregated by hydrophobic interaction.

[1] K. Hukushima, K. Nemoto, J. Phys. Soc.
Jpn. 65 (1996) 1604.

We also developed an improvement of thereplica-exchange [1, 2] and replica-permutation[2] Y. Sugita, Y. Okamoto, Chem. Phys. Lett.

methods [3], which we call the replica sub-314 (1999) 141.

permutation method (RSPM) [4]. Instead of[3] S. G. Itoh, H. Okumura, J. Chem. Theor.

considering all permutations, this method uses

Comput. 9 (2013) 570.

a new algorithm referred to as sub-permutation[4] M. Yamauchi and H. Okumura: J. Comput.

that contains transitions only to the adjacentChem. B 40 (2019) 2694.

temperatures. The RSPM can reduce the248Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Monte Carlo study of the critical properties of noncollinear
Heisenberg magnets: O(3)×O(2) universality class
Yoshihiro Nagano, Kazuki Uematsu and Hikaru Kawamura
Graduate School of Science, Osaka University, Toyonaka 560-0043
The

criticalthetype, characterized by the complex-valued

antiferromagnetic Heisenberg model on thecorrection-to-scaling exponent, ω=0.1+0.4-0.05 +three-dimensional stacked-triangular lattice [1]i 0.7+0.1-0.4. The focus-like nature of the chiral

are studied by means of a large-scale Montefixed point accompanied by the spiral-like

Carlo simulation in order to get insight into therenormalization-group (RG) flow [5] is likely tocontroversial issue of the criticality of thebe the origin of the apparently complicatednoncollinear magnets with the O(3)×O(2)critical behavior. The results are compared andsymmetry [2]. The maximum size studied is

discussed in conjunction with the results of other3843, considerably larger than the sizes studied

numerical simulations, several distinct types ofby the previous numerical works on the model.

RG calculations including the higher-order

Availability of such large-size data enables us toperturbative

examineproperties

calculations and the nonperturbative functional

including the effect of corrections to the leading

RG calculation, and the conformal-bootstrap

scaling. Strong numerical evidence of the

program.

theproperties

detailed

criticalof

massiveandmassless

RG

continuous nature of the transition is obtained.
Our data indicates the existence of significant

References

corrections to the leading scaling. Careful[1] H. Kawamura, J. Phys. Condens. Matter 10,

analysis by taking account of the possible

4707 (1998).

corrections yield critical exponents estimates, α[2] Y. Nagano, K. Uematsu and H. Kawamura,=0.44(3),β=0.26(2),γ=1.03(5),ν=Phys. Rev. B 83, 103704 (2019).

0.52(1),

η=0.02(5), and the chirality exponents βk =[3] H. Kawamura, Phys. Rev. B 8, 4916 (1988).

0.40(3)

the[4] A. Pelissetto, P. Rossi and E, Vicari, Phys.

existence of the O(3) chiral (or O(3)×O(2))

Rev. B 63, 140414(R) (2001); B 65, 020403(R)universality class governed by a new `chiral'(2001).hi

fixed point [3,4]. We also obtain an indication[5] P. Calabrese, P. Parruccini and A.I. Sokolov,

that the underlying fixed point is of the focus-

Phys. Rev. B 66, 180403(R) (2002).

andγk=0.77(6),supporting249

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

Multiscale analysis of polymeric fluids and solids
Takahiro MURASHIMA
Department of Physics, Graduate School of Science,
Tohoku University, Aramaki-Aza-Aoba, Aoba-Ward, Sendai, Miyagi 980-8578
100 tetrahedral elements and the deformationof each element is calculated using the finite
element method. Each element has a molecularsimulation system and molecular simulation iscarried out according to the deformation. The
stresses are then obtained from the molecularsimulation system and passed on to themacroscopic system. On the macro-scale side,
the deformation occurs according to thedetermined stress field. The entire system isdeveloped in time by alternating these macroand micro-scale calculations. Multiscale
simulations are now possible to include
molecular features that are difficult to handle inconventional macro-scale simulations, and to
calculate large deformations that are difficult tohandle in conventional micro-scale simulations.

Polymeric materials are used in many
products and are now an integral part of ourlives. In recent years, however, the use of
microplastics has tended to be discouraged due
to concerns about the damage they can cause to
ecosystems. The development of new materials
with less impact on ecosystems, such as
biodegradable plastics, is expected.
In order to develop more functional
materials, it is essential to control the physicalproperties on the molecular scale and to control
the formability on the macro scale. The
micro-scale properties of polymer materials andmacro-scale mechanical properties are strongly
correlated and cannot be treated independently,
making their understanding difficult. Therefore,
we have been developing multi-scale
simulations that simultaneously combine
micro-scale
molecular simulations andmacro-scale continuum simulations. In thisproject, we have developed a multi-scale
simulation that combines molecular dynamics
simulation and finite element simulation, andthe results have been published in several
papers[1-3].
The method of multi-scale simulation is asfollows. The macro-scale cube is divided into

References[1] T. Murashima, S. Urata, S. Li: Euro. Phys.
J. B ​92​(2019) 221.
[2] T. Murashima: CICSJ Bulletin ​37 (2019)
87.
[3] S. Li, S. Urata, T. Murashima: IACM
Express. 4​6​(2020) 10.

250

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

Dynamical simulation of spin-orbit insulators
Hidemaro Suwa
Department of Physics, University of Tokyo
7-3-1 Hongo Bunkyo, Tokyo 113-0033

Dynamics of correlated quantum materials
is a long-standing problem in condensed matter physics. The recent development of experimental techniques such as inelastic neutronand x-ray scattering has revealed dynamical
properties and excitations with unprecedented
precision. In particular, resonant inelastic Xray scattering has revealed a lot of hiddenfeatures of high-Tc superconductors and spinorbit insulators. It is of great interest to ﬁnd
novel quasi-particles and understand complex
excitations in correlated materials.
Spin-orbit insulators emerge from 5d electron systems with strong spin-orbit coupling.
Coulomb repulsion between electrons in 5d
atomic orbitals is not strong enough to suppress charge degrees of freedom completely.
As a result, the charge gap can be comparable to the magnon-band width, as experimentally measured in strontium iridates. Suchelectron systems can be described by the Hubbard model with an intermediate Coulomb repulsion, which is parameterized by U .
Strontium iridates, which can be approximated by a half-ﬁlled J = 1/2 singleband Hubbard model, provide an excellent
play ground of charge-spin-orbital entanglement. The dimensionality drastically changes
the electronic property: the single layer material Sr2 IrO4 is a Mott insulator with 400
meV charge gap, the bilayer material Sr3 Ir2 O7
is an insulator in the intermediate coupling
regime with 130 meV charge gap, and thethree-dimensional material SrIrO3 is a correlated metal. These systems provide one of thebest opportunities to study the metal-insulatortransition. Particularly, the bilayer system is
expected to show a characteristic excitation
emerging in the crossover regime. Moreover,in the view of engineering, the magnetic energy
scale enhanced in the crossover regime would
be crucial to realistic devices.
To understand the dynamics of spin-orbit
insulators, we have developed an eﬃcient numerical approach that enables large-scale calculation of dynamical quantities at ﬁnite temperatures in a broad-U region of the Hubbard
model. In our approach, sampling of the auxiliary vector ﬁelds from the Boltzmann distribution is combined with real time evolution of
the self-consistent density matrix.
We used ISSP System B through class
C projects (ID: H31-Ca-0057, H31-Cb-0047).
Our code is implemented with MPI parallelization using GPUs. In particular, GPU programming is crucial to fast matrix-vector multiplication and diagonalization. Applying our approach to the bilayer iridate Sr3 Ir2 O7 , we haveclariﬁed excitonic excitations of the longitudinal or Higgs mode appearing in the dynamical
spin structure factor. Interestingly, excitons
exhibit a signiﬁcant softening as the temperature approaches the magnetic transition temperature, while transverse magnons do not.
Our approach is applicable to a broad rangeof the Hubbard model without restrictions of
the lattice, hopping parameters, and electron
ﬁlling. The developed method is expected topush the boundary of dynamical simulation ofcorrelated electron systems and lead to theidentiﬁcation of hidden excitations emerging
from charge-spin-orbital entanglement.

251

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

Computational design of enzymes and ligand-binding
proteins
Kazuhisa OHARA1, Mio SANO1, Shunji SUETAKA1, Nao SATO1, Yoshiki OKA1,
Koji OOKA2, Yuuki HAYASHI1, and Munehito ARAI1,2
Department of Life Sciences, Graduate School of Arts and Sciences, and1

Department of Physics, Graduate School of Science,

2

The University of Tokyo, Komaba, Meguro, Tokyo 153-8902
Proteins have various functions, including

top mutants, created the mutant proteinsbinding and catalysis. Thus. proteins areexperimentally, and measured their activities.

potentially useful for industry and medicine if

We found that 20 out of 34 mutants had thetheir functions are sufficiently improved or

activity higher than the wild-type protein.

modified. However, because the protein folding

Moreover, for the top 2 mutants, we succeeded

problem remains unsolved, development ofin improving the activity ~4-fold compared

computational protein design strategies haswith the wild type. Thus, our strategies forbeen a great challenge. To solve this issue, here

rationally improving activity of an enzyme will

we attempted to use the Rosetta software suite

be applicable to industrial use.

to computationally design various proteins,(2) Second, to develop a protein that can be

including enzymes and ligand-binding proteins.

used as an anti-allergy drug, we designed two

The static version of the Rosetta 3.8 suite wastypes of proteins that inhibit the interaction ofused in the following calculations [1,2]. In

interleukin-33 (IL-33) with the ST2 receptor of

addition, to solve the protein folding problem,

type 2 innate lymphoid cells [4]. One template

we developed a statistical mechanical model offor the inhibitor design was a fragment of theprotein folding that can explain the folding

ST2 protein. Rosetta calculation was used topathways of multi-domain proteins.

reduce the number of disulfide bonds in the(1) First, we aimed to rationally improve afragment, to make it easier to purify the protein

catalytic activity of an enzyme. By optimizingthat was expressed in Escherichia coli. As

theenzyme,expected, we could obtain the mutant of thedihydrofolate reductase [3], and its ligands, we

ST2 fragment having a reduced number ofcomputationally designed many mutants of thedisulfide bonds as well as the ability to bind IL-

enzyme that were expected to improve its

33. Another strategy to design the IL-33–ST2

catalytic activity. Subsequently, we selected 34

inhibitors was to make mutants of IL-33 thatinteractions

betweenan252

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyobinds tightly to one of the two IL-33 binding(5) Finally, to solve the protein folding

sites on ST2 but does not bind to another site.

problem, we developed a statistical mechanicalThe calculations were already done, and themodel of protein folding. The WSME modelexperimental verification is now under way.

has been known to be able to explain the(3) Third, we rationally designed inhibitors

folding pathways of small proteins and isfor the interactions between the KIX domain of

considered to be a promising model to solve theaandprotein folding problem [6]. However, thetranscriptional activators [5]. We designed andmodel could not be applicable to large, multi-

createddomain proteins, and this problem has beentranscriptionalthecoactivatormutantCBP

fragmentsofthetranscription activator MLL, which bind moreunsolvedtightly than the wild type. We also designed theintroducing a new Hamiltonian in the WSME

KIX mutant that can bind more tightly to themodel and computationally calculating thetranscription activator c-Myb than the wild-

closed form of the partition function, we

type KIX. Since the protein-protein interaction

succeeded in predicting the folding pathways of

between KIX and c-Myb is responsible for themulti-domain proteins. Thus, our model isblood cancer, leukemia, our results may bepromisinguseful for developing anti-cancer drugs.

description of protein folding mechanisms.

for

more

for

a

than

20

unifiedyears.

By

theoretical(4) Fourth, we attempted to modify asubstrate specificity of a ligand-binding proteinusing the Rosetta software. Here, an ATP-

Referencesbinding protein was designed to bind GTP[1] B.J. Bender, et al.: Biochemistry, 55 (2016)

instead of ATP. For experimental verification,

4748-4763.

we selected more than 20 mutants that were

[2] R. F. Alford, et al.: J Chem Theory Comput.

predicted by Rosetta to bind GTP more tightly

13 (2017) 3031-3048.

than ATP. However, the experimental results[3] M. Arai, et al.: J Mol Biol. 410 (2011) 329-

showed that none of the mutants acquired a342.

new function to bind a different ligand. It might[4] X. Liu, et al.: Proc Natl Acad Sci USA. 110

be possible that dynamics of both the protein(2013) 14918-14923.

and ligand should be taken into account more

[5] R. N. De Guzman, et al.: J Mol Biol. 355

extensively to accommodate the difference in(2006) 1005-1013.

the size and shape of ligands, although it[6] M. Sasai, et al.: Biophys Physicobiol. 13

requires much more calculations.

(2016) 281-293.

253

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

Event-chain Monte Carlo method applied to
non-equilibrium relaxation analysis
Y. Ozeki and Y. Osada
Graduate School of Informatics and Engineering, The University of Electro-Communications

The event-chain algorithm is investigated for
its applicability to the nonequilibrium relaxation(NER) method and the estimation of the dynamical exponent z, which indicates an eﬃciency to perform Monte Carlo (MC) simulations. The eventchain method has introduced for the MC simulation with a multi-spin-ﬂip algorithm [1,2]. It shows
an eﬃcient performance for the relaxation even in
frustrated systems and gives estimation of the dynamical critical exponent z ∼ 1 [1], which indicatesa faster dynamics as compared with those with asingle-spin-ﬂip algorithm such as the the Metropolis one.
The NER method is an eﬃcient numerical technique for analyzing equilibrium phase transition [3].
It provides the critical temperature and critical exponents accurately for second-order transition systems, and has been used successfully to study various problems, including frustrated and/or random
systems. In the NER analysis, the equilibration
step is not necessary. Simulation is made only up
to steps when the asymptotic behavior indicatesthe equilibrium state. Thus, one can analyze large
systems as compared with equilibrium simulations.
This advantage becomes more eﬀective for slowrelaxation systems.
We analyze the classical XY model in three dimensions;∑
Si · Sj ,(1)
H = −J

as those in the thermodynamic limit. In otherwords, the ﬁnite size eﬀect must be avoided in
the observed regime of relaxations. Thus, ﬁrst,
we check the size-dependence for the relaxation of|m(t)|. Calculations are carried out at around the
expected transition temperature T / = 2.202 on 43
to 350 − 3 simple cubic lattices up to the observation time of 200 MCSs. Note that the temperatureis measured in the unit of J/kB hereafter. About
840 samples are taken for statistical averaging. The
result is plotted in Fig. 1. It is conﬁrmed that thesize dependence disappears for L ≥ 300 up to 200
MCSs. It is remarked that the relaxation simulated by the event-chain algorithm is much fasterthan that by the Metropolis one, which is an advantage. Furthermore, the appearance of the size
dependence also becomes faster. In the Metropoliscase, it does not appear at least in 104 MCSs for
L = 200 lattice, while it does at 100 MCSs in theevent-chain case. The latter property is an disadvantage, since we must simulate on larger lattices.
Anyway, in the present study, we will use L = 300
lattice for the analysis of phase transition.

⟨ij⟩

where summation for ⟨ij⟩ is taken over all nearestneighboring sites on a simple cubic lattice. The initial state of relaxation is prepared as the all-aligned
state (1, 0). Since the relaxation is much fast in thesimulation with the event chain analysis, we cannot
use the magnetization ⟨Six ⟩ as the dynamical orderparameter, which is used in the NER analysis with
the Metropolis algorithm. Thus, we estimate therelaxation of the absolute value of the magnetization,|m(t)| = (⟨Six ⟩2 + ⟨Siy ⟩2 )1/2 .
(2)In the standard NER analysis, the asymptotic
behaviors of relaxations are used to analyze thephase transition, and are necessary to be regarded

254

Figure 1: Size dependence for the relaxation of|m(t) at T = 2.202 for the 3D XY model.
The relaxation of the order parameter |m(t|) arecalculated for 2.182 ≤ T ≤ 2.222 on the 3003 simple cubic lattice up to an observation time of 200

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of TokyoY = τλm(t,J)

10

1

0.1
0.1

1

X = t/τ

Figure 3: Scaling plot for the data in2.204 ≤ T ≤2.222 in Fig. 2
10

MCSs. About 864 samples are taken for statistical
averaging. The estimations are plotted in Fig. 2.
Using the improved dynamical scaling analysis for
the NER data [4], we analyze the dynamical scalingform,|m(t, T )| = τ −λ Y [t/τ ],(3)

Y = τλm(t,J)

Figure 2: Relaxation of |m(t) for the 3D XY modelin 2.182 ≤ T ≤ 2.222 with the interval ∆T = 0.002.

1

0.1
0.1

1

X = t/τ

where τ (T ) is the relaxation time and λ is a dynamical critical exponent. In Fig. 2, the relaxation
curves show downward trend in 2.204 ≤ T ≤ 2.222,
which indicates the paramagnetic phase. Thus, applying the scaling form (3), the data for these temperatures collapse as in Fig. 3 with Tc = 2.20216
and z = 1.461. Furthermore in Fig. 2, the relaxation curves show upward trend in 2.182 ≤ T ≤2.200, which indicates the ferromagnetic phase.
Thus, applying the scaling form (3), the data forthese temperatures collapse as in Fig. 4 with
Tc = 2.20183 and z = 1.478.
The estimated transition temperature Tc = 2.202
is consistent with Tc = 2016 obtained previously[5]. The estimated dynamical exponent z ∼ 1.47
is greater than z ∼ 1 estimated previously forthe event-chain algorithm [1], while it keeps stilla smaller value than z ∼ 2 which is expected forthe short-range algorithm such as the Metropolis
one. Consequently, the event-chain algorithm canbe used in the NER analysis with an eﬃcient relaxation performance, and would be applicable tovarious slowly relaxing problems such as frustrated
and/or random systems.

255

Figure 4: Scaling plot for the data 2.182 ≤ T ≤2.200 in Fig. 2

References[1] Y. Nishikawa et. al., Phys. Rev. E92 063306(2015)[2] M. Michel et. al., Europhys. Lett. 112
20003 (2015)[3] Y. Ozeki and N. Ito, J. Phys.: Math.
Theor. 40 R149 (2007).
[4] Y. Echinaka and Y. Ozeki, Phys. Rev. E94
043312 (2016).
[5] A. P. Aloysius and M. Harsenbusch,
Physica A201 593 (1993)

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

Entanglement, Quantum Phase Transition, andQuantum Dynamics in Quantum Many-body Systems
Synge Todo
Department of Physics, University of Tokyo, Tokyo 113-0033, Japan
Institute for Solid State Physics, University of Tokyo, Kashiwa, 277-8581, Japan

We have developed various novel numerical methods for quantum/classical many-bodysystems with strong correlations. An eﬃcientquantum Monte Carlo method has been developed to calculate the entanglement entropyof quantum spin systems. The entanglement
entropy is an essential index that characterizes quantum correlations and the topologicalorder in quantum many-body systems. Our
method works directly at absolutely zero temperature, and thus it is entirely free from systematic errors coming from extrapolation to
the zero-temperature limit, unlike the previous
methods.
We also have developed a nonlocal-update
quantum Monte Carlo method for quantum
dimer models on general lattices. By implementing the technique recently suggested forquantum Monte Carlo simulation on quantum
dimer model with the help of priority queue,
which can help ﬁnd the next operator for updating eﬃciently, we can simulate on the ﬁnite temperature region of the quantum dimer
model. By the ﬁnite-size scaling, the phaseboundary between low-temperature columnar
phase with the maximum number of parallel
neighboring dimers, and disordered phase athigh temperatures can be decided in a broad
region of dimer-dimer interaction potential,
and the value of which is both consistent of theclassical limit, where a KT phase transition is
known, and also tends to a constant slope near
the RK-point.

The tensor renormalization group method
can eﬃciently calculate physical quantities oflarge-scale classical/quantum systems. However, the existing methods, such as the tensor
renormalization group (TRG) and the higherorder tensor renormalization group (HOTRG),
suﬀer from the huge computational cost in
higher dimensions. To solve this problem, we
have developed the anisotropic tensor renormalization group (ATRG) that can dramatically reduce the computational complexity
in high-dimensional systems such as threedimensional quantum systems. In addition,
we further improve the accuracy of ATRG by
combining with techniques such as the tensor
renormalization group method that incorporates bond weights, etc.
Spins coupled with photon degrees of freedom in the cavity exhibit a non-equilibriumphase transition with bistability in the number of photons. We have performed simulation
based on the quantum master equations. Fromthe eigenvalues and eigenstates of the time evolution operator, we discussed the correspondence to the ﬁrst-order phase transition in theequilibrium system. We have calculated thehysteresis loop under the periodic modulationof laser intensity, from the viewpoint of the
Floquet operator, and clariﬁed that a dynamical phase transition phenomenon appears withrespect to the period of laser intensity.

In the meantime, the real-space renormalization group methods based on the tensor networks have been widely used in recent years.

256

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

Superconducting dome in electron doped FeSeintercalates
Harald JESCHKE
Research Institute for Interdisciplinary Science, Okayama University
3-1-1 Tsushima-naka, Kita-ku, Okayama 700-8530

The van-der-Waals gap of iron chalcogenide
superconductors can be intercalated with a variety of inorganic and organic compounds thatmodify the electron doping level of the ironlayers. In Lix (C3 N2 H10 )0.37 FeSe, a dome inthe superconducting transition temperature Tc
has been reported to occur in the doping rangeof x = 0.06 to x = 0.68 [1]. We used massive parallel computer simulations to address
magnetic and superconducting properties ofLix (C3 N2 H10 )0.37 FeSe [2]. With a combinationof density functional theory and spin fluctuation theory, we capture the evolution of superconducting transition temperatures theoretically. We clearly demonstrate how the changing electronic structure supports an increasing
superconducting Tc (Figure 1). The suppression of Tc at high doping levels can, however,
only be understood by analyzing the magnetic
tendencies. We determine the Heisenberg exchange couplings between the iron moments;
superconductivity is realized out of these underlying interactions. We can show that withdoping, they evolve from stripe-type at low
doping to bicollinear at high doping. With increasing carrier density, magnetic interactions
mutate from FeSe-like to FeTe-like, explaining
a Tc suppression at high doping.

Figure 1: Structure, band structure andleading eigenvalues of the gap equation forLix (C3 N2 H10 )0.37 FeSe.
[2] M. Shimizu, N. Takemori, D. Guterdingand H. O. Jeschke, Importance of Fermi
surface and magnetic interactions for thesuperconducting dome in electron doped
FeSe intercalates, under review (2020).

References[1] R. J. Sun et al., Phys. Rev. B 98, 214508(2018).

257

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

Properties of electron doped iron based
superconductor LaFe2As2
Harald JESCHKE
Research Institute for Interdisciplinary Science, Okayama University
3-1-1 Tsushima-naka, Kita-ku, Okayama 700-8530

For more than a decade after the discovery of Fe-based superconducting pnictides it
seemed that superconductivity existed in a relatively narrow range of dopings away from thenominal Fe2+ valency. One of the few examples of strongly overdoped (up to Fe2.5+
valency) pnictides was provided by the 1111
material LaFeOAs, with up to 50% of O2−
replaced by H− .
Recently, another compound with formally Fe2.5+ has been synthesized [1], LaFe2 As2 . It was found experimentally that the material can exist in two
distinctly different crystallographic phases,“collapsed-tetragonal” (CT) and “uncollapsedtetragonal” (UT). We used massive parallel computer simulations to address structural, magnetic, and electronic properties of
LaFe2 As2 [2]. We show that the UT phase, as
opposed to the CT one, bears a strong shortrange magnetism of the stripe type driven
by the next-nearest-neighbor exchange. The
structural changes are triggered by magnetic
collapse, as in CaFe2 As2 . Even though the orbitals relevant for the low-energy physics are
not the usual dxz and dyz , the electron pockets do not exclude the usual spin-fluctuation
driven mechanism with an overall s symmetry (Figure 1). Importantly, La in this compound assumes a non-integer valence closer to
2.7+ than to 3+, corresponding to doping of∼ 0.35e, rather than 0.5e. It is thus overdoped,
but not dramatically.

Figure 1: Susceptibility and Fermi surface of
LaFe2 As2 .

References[1] A. Iyo et al., J. Phys. Chem. Lett. 10,1018 (2019).
[2] I. I. Mazin, M. Shimizu, N. Takemori, and
H. O. Jeschke, Novel Fe-Based Superconductor LaFe2 As2 in Comparison with Traditional Pnictides, Phys. Rev. Lett. 123,
267001 (2019).

258

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

Novel phases and dynamics in frustrated spin systems
Tsuyoshi OKUBO
Department of Physics, The University of Tokyo
7-3-1 Hongo, Bunkyo-ku, Tokyo 133-0033

Frustrated interactions in spin systems often
induce novel phenomena, e.g., non-magnetic
ground states including spin liquid and topological orders. Conventionally, geometricallyfrustrated spin systems, such as the kagomélattice and the pyrochlore lattice Heisenberg
antiferromagnets have been investigated actively. Recently, a diﬀerent kind of frustrated
interaction introduced by the strong spin-orbit
coupling has attracted much interest. For example, in Na2 IrO3 and α-RuCl3 , the dominantmagnetic interaction induced by the strong
spin-orbit coupling is considered to be the Kitaev interaction: it is an Ising like anisotropicinteraction with diﬀerent spin component depending on the bonds. When we consider theS = 1/2 quantum spin model on the honeycomb lattice with only the Kitaev interaction (Kitaev model), its ground state is a spin
liquid state without any spontaneous symmetry breakings [1]. This fact indicate real compounds might have a spin liquid ground states.
However, real compounds contain other interactions such as the Heisenberg interaction, further neighbor interactions, or the oﬀ-diagonalinteractions, and they can induce magnetic
long range orders at low temperature.
In this year, in order to clarify eﬀects of additional interactions in the Kitaev materials,
we numerically investigated the ground stateand ﬁnite temperature properties of the models with the Kitaev interaction and the Γ interactions (the Kitaev-Γ model) by means of
tensor network methods. The Γ term is anoﬀ-diagonal interaction and the Hamiltonianof the model is given by
∑
H=
Hγ ,

(1)

γ∈x,y,zwhere for γ = z,∑
Hz =[KSiz Sjz + Γ(Six Sjy + Siy Sjx )].

(2)⟨i,j⟩zFor γ = x, y, we consider similar interaction
with cyclic rotation of x, y, z components. We
consider the honeycomb lattice and γ = x, y, z
corresponds to three directions of the bonds.
This model is considered to be an eﬀectivemodel for α-RuCl3 [2, 3].
Firstly, we investigated the phase diagram of
the Kitaev-Γ model by using a two-dimensional
inﬁnite tensor network approach. We represented the ground state wave function as products of small tensors, and optimize each tensor element so that the energy expectation
value becomes smaller. By assuming severaltransnational symmetries of the ground state,
we can represent inﬁnitely large system by using ﬁnite numbers of tensors. As a networkof tensors, we used the inﬁnite tensor productstates (iTPS), which is also called as the inﬁnite projected entangled-pair states (iPEPS).
As the optimization algorithm, previously, we
have used the imaginary time evolution. However, it have become clear that the imaginary
time evolution often trapped at local minima
and it was very hard to obtain novel spin liquid state from random initial states. Thus, thisyear, we used a diﬀerent type of optimization,so called variational optimization [4]. We conﬁrmed that by using the variational optimiza-

259

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyotion, we correctly obtained the Kitaev spin liquid state at the pure Kitaev model (Γ = 0)even from the random initial states. Then weapplied the variational optimization to ﬁnite
Γ and draw the phase diagram. Although theobtained variational energies were largely improved from those obtained by the imaginary
time evolution, the phase diagram was almost
unchanged; except for the vicinity of the pure
Kitaev model, the ground states were magnetically ordered states. This observation was contrasted with the phase diagrams obtained by
the exact diagonalization (ED) [2] or by theinﬁnite density matrix renormalization group(iDMRG) [3]; they concluded much wider spin
liquid phases. In order to clarify nature of thephase diagram further, we might need larger
scale calculation with iTPS.

[3] M. Gohlke, G. Wachtel, Y. Yamaji,
F. Pollmann, and Y. B. Kim: Phys. Rev.
B 97 (2018) 075126.
[4] P. Corboz:
035133.

Phys. Rev. B 94 (2016)[5] P. Czarnik, J. Dziamaga, and P. Corboz:
Phys. Rev. B 99 (2019) 035115.
[6] A. Kshetrimayum, M. Rizzi, J. Eisert, and
R. Orús: Phys. Rev. Lett. 122 (2019)070502.
[7] J. Nasu, M. Udagawa, and Y. Motome:
Phys. Rev. B 92 (2015) 115122.

Secondly, we investigated ﬁnite temperature
properties of the model. In order to calculate
physical quantities at a ﬁnite temperature, we
represented the density matrix of the system intensor network representations. We employed
two types of the tensor network representations: the local puriﬁcation, e.g., Ref. [5], andthe direct tensor product operator (TPO) [6].
By applying them to the pure Kitaev model,
we found that the former representation failed
to reproduce the expected lower temperature
speciﬁc heat peak [7], while the latter TPO
representation reproduced it, at least qualitatively. This observation indicated that thedirect TPO representation was more eﬃcientthan the local puriﬁcation. The applicationsof the developped technique to the model withﬁnite Γ is on going.

References[1] A. Kitaev: Ann. Phys. 321 (2006) 2.
[2] A. Catuneanu, Y. Yamaji, G. Wachtel,
H.-Y. Kee, and Y. B. Kim: npj Quantum
Materials 3 (2018) 23.

260

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

Thermal effects on quantum frustrated magnetisms
Tokuro Shimokawa
Okinawa Institute of Science and Technology Graduate University,
Tancha, Onna-son, Kunigami-gun, Okinawa 904-0495
moderate

We know well that quantum frustrated

temperatures

exhibitingamagnetisms often exhibit rich many-bodycharacteristic ring-like structure in the equal-

physics, but we also know difficulties totime spin structure factor. These obtained

investigatetheresults

fluctuations

numerically,effects

of

the

hence

thermalfreedom in the intermediate plaquette groundstate, we found the signature of the finite-

investigating the finite-temperature properties

temperature phase transition associated with theof S=1/2 bilayer-breathing-kagome (BBK)

breaking of the mirror symmetry at moderate

Heisenberg model for Ca10Cr7O28 and S=1/2

temperatures via specific heat and correlation

Shastry-Surtherland (SS) Heisenberg model for

functions calculations. We also investigate the

SrCu2(BO3)2.

joinedrelationship the transition and the existence ofrecently the family of quantum spin liquid

the degenerated excited states by the exact

candidates [3] in spite of the complexity of its

diagonalization method. [6]former

[1-2]

explaining

for

The

methodfor

The SS model: Owing to the Ising degrees of

We use a state of art numerical method,
typicality

reasonable

qualitatively the nature of the Ca10Cr7O28. [5]thermally

induced novel phases still remain unexplored.
quantum

aremagnetHamiltonian. The latter magnet has been known
well, but recent experiments under highReferences

pressures [4] could reveal successfully the[1] A. Hams and H. De Raedt, Phys. Rev. E

existence of the intermediate plaquette state, so

62, 4365 (2000).

the theoretical studies at finite temperatures are[2] S. Sugiura and A. Shimizu, Phys. Rev.

required.

Lett. 111, 010401 (2013).

The BBK model: We computed the several[3] C. Balz, et al, Nat. Phys. 12, 942 (2016).

physical quantities such as specific heat,[4] M. E. Zayed et al, Nat. Phys. 13, 962 (2017).

magnetization curve, and equal-time spin

[5] T. Shimokawa, R. Pohle, N. Shannon, in

structure factor at finite-temperature with the

preparation.

quantum typicality method and found the[6] T. Shimokawa, in preparation.

possible existence of a spiral spin liquid state at261

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

Numerical studies on magnetization process of theKitaev spin liquid
Kota IDO
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581

The Kitaev model is one of the prominent examples of quantum many-body systems
where quantum spin liquids emerge[1]. In the
Kitaev model, the nearest-neighbor spins interact with each other via bond-dependent Isingtype interactions. Although the Kitaev modelseems to be an artiﬁcial model, it has been proposed that the Kitaev model can be realized inreal materials[2]. Inspired by the proposal, it
is now one of the hot topics in the condensed
matter physics how to synthesize and characterize the Kitaev-like materials.
The Kitaev model under a magnetic ﬁeld has
been intensively studied, because the magnetization processes oﬀer useful information forcharacterizing the nature of the quantum spin
liquid in the Kitaev model. Recently, it wasreported that an intermediate gapless state
appears in the antiferromagnetic coupling Kitaev honeycomb model under a magnetic ﬁeld.
For example, exact diagonalization (ED) calculations showed that the magnetization curve
has two abruptly changes by increasing theﬁeld[3, 4]. This feature was observed by othernumerical methods such as the density matrix
renormalization group (DMRG)[5] and Majorana mean-ﬁeld (MF) theory[6]. However, ED
and DMRG studies were done in small system size and quasi-one-dimensional system on
cylindrical boundary condition, respectively.
It is desired to perform accurate analysis basedon a numerical calculation method applicableto the large two-dimensional system.
In this study[7], by using the variational

Monte Carlo (VMC) method, we investigate
whether the intermediate state exists as theground state of the Kitaev model in two dimensions in a [001] magnetic ﬁeld. We used a generalized Bardeen-Cooper-Schrieﬀer trial wavefunction with Jastrow factor, which can exactly represent the ground state of the pure Kitaev model fermionized by the Jordan-Wignertransformation[8, 9]. Our benchmark for asmall system show that this trial wavefunction well reproduces the exact magnetization
process. In applications to the antiferromagnetic Kitaev model with large system sizes, adouble-peaked structure is found in the magnetic susceptibility, which signals the existenceof the intermediate phase. In the intermediate strength of the magnetic ﬁeld, we ﬁnd thatthe local Z2 gauge ﬁeld takes non-integer values, which is consistent with the spin liquid
phase with ﬂuctuating ﬂux for [111] case. Our
results also show that the momentum distribution of the complex fermions changes topologically at the ﬁrst transition point, suggesting that this is a topological phase transition
resembling the Lifshitz transition. We cannot ﬁnd these properties in the ferromagnetic
Kiteav model. Our results show that in theantiferromagnetic Kitaev model, the intermediate state is stable against many-body correlations beyond the MF approximation.
Our VMC results also suggest that it is a
useful and accurate way to analyze the Kitaev
model in a magnetic ﬁeld by using fermionization of the localized spin diﬀerent from the262

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

spinon representation. Using other types of thefermionizations, it would be possible to accurately analyze the magnetization process of theKitaev model with additional interactions suchas the Heisenberg term and the Γ term within
the framework of VMC. Our study oﬀers a ﬁrm
basis for such extended treatments.
In addition, we have proposed how to calculate the charge (density) dynamical structure
factors for the ground states of correlated electron systems based on the VMC method[10].
It is a future issue to study complex fermion
dynamics of quantum spin liquids by using our
proposed method.

References[1] A. Kitaev: Ann. Phys. 321 (2006) 2.
[2] G. Jackeli and G. Khaliullin: Phys. Rev.
Lett. 102 (2009) 017205.
[3] Z. Zhu, I. Kimchi, D. N. Sheng, and L.
Fu: Phys. Rev. B 97 (2018) 241110(R).
[4] C. Hickey and S. Trebst: Nat. Commun.
10 (2019) 530.
[5] M. Gohlke, R. Moessner, and F. Pollmann: Phys. Rev. B 98 (2018) 014418.
[6] J. Nasu, Y. Kato, Y. Kamiya, and Y. Motome: Phys. Rev. B 98 (2018) 060416(R).
[7] K. Ido, and T. Misawa: Phys. Rev. B 101(2020) 045121.
[8] H.-D. Chen and Z. Nussinov: J. Phys. A
41 (2008) 075001.
[9] X.-Y. Feng, G.-M. Zhang, and T. Xiang:
Phys. Rev. Lett. 98 (2007) 087204.
[10] K. Ido, M. Imada, and T. Misawa: Phys.
Rev. B 101 (2020) 075124.

263

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

A new block-spin transformation in Potts models
Yusuke TOMITA
College of Engineering, Shibaura Institute of Technology
Saitama, Saitama 337-8570

Magnetic susceptibilities for 2D Potts models through the new block-spin transformation are shown in Fig. 1. The original system size is 8192. Slopes of lines in Fig. 1 are
drawn by using exact critical exponents, γ/ν’s.
The result indicates the new transformation
rule works well for each q-state Potts model.
While the newly proposed transformation utilizes the fractal property of Kasteleyn-Fortuin
cluster on critical points, I also examined
the renormalization properties at oﬀ-criticalpoints. Figure 2 shows numbers of clusters forthe 2D Ising model for T = 0.95 × Tc , Tc , and1.05 × Tc . At the critical point, the number ofclusters is well ﬁtted by a power function withexponent 4/(2 − β/ν). On the other hand, thenumbers of clusters at oﬀ-critical points suddenly drop at some block-spin transformation
steps. A renormalized cluster expands almost
whole the system at a some ﬁnite renormalized
step when the temperature is lower than thecritical temperature, while all of the clusters
disappear at a some ﬁnite renormalized stepwhen the temperature is higher than the critical temperature. Therefore, we can estimate
critical points by the linearity of the numberof clusters.
107

q⟨1

106

q⟨2

q⟨4

q⟨3

⟨χ⟩

105
104
103
102
101
100

213211 29 27 25 23 21213211 29 27 25 23 21213211 29 27 25 23 21213211 29 27 25 23 21

L

Figure 1: Plot of the susceptibilities for the
2D Potts models. Slopes of lines are drawn by
using exact critical exponents, γ/ν

107

T ⟩ TcT ⟩ 0⟨95 × Tc
T ⟩ 1⟨05 × Tc

106
105
104

⟨nc ⟩

In usual block-spin transformation, the majority rule is mainly used as a renormalization
scheme. Many studies using the block-spintransformation have shown that spin conﬁgurations are well renormalized well by the majority rule. But the legitimacy of the majority
rule has not been understood. In this study,
I propose a block-spin transformation rule,
which utilizes the fractality of the KasteleynFortuin cluster on the critical point. The
new block-spin transformation is constructedto conform Kasteleyn-Fortuin clusters to theusual ﬁnite-size scaling form, L−β/ν .

103
102
101
100
10−1

213

211

29

27

L

25

23

21

Figure 2: Plot of the numbers of clusters forthe 2D Ising (q = 2) models. A slope of line isdrawn by using exact critical exponents, 4/(2−
β/ν)

264

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

Screening for Thermal Functional Materials using Materials
Informatics
Junichiro SHIOMI
7-3-1, Hongo, Bunkyo-kun, The University of Tokyo, Tokyo 113-8656
Controllingthethermalconduction, the overlap results in the higher

boundary

conductance (TBC) between copper and carbon

value

in

TBC

compared

crystals is important since it can bottleneck thecopper/diamond interface.

with

the

thermal conductivity when reinforcing copper

To improve the TBC at the copper/diamondwith carbon-crystals fillers, namely diamond orinterface, the covalently bonded SAMs withgraphite, to develop heat sinks and spreadersdifferent chain lengths are installed at theneeded for the thermal management. In the work

interface. The resultant values of TBC withof this year [1], we demonstrate the spectral

SAMs of various chain lengths are plotted in Fig.

control of TBC between cooper and carbon

2. The data at n=0 refers to the TBC at

crystal by self-assembled monolayers (SAMs)

copper/diamond and copper/graphite interfaces

using the non-equilibrium molecular dynamics

without SAMs. For the chain length studied insimulation.

this work, the value of TBC increases
formonotonically with the chain length for chaincopper/diamond and copper/graphite interfaces

length shorter than C15 and reaches a plateau

are 9.06 ± 0.16 MW/m2K and 81.76 ± 14.85

afterward. The highest value of TBC approaches

MW/m2K, respectively. The result indicates theto that of bare copper/graphite interface.

The

resultantvaluesofTBC

TBC at copper/graphite interfaces is almost anorder of magnitude greater than that at
copper/diamond interface. In order to interpretthe behavior, the vibrational density of states(VDOS) of copper/diamond and copper/graphite
interfaces are calculated, as shown in Fig. 1. The
result indicates the out-of-plane VDOS ofgraphite surface provides the overlap in phonon

Fig. 1: Vibrational density of states of various

vibrational modes with that of copper surface. In

components of (a) copper/diamond and (b)view of elastic energy transport, which is usuallycopper/graphite interfaces.

the main contribution to interfacial heat265

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

useful to improve the thermal conductivity ofmetal/carbon-crystal composite materials in theindustry applications.
In summary, we have demonstrated (1) thepoor TBC at cooper/diamond interface is mainly
due to the mismatch of VDOS between cooper
Fig. 2: Thermal boundary conductance as aand diamond (2) the installation of SAMs atfunction of chain length of SAMs. n=0 refer tocooper/diamond interface can improve the TBC

the interfaces without SAMs.

significantly.

The TBCshowsastrongdependence on chain length that the TBC
increases with the increasing chain length for

Both the spectral and structural analysis arethechain length shorter than C15 and reaches aunderlying mechanism. The analysis identifies

plateau afterward. The materials informatics

that this is because the SAMs becomes softer

technique will be helpful for the design of SAMs

with increasing chain length due to disordering

toof the collective SAMs structure, enhancing themetal/carbon-crystal composite materials.

implemented

tofurtherunderstandenhancethethermalconductivityofspectral overlap with copper at low frequency
vibrational modes. In the future, the materials

References

informatics techniques can be conducted to find[1] S.-W. Hung, S. Hu, and J. Shiomi, ACS Appl.

out the optimum design of SAMs to enhance theElectron. Mater. 1 (2019) 2594.

interfacial heat transfer. The obtained results are

266

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

Wave function analyses in disordered electron
systems
TOMI OHTSUKI1
TOHRU KAWARABAYASHI2
KEITH SLEVIN3
KOJI KOBAYASHI4
TOMOHIRO MANO1
1) Dept. Phys., Sophia University, Chiyoda-ku, Tokyo 102-8554, Japan
2) Dept. Phys., Toho University, Miyama 2-2-1, Funabashi 274-8510, Japan
3) Dept. Phys., Osaka University, Toyonaka, Osaka 560-0043, Japan
4) IMR, Tohoku University, Sendai 980-8577, Japan

Recent discoveries of Weyl semimetal(WSM) have inspired extensive research ofthese novel topological materials, especiallyfrom the point of the new universality classes,
which has discrete symmetries such as particlehole symmetry. Here we have studied the Anderson transition in chiral orthogonal class aswell as class BDI, the latter showing the nodalline Weyl semimetal phase[1]. We have alsostudied the various quantum phase transitions
by analyzing the wave functions via convolutional neural network (CNN) [2]. The typicalwave functions in disordered quantum systems
are shown in Fig. 1. The method is especially powerful when analyzing the quantum
percolation, where the transfer matrix method
breaks down due to the irregular lattice structure. We have trained CNN to capture the features of wave functions in topological as well
as nontopological phases, and let CNN draw
the phase diagram (Fig.2) [3]. In addition, we
have calculated the distribution of Kondo temperature at the Anderson transition from thelocal density of states obtained by kernel polynomial method[4].
References
1. X. Luo, B. Xu, T. Ohtsuki, R. Shindou: Physical Review B 101, 020202(R), (2020).
2. T. Ohtsuki, T. Mano: Journal of the PhysicalSociety of Japan 89, 022001 (2020).
3. T. Mano, T. Ohtsuki: Journal of the PhysicalSociety of Japan 88, 123704 (2019).

4. K. Slevin, S. Kettemann, T. Ohtsuki: The European Physical Journal B 92, 281 (2019).

Figure 1: Typical 3D wave functions in real
space in various phases. Taken from [2].

Figure 2: Phase diagram of quantum percolation in 3D topological insulator. OI (ordinary insulator), WTI (weak topological insulator), STI (strong topological insulator), andDM (diﬀusive metal).(ijk) indicates the weakindices.

267

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

Numerical study on low-energy statesof quantum spin systems
Hiroki NAKANO
Graduate School of Material Science, University of Hyogo
3-2-1 Kouto, Kamigori-cho, Ako-gun, Hyogo 678-1297, Japan

It is diﬃcult to precisely estimate physicalquantities for many-body problems. Particularly, we often have to tackle such manybody problems in condensed matter physics.
Quantum spin systems are typical many-bodyproblems. Under circumstances, numerical approaches have widely and eﬀectively been employed to examine the systems. Many computational investigations have been carried out
and provided us with useful information anddeeper understanding of the target systems.
It is well known that within the field of
quantum spin systems, three methods are effectively used; the numerical diagonalizations,the quantum Monte Carlo (QMC) simulations,and the density matrix renormalization group(DMRG) calculations. Each of them has advantages and disadvantages at the same time.
In the QMC simulations, large systems canbe treated irrespective of their spatial dimensions although it is diﬃcult to obtain precise
evaluations due to the negative sign problem
when the systems include frustrations. On theother hand, the DMRG method is very useful when spatial dimension of a target system
is one irrespective of whether frustrations ispresent or absent in the system. However,this method is still under development for thecases when the spatial dimension is larger thanone. The numerical diagonalization method
can be applied irrespective of the spatial dimension and the presence of frustrations. Unfortunately, this method has a serious weak
point that it can treat only very small systemsizes. As one way to overcome this disadvantage, we successfully developed a hybrid-typeparallelized code of Lanczos diagonalization[1].
If we use this Lanczos-diagonalization code, it268

is possible to treat various large systems thathave not been previously treated yet within
this method. Therefore, we investigate various quantum spin systems by this method asa primary approach.
In the project in 2019, we tackled two systems. The first one is the S = 1/2 Heisenberg antiferromagnet on the kagome latticewith a specific distortion[2, 3]. We investigated a family of materials with Ti3+ S = 1/2
ions forming distorted kagome lattices. Employing all electron density functional theory
technique, the Heisenberg Hamiltonian was established for each material; finite-size clusters
for the Hamiltonians were numerically diagonalized. Our calculations successfully provided
us with the magnetization processes, which enable us to analyze experimental measurements.
The comparisons were quite good; our theoretical calculations succeeded in predicting the
behavior in fields higher than the maximum
field which the present experiments can reachpresently.
The second one is the large-S Heisenberg
antiferromagnet in one dimension[4]. Particularly, we focused our attention on the cases of
S = 5 and S = 6. The case of S = 5 had beenstudied before in Ref. [1] where systems up to
10 sites were treated. We successfully carriedout our huge-scale parallel calculation in the K
computer using approximately 97% of the total computational nodes without any troubles
and obtained novel results up to 14 sites of
the target system. Note also that the Hilbert
space dimension of 14-site cluster for the S = 5
chain is larger than the case of 20-site cluster
for the S = 2 chain[5]. From the new results,
our analysis gave a more precise estimation of

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

the S = 5 Haldane gap. In this study, we also
attempted to estimate the Haldane gap for thecase of S = 6, which had not been examined
before. From our new estimates of the Haldane gaps for large S, the asymptotic formulaof the Haldane gap was examined; we successfully determined the coeﬃcient in the formulamore precisely.
Our studies contribute to our deeper understandings of the various antiferromagnets andthe nontrivial eﬀect of quantum nature andfrustration in magnetic materials.

References[1] H. Nakano and A. Terai: J. Phys. Soc.
Jpn. 78 (2009) 014003.
[2] H. O. Jeschke, H. Nakano, and T. Sakai:Phys. Rev. B 99 (2019) 140410(R)(1-6)[3] R. Shirakami, H. Ueda, H. O. Jeschke,
H. Nakano, S. Kobayashi, A. Matsuo, T.
Sakai, N. Katayama, H. Sawa, K. Kindo,
C. Michioka, and K. Yoshimura: Phys.
Rev. B 100 (2019) 174401(1-7)[4] H. Nakano, N. Todoroki, and T. Sakai: J.
Phys. Soc. Jpn. 88 (2019) 114702(1-4)[5] H. Nakano and T. Sakai: J. Phys. Soc.
Jpn. 87 (2018) 105002(1-2)

269

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

Efficient Sampling Simulation of
the Soft Modes Significantly Contribute to Protein Properties
Kazuhiro TAKEMURA, Duy TRAN, Hiroaki HATA, Kenichiro TAKABA,
Mohamed Marzouk SOBEH, and Akio KITAOSchool of Life Science and Technology, Tokyo Institute of Technology
Ookayama, Meguro, Tokyo 152-8550
We have developed a method to evaluate

for efficiently selecting near-native decoys

binding free energy differences of complex

using Best Arm Identification (BAI), called

models generated by docking prediction

evERDock BAI[3].

using the all-atom molecular dynamics
simulation and solution theory in the energy
representation (ER), termed evERdock(evaluation with ER method of docking
generateddecoys)[1].

Byevaluating

binding free energy difference with themethod, we previously showed that “nearnative” models similar to crystal structure
are successfully selected as the lowest
energystructuresinseveralcomplexsystems. The method requires relatively
short MD simulations (2 ns) and cancalculate binding free energy difference
among several hundreds of complex models.
We also reported that refinements ofdocking generated decoys by relatively long
MD simulation and averaging the binding
free

energies

obtainedfrommultiplesnapshots further improved the evaluationof evERdock[2]. However, such treatments
required additional computational costs. To

Fig. 1 Flows of (A) uniform sampling and (B)address this problem, we propose a method

BAI (UCBE) based evaluation.

270

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

Let the number of evERdock runs decoys

successfully reduced computational cost by a

are N be M, respectively. If we conduct

factor of 2.20, 1.22, and 4.05 (values in

evERdock runs for all decoys equally

parentheses in Table 1) for 1AY7, 1PPE, and(uniform sampling, US), all decoys are

3SIC, respectively.

evaluated M/N times (Fig.1A). Fig. 1B shows

As shown above, we successfully improved

the flow chart of evERdock-BAI using UCBE

the computational efficiency of the evERdock

as a BAI algorithm. After all decoys areevaluation of protein-protein complex models

evaluated once, the algorithm repeatedly

by combining with BAI algorithm.

calculates UCBE scores, selects the decoy
with

the

highestscore,

andconducts

evERdock run for the selected decoy.
Promisingdecoysfor

exploitationareevaluated by the first term of the score(binding free energy, −∆𝐺 ). The second termof the score contains the number of evERdock
runs for the decoy i (ni) in its denominator.
Fig. 2 Structures of target complexesThus, the score decreases as the number of
evERdock runs for a decoy i increases. As a

Table 1 Summary of results

result, unpromising decoys are sometimes

Complex
NdecoysNnear-native
N95% (US)
N95% (UCBE)

selected, and exploration is realized. Such abalance between exploitation and exploration
is essential to achieve an efficient and

1AY7
11
2
33
15 (2.20)

1PPE
11
4
33
27 (1.22)

3SIC
21
1
336
83 (4.05)

accurate prediction of the correct decoy.
WeappliedevERdock-BAI

to

three

References

protein-protein complex systems (Fig. 2).

[1] K. Takemura, N. Matubayasi, and A.

Decoys are generated by rigid-body dockings.

Kitao, J. Chem. Phys. 148, 105101 (2018).

The number of decoys (Ndecoys), the number of[2] A. Shinobu, K. Takemura, N. Matubayasi,

near-nativeand A. Kitao, J. Chem. Phys. 149,

decoys(Nnear-native),

and

the

number of total evERdock runs to achieve a

195101 (2018).

prediction accuracy 95 % (N95%) using US[3] K. Terayama, A. Shinobu, K. Tsuda, K.

and UCBE are shown in Table 1. For all three

Takemura, and A. Kitao. J. Chem. Phys. 151,

cases, N95% using UCBE were smaller than

215104 (2019).

those using US. The UCBE algorithm271

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

Pseudo-fermion functional renormalization group
approach for frustrated qunatum spin systems
Kiyu FUKUI, Yukitoshi MOTOME, and Yusuke KATO
Department of Physics, The University of Tokyo, Bunkyo, Tokyo 113-0033
Department of Applied Physics, The University of Tokyo, Bunkyo, Tokyo 113-8656
Department of Basic Science, The University of Tokyo, Meguro, Tokyo 153-8902

Frustrated quantum spin systems have beenstudied intensively for very long time [1, 2].
It is impossible to solve them exactly mathematically in general. Therefore, numericalapproaches to investigate frustrated quantum
spin systems are indispensable. Numerical approaches used widely to study them are, forexample, exact diagonalization (ED), quantum
Monte Carlo (QMC), density matrix renormalization group (DMRG). Each method has pros
and cons, of course. Since 2010s, a new method
called pseudo-fermion functional renormalization group (PFFRG or pf-FRG) has been usedgradually [3]. Although PFFRG overestimatesordering tendencies to magnetic orders [4],it can treat large-size systems and stronglyfrustrated systems and detect quantum paramagnetism due to strong frustration.

nian (1 = Heisenberg Hamiltonian with timereversal symmetry), respectively. Typically,
NL = 10-20, Nω = 50-100, and NΛ = 200500. Therefore, it is necessary to exploit parallelization. We submitted jobs parallelized by
openMP+MPI to F144cpu queue using 64-144
nodes and O(1000) MPI processes. One PFFRG calculation with spin susceptibility calculation with vertices obtained by PFFRG takes
one to several days. In addition, we tried to use
GPU nodes by openACC directives in order toaccelerate our PFFRG calculation.

In order to study quantum paramagnetism,
particularly quantum spin liquid, we made up
a PFFRG numerical code and tested it by reproducing previous studies. At first, we madePFFRG codes for SU(2) symmetrical systems(i.e. systems with Heisenberg-type interactions) [3] and benchmarked it. Then we extended our codes to systems with XXZ-type [4]and Kitaev-type [5] interactions.

[2] H. T. Diep, Frustrated Spin Systems(World Scientific, Singapore, 2013) 2nded.

Computational complexity of PFFRG calculation is scaled by O(NL2 · Nω4 · NΛ · Nσ ). NL ,Nω , NΛ , and Nσ are symmetry-reduced systemsize, the number of frequency mesh points, thenumber of energy cut-oﬀ grid points, and relative complexity due to symmetry of Hamilto-

272

References[1] C. Lacroix, P. Mendels, and F. Mila,
Introduction to Frustrated Magnetism(Springer, New York, 2011).

[3] J. Reuther and P. Wölfle, Phys. Rev. B
81, 144410 (2010).
[4] S. Göttel, S. Andergassen, C. Honerkamp,
D. Schuricht, and S. Wessel, Phys. Rev. B
85, 214406 (2012).
[5] J. Reuther, R. Thomale, and S. Trebst,
Phys. Rev. B 84, 100406(R) (2011).

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

Precise measurement of crystal nucleation rate in a
hard-sphere colloidal system
Michio Tateno and Hajime Tanaka
Department of Fundamental Engineering, Institute of Industrial Science, University of Tokyo
4-6-1 Komaba, Meguro-ku, Tokyo 153-8505

Hard-sphere colloids have been intensively
studied as one of the most fundamental systems that show a liquid-crystal transition [1].
In 2001, Auer and Frenkel pointed out [2]a huge mismatch in the crystal nucleation
rate between results of light scattering experiments and numerical prediction based onMonte Calro simulations using the Umbrella
Sampling method. Revealing the physical origin behind this huge discrepancy is one of
the central issues in this field [1]. Hydrodynamic interactions (HI) among colloids are
viewed as one of the possible factors to cause
the discrepancy since simulation results used
for the comparison completely neglect HI. So
far, two simulation studies based on different
models incorporating HI have been performed
to address this problem. However, these studies reported conflicting results: one concluded
that HI might speed up crystal nucleation [3],whereas the other claimed that HI slows down
it [4].
In this project, we study the crystallization
kinetics of hard-sphere colloids by two simulation methods with HI (Fluid Particle Dynamics method: FPD [5]) and without (Brownian
Dynamics method: BD). By comparing two
results, we can investigate the impact of HI
on the crystal nucleation process. By combining the FPD method with fluctuating hydrodynamics, we have shown that FPD simulations can almost perfectly reproduce the motion of colloids under thermal noise at low volume fractions while retaining statistical mechanical consistency between colloids and asolvent. Furthermore, we have confirmed that
this simulation approach can precisely reproduce the phase separation process of colloidal
suspensions observed by confocal microscopy
measurements without any adjustable parameters [6]. However, crystal nucleation is a rare
event, and thus, the resulting computational
cost is quite high. We have overcome this problem by speeding up our code, utilizing ISSP’s
GPGPU implementation services [7].
According to Ref. [2], the theoretical valuesof I change mainly by a small change in the volume fraction φ: for example, I at φ = 52.8 %differs from that at 53.4% by at least 4 ordersof magnitude. Thus, in order to study the nucleation rate accurately, we need to estimate
the volume fraction φ precisely. To realize this,
we first check the volume fraction dependenceof the virial pressure. In Fig. 1a, we show theresults obtained by BD and FPD together withtheoretical predictions (the Carnahan-Starling
equation of a state in the liquid phase and the
Speedy equation of a state in the solid phase).
Here we can see almost perfect agreement between them, supporting the validity of our numerical simulation methods from the thermodynamic point of view.
Now we turn our attention to a dynamic aspect. Here we examine the volume-fraction dependence of the long-time self-diffusion constant, DL as one of the most fundamental
quantities to characterize colloidal dynamics.
Figure 1b shows the results for FPD and BD,
where we can see that DL deviates from thediffusion coefficient for a free isolated particle

273

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

Figure 1: (a) Virial pressure in equilibrium
states. Cross and square symbols represent theresults obtained by BD and FPD simulations,
respectively. Gary curve is the theoretical prediction for a hard-sphere system. (b) Volume
fraction dependence of the long-time diffusion
constant DL scaled by that in the dilute limit
D0 . Cross and circle symbols represent the results obtained by BD and FPD simulation, respectively. The light blue square symbols areexperimental data obtained by dynamic lightscattering measurements for hard-sphere-likecolloids with 5 % polydispersity [8].
D0 for large φ, but the deviation is more significant for FPD than for BD, which depends
on the presence or absence of hydrodynamic
interactions. For example, we can see that
DL in BD is approximately four times larger
than that in FPD at φ = 0.54. In the figure, we also compare the simulation resultswith experimental data obtained by dynamic
light scattering measurements for hard-spherelike colloids with 5 % polydispersity [8]. Here
we can see that there is a mismatch between
FPD and BD, but there is a good agreement
between FPD and the experimental results atleast at intermediate volume fractions. To determine the nucleation rate I, we used an approach to monitor the temporal change of acluster with the largest size. We compare I at
0.54 ≤ φ ≤ 0.545 between BD and FPD, by
using DL and the diameter of colloids as spatiotemporal units, and find that all nucleation
rates collapse onto a single curve.

time diffusive dynamics of the colloids is more
rapidly slowed down by HI with increasing volume fraction. We find that this long-time diffusion controls the nucleation rate and that it
is possible to account for most of the effectsof HI by rescaling with this time scale. The
research results are summarized in Ref. [9].

References[1] T. Palberg: J. Phys. Cond. Matt. 33(2014) 333101[2] S. Auer and D. Frenkel:(2001) 1020

Nature 409[3] M. Radu and T. Schilling:(2014) 26001

EPL 105[4] D. Roehm, S. Kesselheim and A. Arnold:
Soft matter 10 (2014) 30[5] H. Tanaka and T. Araki: Phys. Rev. Lett.
85, (1999) 1338; A. Furukawa and H.
Tanaka: Phys. Rev. Lett. 104, (2010)
245702; A. Furukawa, M. Tateno and H.
Tanaka: Soft Matter 14 (2018) 3738
[6] M. Tateno and H. Tanaka: npj Comput.
Mater. 5 (2019) 40[7] M. Tateno, K. Takae and H. Tanaka,“GPGPU implementation of Fluid Particle Dynamics method”, Activity Reportin 2017[8] P. N. Segre and S. P. Meeker, P.N. Pusey
and W. C. K. Poon: Phys. Rev. Lett. 75(1995) 958[9] M. Tateno, T. Yanagishima, J. Russo, and
H. Tanaka: Phys. Rev. Lett. 123, (2019)258002In summary, to examine the influence of HI
on crystallization kinetics of hard-sphere colloids, we precisely compare simulation resultswith and without HI. We show that the long-274

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

Impact of hydrodynamic interactions on colloidal
dynamics near a critical point
Michio Tateno and Hajime Tanaka
Department of Fundamental Engineering, Institute of Industrial Science, University of Tokyo
4-6-1 Komaba, Meguro-ku, Tokyo 153-8505

The behaviors of simple liquids near a critical point are well-understood both from a theoretical and experimental point of view. The
density fluctuation decays exponentially withsingle relaxation time. Contrary to this knowledge, several studies on colloidal suspensions
using light scattering experiments report thattime correlation function of concentration fluctuations in the vicinity of the critical pointshow nonexponential decay [1].
The central aim of this project is to reveal
the physical mechanism underlying the nonexponential decay. As a possible candidate tobring about this unconventional dynamical behavior in colloidal suspensions, we focus on hydrodynamic interactions (HI), which is an essential difference between simple liquids andcolloidal suspensions. We use two different
simulation methods; the Brownian dynamics(BD) method neglecting HI and Fluid particle
dynamics (FPD) method incorporating HI [2].
We compare the results from two simulation
methods and aim to examine the effect of HI
on dynamical critical phenomena.
In order to capture dynamical critical phenomena numerically, a long time simulation isrequired, since the relaxation time for density
fluctuation grows dramatically as approaching
to a critical point. However, FPD method isnumerically costly because it is based on the direct computation of Navier-Stokes equations.
To overcome this, we implemented GPU on
our FPD codes utilizing a service provided by
ISSP. To further suppress numerical cost, weFigure 1: Temperature-volume fraction phase
diagram. The cross symbols show the gas andliquid phases. The black line represents thetieline, which is obtained from asymptotic behavior of median between two phases (square)
towards the vertex of binodal line. We fit aparabolic function to the cross symbols withvolume fraction and tieline as two axis, anddetermine the critical point (circle).

adopt a rather short-range interaction. Morespecifically, we use the modified Lennard-Jones
potential with cutoff length rc = 1.7σ: U (r) =4((σ/r)24 − (σ/r)12 + c2 (r/σ)2 + c0 ) for r < rc
and U (r) = 0 for otherwise. Here c2 and c0 are
constants to satisfy U (rc ) = 0 and U 0 (rc ) = 0.
This short-range nature of intercolloid potential suppresses the correlation length of density fluctuations and allows us to perform simulations with a rather small system size (L =17.3σ).

275

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

Figure 2: The correlation length ξ (a) and susceptibility χ (b) as a function of critical parameter  := Tc /(T −Tc ). The red line represents a
power-law function with critical exponent (approximately 0.62 and 1.24 for ξ and χ, respectively)

correlation function decays exponentially witha single time τξ scale, and τξ is proportional
to the cube of ξ, which coincides with the caseof dynamical critical phenomena for simple liquids.
In summary, we performed numerical simulations for colloidal suspensions with and without HI and examined the effect of HI on the dynamical critical phenomena. We obtained theresults that HI does not influence both staticand dynamical behaviors. However, due to theconstraint by numerical cost, the critical parameter that we accessed is limited (0.05 < ),and the anomaly in the time correlation function of concentration fluctuations is reported
in a smaller . To make a definitive conclusion
for this phenomenon, we realize a much larger
simulation.

References[1] J. Rouch, et al.: Phys. Rev. Lett. 71,(1993) 1947; B. M. Fine, et al.: Phys. Rev.
Lett. 74, (1995) 198; B. M. Fine, et al.:J. Chem. Phys., 104, (1996) 326[2] H. Tanaka and T. Araki: Phys. Rev. Lett.
85, (1999) 1338

We study the static behavior of colloidal suspensions near a critical point. We perform
simulations for states on tieline (see the blackline in Fig. 1) approaching to a critical point(red circle point). Figure 2 shows the correlation length ξ and susceptibility χ, which is determined by fitting the static structure factor
with the Ornstein-Zernike form. Here we cansee that the two results obtained from BD and
FPD show an excellent agreement with eachother, and both follow to power-law function
with the critical exponent. These trends are
consistent with a well-known fact that HI does
not influence the thermodynamic behaviors of
colloidal suspensions. Next, we compute theauto-correlation function of density fluctuationand evaluate the dependence of lifetime on thecritical parameter . We find that the auto-

276

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

Coarsening mechanism of network-forming phaseseparation in colloidal suspensions
Michio Tateno and Hajime Tanaka
Department of Fundamental Engineering, Institute of Industrial Science, University of Tokyo
4-6-1 Komaba, Meguro-ku, Tokyo 153-8505

Phase separation is one of the most fundamental physical phenomena in nature to form aspatially inhomogeneous structure, whose basic coarsening laws were well-established in
the 20th century. However, the conventional
mechanisms are valid only near the criticalpoint, and the coarsening behavior under deep
quench remains elusive. We noticed that for
network-forming deeply-quenched phase separation, unusual coarsening behavior of network
structures had been reported in experimentsfor various soft matter such as colloidal suspensions [1], protein solutions [2] and lyotropic
liquid crystal [3]: their characteristic network
size ` commonly grows as ` ∼ t1/2 while retaining their connectivity. The goal of this projectis to reveal the physical principle behind thegrowth exponent 1/2.
To address this problem, we used a hydrodynamic simulation model for colloidal
suspensions, Fluid Particle Dynamics (FPD)
method [5], which is based on the direct computation of the Navier-Stokes equation, andstudied network-forming process during gasliquid phase demixing. To capture the coarsening of network structure without suffering froma finite-size effect, we developed GPU-MPI hybrid codes utilizing a service provided by ISSP
and realized a large scale simulation. Since ourinterest is in phase separation dynamics undera deep quench, we neglect thermal noise.
We show in Fig. 1 the time evolution ofphase-separation structures in the colloidal
suspension. Here we can clearly see that aspace-spanning network structure is spontaneously formed in the early stage, and its characteristic length scale grows afterward. By
computing the temporal change of the first moment of the structure factor, we confirmed thatthe characteristic length of network structure `follows a power-law growth with exponent 1/2
while retaining self-similarity.
In the absence of thermal noise at zero temperature, the coarsening cannot be due to thermal activation but should be of purely mechanical nature. Then, we focus on the slow elastic
dynamics of dense colloidal aggregates and analyze the strain field  by coarse-graining thedisplacements of colloidal particles. Here inparticular, we consider the time evolution of
the distribution function, P (), from a reference time t0 to an arbitrary time t0 + t0 . We
find that the distribution functions are collapsed onto a master curve after scaling theargument  by t0 /t0 (see Fig. 2), indicating thatthe elastic relaxation of network structure controls the speed to the coarsening process.
Here we note that the exponent 1/2 is observed only in the presence of solvent component: for example, the Brownian dynamics
method completely neglecting the motion of asolvent does not reproduce the exponent 1/2(see, e.g., Ref. [4]). This fact implies that dynamic coupling between densely-packed colloid
and a solvent saturated in the colloid-rich domain plays a crucial role. Biot’s poroelastic
theory is known as a model to describe themotion of fluid saturated elastic materials. By

277

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

Figure 1: Coarsening behavior of network-forming phase separation of colloidal suspensions.
Particles are colored to distinguish front particles from back ones.
and the domain grows with the growth exponent of 1/2. We revealed that the exponent is
a consequence of slow elastic relaxation of the
network-forming dense phase, whose limiting
process is permeation flow of the solvent.

References[1] H. Tanaka, et al., J. Phys.: Condens.
Matter, 17, (2005) L143: A. Bailey, etal.: Phys. Rev. Lett. 99, (2007) 205701:
M. Tateno and H. Tanaka, npj Comput.
Mater., 5, (2019) 40

Figure 2: P () after scaling  by t0 /t0 . The
data are sampled from the data whose reference time is t0 = 31.0 (blue), 46.5 (green),
62.0 (brown).
Then, cross, triangle andsquare symbols represent the data at t0 /t0 =0.001, 0.01, 0.02, respectively.

performing a scaling analysis on this theory, we
obtain a coarsening law, ` ∼ (Dp t)1/2 , where
Dp is a poroelastic diffusibility. See Ref. [6] forfurther discussion on the validity of the coarsening law.

[2] H. Tanaka and Y. Nishikawa: Phys. Rev.
Lett., 95, (2005) 078103[3] Y. Iwashita and H. Tanaka: Nat. Phys.,
5, (2006) 147[4] A. Furukawa and H. Tanaka: Phys. Rev.
Lett. 104, (2010) 245702.
[5] H. Tanaka and T. Araki: Phys. Rev. Lett.
85, (1999) 1338
[6] M. Tateno and H. Tanaka: Nat. Commun., under review.

In summary, we performed a large-scale hydrodynamic simulation for colloidal suspensions and studied the dynamic process of
network-forming phase demixing. We foundthat the coarsening process has self-similarity,

278

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

Structural origin of glassy slow dynamics
Hua Tong and Hajime Tanaka
Department of Fundamental Engineering, Institute of Industrial Science, University of Tokyo
4-6-1 Komaba, Meguro-ku, Tokyo 153-8505

Despite the use of glasses for thousandsof years, the nature of glass and the glasstransition remains probably “the deepest andmost interesting unsolved problem” in condensed matter physics and materials science.
In contrast to crystallization where the solidity emerges as a result of the formation of
long-range periodic order, the drastic dynamical slowing down toward the glass transition isaccompanied by little change in the geometricstructure, as measured by two-point correlators accessible through diffraction and scattering experiments [1, 2]. Therefore, glasses arewidely thought as the epitome of completely
disordered state of materials and the glasstransition is described as a purely dynamical
phenomenon in kinetically constrained models.
The physical scenarios which posit a growing
static order and hence a thermodynamic origin behind the slowing down of glassy dynamics have regained popularity since the discovery of the so-called dynamic heterogeneity in
1990s. The spatially correlated domains whichmove significantly faster or slower than the average are proposed to be the long sought-after
cooperatively relaxing regions (CRRs), whichis the core concept of Adam-Gibbs theory ofglass transitionand its modern version, i.e., therandom first-order transition (RFOT) theory.
Meanwhile, a growing correlation length characterizing the extent of heterogeneous dynamics is also suggestive of a similarity between
glass transition and the critical phenomena [3–5]. However, since the dynamic heterogeneity
is only accessible through dynamic variables[6, 7], the crucial physical mechanism is still
not established, namely a quantitative characterization of the glassy structural order andits link to dynamics, which precludes a decisive279

underpinning to the (thermodynamic) natureof the glass transition.
To look for key structural features responsible for glassy dynamics, one typical approach is
to consider specific physical aspects of the local
atomic environment, e.g., free volume, potential energy, and symmetry. The free volume
approach [8, 9] and the inherent potential energy based on the potential energy landscape(PEL) formalism enjoyed early success showing a clear macroscopic correlation (that is,for globally averaged quantities) with dynamics, but not microscopically at a particle level[10–12]. Therefore, more efforts have been devoted to the identification of locally favored
structures (LFSs) based on symmetry considerations. For instance, icosahedral [13, 14],
crystal-like orders [3, 4, 15], or more complicated topological clusters [16] are identified
and suggested as the origin of slow dynamics
in different glass-forming liquids. However, thegeneral relevance of LFSs in glassy slow dynamics remains controversial due to its strong
and sensitive system dependence. This may
arise from the system-dependent nature of locally low free-energy configurations [4]. Another approach developed by Cubuk et al. recently is to define some structural quantity,“softness”, based on machine-learning methods [17, 18]. The strong correlation observed
between softness and structure relaxation suggests that important structural features arecaptured. However, softness is defined in ahigh-dimensional space with more than 100
structure functions, prohibiting a clear identification of glassy structural order. It is desirable to unveil the structural characteristics ofglass-forming liquids from a physically inspired
perspective.

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

Recently in ref. [12], we constructed a setof structural order parameters, namely Θ in
2D and Ω in 3D , to capture sterically favored structures with high local packing capability in hard-sphere-like glass formers. In
the absence of obvious density fluctuations beyond the particle scale, such sterically favored
structures provide more room for particle motions through better arrangement, and hence
higher correlational (or, vibrational) entropy
and lower free energy [4]. This is analogous tothe crystallization process of hard spheres. We
showed at a picture level that the mobility develops systematically following the structuralorder parameter field and confirmed a common
structural origin for both fast β and slow α dynamics.
In this study, we move a crucial step forward to study quantitatively the formation of

structural order and its relation to slow glassy
dynamics. We find that the degree of glassy
structural order follows a linear scaling law
with temperature, which initiates around theonset temperature Ton of the super-Arrhenius
dynamics and is ceased by the dynamical glasstransition Tg . This provides us with a direct
quantitative relation between structural orderand structure relaxation time τα , and therefore represents a key thermodynamic characteristic of the glass transition. More interestingly, based on a nonlocal excitation scenario,
we confirm that such an intimate structuredynamics correlation is valid even microscopically at a particle level. Therefore, our results may provide an essential piece in the microscopic theoretical description of the longstanding glass transition problem from a structural perspective.

[1] J. C. Dyre, Rev. Mod. Phys. 78, 953 (2006).
[2] K. Binder and W. Kob, Glassy materials anddisordered solids: An introduction to their statistical mechanics (World Scientific, 2011).
[3] H. Tanaka, T. Kawasaki, H. Shintani, andK. Watanabe, Nature Mater. 9, 324 (2010).
[4] H. Tanaka, Eur. Phys. J. E 35, 113 (2012).
[5] J. S. Langer, Phys. Rev. E 88, 012122 (2013).
[6] M. D. Ediger, Annu. Rev. Phys. Chem. 51, 99(2000).
[7] L. Berthier and G. Biroli, Rev. Mod. Phys.
83, 587 (2011).
[8] M. H. Cohen and G. S. Grest, Phys. Rev. B
20, 1077 (1979).
[9] F. W. Starr, S. Sastry, J. F. Douglas, and S. C.
Glotzer, Phys. Rev. Lett. 89, 125501 (2002).
[10] A. Widmer-Cooper and P. Harrowell, Phys.
Rev. Lett. 96, 185701 (2006).
[11] A. Widmer-Cooper and P. Harrowell, J. Non-

Cryst. Solids 352, 5098 (2006).
[12] H. Tong and H. Tanaka, Phys. Rev. X 8,011041 (2018).
[13] F. C. Frank, Proc. R. Soc. Lond. A 215, 43(1952).
[14] P. J. Steinhardt, D. R. Nelson, and
M. Ronchetti, Phys. Rev. B 28, 784 (1983).
[15] M. Leocmach, J. Russo, and H. Tanaka, J.
Chem. Phys. 138, 12A536 (2013).
[16] C. P. Royall and S. R. Williams, Phys. Rep.
560, 1 (2015).
[17] E. D. Cubuk, S. S. Schoenholz, J. M. Rieser,
B. D. Malone, J. Rottler, D. J. Durian,E. Kaxiras, and A. J. Liu, Phys. Rev. Lett.
114, 108001 (2015).
[18] S. S. Schoenholz, E. D. Cubuk, D. M. Sussman, E. Kaxiras, and A. J. Liu, Nature Phys.
12, 469 (2016).

280

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

Universality of bulk-edge correspondence bynumerical methods
Y. Hatsugai
Department of Physics, University of Tsukuba
1-1-1 Tennodai, Tsukuba 305-8571, Ibaraki, JapanOne of the key concepts of the topological phases is the bulk-edge correspondence, that relates
non-trivial topology of the bulk to the low energy physics near boundaries or impurities[1]. Nontrivial bulk topology is hidden for low energy physics of the bulk and most of the characteristic
and experimental features appear as generalized edge states living near the boundaries. Although
its historical example is a quantum phenomenon (quantum Hall eﬀect), its validity is much
more universal and the concept is applicable even for non-quantum phenomena such as photonic
crystals of classical electromagnetic ﬁeld and even classical mechanics of periodic systems. Usinga formal analogy among the various systems, one can deﬁne a Berry connection for the bulk
(even for the classical systems). Although there is no physical, at least, direct meaning forthe topological quantities based on the Berry connection, the edges states predicted by thebulk-edge correspondence are physical observables. Such a hidden Berry connection of thebulk and topological numbers support the edge states of various physical phenomena. This is
the universality of the bulk-edge correspondence, which we numerically try to establish in the
project.
The corner states of the higher order topological phases are such generalized edge states whichcharacterize the speciﬁc bulk topology reversely. We have shown such examples supported by
the quantized ZN (n ∈ Z) Berry phases as a bulk topological invariants for a quantum three
dimensional model [2] and classical mechanics on a Kagome lattice[3]. We have further proposeda higher order topological Mott insulator that hosts spin-charge separated corner states anddemonstrated its validity numerically[4].
We also performed a topological characterization of the J1 − J2 spin ladder for higher integer
spins using a quantized Berry phases[5]. It justiﬁes a generalized valence bond solid state as
the adiabatic limit of the ground state obtained numerically. Then sequential topological phasetransitions associated with the valence bond reconstruction are clearly demonstrated.
We further discussed several other topics related such as non hermitian topological phenomena[6,7, 8, 9] and obtained new insight for the adiabatic heuristic argument[10] that is a conceptual/historical background of topological characterization.

References[1] Yasuhiro Hatsugai, Phys. Rev. Lett. 71, 3697 (1993).
[2] Hiromu Araki, Tomonari Mizoguchi and Yasuhiro Hatsugai: Phys. Rev. Research 2,
012009(R) (2020).
[3] Hiromasa Wakao, Tsuneya Yoshida, Hiromu Araki, Tomonari Mizoguchi, and Yasuhiro
Hatsugai: Phys. Rev. 101, 094107 (2020).
[4] Koji Kudo, Tsuneya Yoshida, and Yasuhiro Hatsugai : Phys. Rev. Lett. 123, 196402 (2019).

281

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

Figure 1: (a) Schematic lattice structure of the three dimensional BBH model, (b) corner statesand (c) Z4 Berry phases (See. Ref.[2] )Figure 2: (a) J1 − J2 Heisenberg spin chain of integer S(≥ 1) and (b) Z2 Berry phase of S = 4(8 site system with periodic boundary condition) (See Ref. [5] )[5] Shota Fubasami, Tomonari Mizoguchi and Yasuhiro Hatsugai: Phys. Rev. B 100, 014438(2019).
[6] Tsuneya Yoshida, Koji Kudo, and Yasuhiro Hatsugai: Scientiﬁc Reports 9, 16895 (2019).
[7] Tsuneya Yoshida and Yasuhiro Hatsugai: Phys. Rev. B 100, 054109 (2019).
[8] Tsuneya Yoshida, Robert Peters, Norio Kawakami, Yasuhiro Hatsugai: arXiv:2002.11265[9] Tsuneya Yoshida, Tomonari Mizoguchi, and Yasuhiro Hatsugai: arXiv:1912.12022[10] Koji Kudo and Yasuhiro Hatsugai : arXiv:2004.00859.

282

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

Kinetics of phase transition and polyamorphism
Kazuhiro Fuchizaki, Kei Watanabe, and Hiroki Naruta
Department of Physics, Ehime University, Matsuyama 790-8577

In FY2019, we have undertaken two subjects, as reported brieﬂy below. The ﬁrst task
has been almost completed, and its full reportwill soon be submitted elsewhere. The second
one is still at a preliminary stage and will be
ongoing in FY2020.
Growth of density fluctuations at the beginningof a phase separation
The concept of spinodal has been examined
for an Ising system [1]. However, no explicit
discussion has been found for a ﬂuid system,
which is believed to belong to the same unitary class. Ushcats has suggested [3], extending the equation of state beyond the thermodynamically stable region, that the spinodal line of the modiﬁed Lennard-Jones (mLJ)
system [2] approaches the liquid-gas coexistence (binodal) line. Our prediction based on
the coarse-graining limit for the coarse-grained
free energy, estimated by applying the blockspin coarsening technique, supports Ushcats’
conjecture [4]. It is then highly anticipatedthat the “dynamics” of a system also depends
on a coarsening level.
We conducted isothermal–isobaric molecular
dynamics (MD) simulations for an mLJ system
consisting of 4 × 106 particles to examine thedynamical process of liquid–gas phase separation. Length, time, and temperature are ex√pressed in units of σ, mσ 2 /ϵ, and ϵ/kB [5],respectively, where σ and ϵ are the fundamental length and energy scales of the mLJ potential function and kB denotes Boltzmann’s constant. The critical temperature, pressure, anddensity of the mLJ system are Tc = 1.0762(2),pc = 0.09394(17), and ρc = 0.331(3), respectively. We treated two systems with diﬀer-

283

Figure 1: The peak position of SL is plotted against time with the same color for every
coarsening scale of L. Those at ρc (0.02) aredelineated by solid (dotted) lines.

ent density; 0.02 and ρc . The ﬂuid systemwell-equilibrated at 20Tc was instantaneously
brought to 0.7Tc to undergo phase separation.
The lower- and critical-density systems were
expected to be quenched in the metastable andunstable regions, respectively, judging fromthe shape of the coarse-grained free energy.
The term, “local density,” has been quite
often used in mesoscopic descriptions. However, it has been far unclear what is meantby “local.” In this case, we can deﬁne the local density unambiguously using the length L,
the coarse-graining scale based on which thecoarse-grained free energy was deﬁned [4].
Based on this well-deﬁned “local” density
and its ﬂuctuations around the average density, we deﬁned the structure factor SL . Figure 1 shows the time evolution of SL after
quenching. On the short-time side (t < 2),
L predominantly determines the peak positionof SL , whereas it approaches the value determined by the average density on the long-timeside (t > 2).

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

To examine in detail the L-dependence of
the wavelength at which structural destabilization occurs, SL was scaled by its peak
position and height. The dynamical process
in another state, obtained by quenching thecritical-density system to 0.97Tc , was also examined. We could identify in the short-time
region the new universal process (Figure 2), in
which phase separation proceeds in a statistically self-similar fashion irrespective of the average density, quenching temperature as wellas the coarsening length. We name this stage
a “super-early” stage.

assertion.
An isochoric–isothermal MD simulation generated an equilibrated mLJ ﬂuid consisting of2048 particles. We employed such a small number of particles because Goddard et al. used
for a similar size for an LJ ﬂuid [6]. Similarﬁnite-size eﬀects would then be involved in the
results. We chose for this purpose three stateswith density 5.00 × 10−2 , ρc , and 0.600 along
the isotherm 1.001Tc , where Tc is the aforementioned critical temperature of the mLJ system [6]. The density is again expressed interms of the mLJ units.
The entropy of the three states was also calculated using the λ-integration [7], which is too
costly but can give the numerically exact result. The entropies based on the method of
Goddard et al. could not well reproduce theλ-integration based ones. We should undertake a critical review of their basic idea.

References[1] K. Kaski, K. Binder, and J. D. Gunton,
Phys. Rev. B 29, 3996 (1984).
Figure 2: The scaled structure factors at three
diﬀerent states fall into a single curve whenplotted against scaled wavelength, implying
the existence of universality in a phase separation process right after a quench.

[2] Y. Asano and K. Fuchizaki, J. Chem.
Phys. 137, 174502 (2012).

Estimating entropy of a fluid
A method proposed by Goddard et al. [6] hasbeen widely used for evaluating entropy of adisordered system even with such a kind of network structure. The entropy is generally estimated from the density of states (DOS), which
is obtainable through Fourier inversion of thevelocity autocorrelation. Goddard et al. haveclaimed that a DOS of a ﬂuid state can be constructed from those of solid and gas states [6].
This claim does not have any physical basis.
Although an LJ ﬂuid was used to testify theirclaim [6], we employ as another simple liquid
an mLJ system to thoroughly examine their[5] K. Okamoto and K. Fuchizaki, J. Phys.
Soc. Jpn. 86, 034003 (2017).

284[3] M. V. Ushcats, J. Chem. Phys. 141,
101103 (2014).
[4] K. Fuchizaki and K. Watanabe, J. Phys.
Soc. Jpn. 87, 114006 (2018).

[6] S. Lin, M. Blanco, and W. A. Goddard,
J. Chem. Phys. 119, 11792 (2003).
[7] J. Q. Broughton and G. H. Gilmer, J.
Chem. Phys. 79, 5095 (1983).

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

Ground-State Phase Diagram of the S=1/2 Kitaev-G model ona Honeycomb Lattice
Takafumi SUZUKI
Graduate School of Engineering,
University of Hyogo, Shosha, Himeji, Hyogo 670-2280
We have studied the ground-state phase

Next,

we

have

studiedlow-energy

diagram of the S=1/2 Kitaev-G model on a

excitations of an S=1 Ising spin chain with thehoneycomb lattice [1] and quantized DS=2

negative single-ion anisotropy in magnetic

excitation spectra by confinement in an S=1

fields. We have calculated the dynamical

spin chain [2].

structure factor with the infinite time-evolving
bock decimation (iTEBD) method [4]. We haveFirst, we have investigated the ground-statefound that when the transverse magnetic field isphase diagram of the S=1/2 Kitaev-G model on

absent, both a DS=2 excitation continuum anda honeycomb lattice with the numerical exact

one-magnon isolated mode appear in the low-diagonalization method, DMRG and cluster-

lying

series expansion method [3]. By changing theexcitation.

continuum

magnitude of interactions on each bond, the

is

This

DS=2

excitation

quantized

by

the

weak

longitudinal magnetic field. The quantized DSmodel connects with three characteristic points,=2 excitation spectra originate from thewhere the model is described by spin chains,confinement of two domain walls and theisolated dimers, or an isotropic honeycomb-

excitation energies of the quantized DS=2

lattice model. The obtained results havespectra are well explained the negative zeros inclarified that in the spin chain limit, two

the Airy function.

Tomonaga-Luttinger (TL) liquid phases andReferences

three non-TL liquid phases appear when the[1] Y. Yamada, et al., in preparation.

ratio of the Kitaev interaction and the G[2] T. Suzuki and S. Suga, J. Phys. Soc. Jpn. 88

interaction changes. In the dimer limit, the(2019) 053702.

dimerized state with the lowest energy is |tx⟩ =[3] J. Oitmaa, et al., Series Expansion Methods(| ↑↑⟩−i| ↓↓⟩)/ 2, when the positive G interactionforand the negative Kitaev interaction are

Strongly

Interacting

LatticeModels,Cambridge University Press, 2006.

considered. We have found that this |tx⟩ state[4] H. N. Phien, et al., Phys. Rev. B 86 (2012)

can survive up to the isotropic point.

245107.

285

Figure 2 shows a two-component fluid under thermal conduction separated by a semipermeable membrane. The osmotic pressure
is determined from the balance of the chemical potentials in the isothermal state, but the286
sha1_base64="nrZYJsaR2ltZo5JYKWOP8wCMXg4=">AAAC+nichVJNLwNRFD0d3/VVbCQ2jYZYNbdIiFXDRqwoRYI0M+PhpfNlZtqkyh+wsUSsSCwaW//Axh+w8BPERkJiY+HOdBJB6J3Me/eee899751czTGk5xM9xpSm5pbWtvaOeGdXd09voq9/1bNLri7yum3Y7rqmesKQlsj70jfEuuMK1dQMsaYV54L8Wlm4nrStFb/iiC1T3bXkjtRVn6H8QqG6f1RIpChNoSV/O5nISSGyRTvxgk1sw4aOEkwIWPDZN6DC428DGRAcxrZQZcxlT4Z5gSPEmVviKsEVKqNFXnc52ohQi+OgpxeydT7F4N9lZhIj9EA1eqV7uqEn+vizVzXsEdylwrtW5wqn0Hs8uPzekGXy7mPvi/UPQ+ObqRwLZvz9tqDObJBvVCE47/CtBPYbqOhjB9OhepLVdEIk0FWvv6h8cPq6PJMbqY7SFT2zopf0SHesqVV+06+XRO4CcR6JzM8B+O2sjqczE2lamkxlZ6PhaMcQhjHGEzCFLOaxiDyfK3GCM5wrh8qVUlNu6qVKLOIM4Jspt58wQKq1</latexit><latexit

sha1_base64="Xe30dIM3EgGsd+RFucw669WYRvg=">AAADCXichVJNSxxBEH2OidFNohu9CF4ki5JLlloTMHha1INHv1YFV5aesdYMzpczvYvrsH/A/IAcBCGBEENOuYlXL/kDOfgTQi4BA148WDM7EBLRrWG6q17Vq+5+lBk4dqSJLnqM3gcP+x71D+QeP3k6OJR/NrwW+Y3Q4orlO364YaqIHdvjira1wxtByMo1HV43d+eS/HqTw8j2vVXdCnjLVTueXbctpQWq5YdXX1brobLi6jw7WrXjqXYtX6AipTZ+2yllTgGZLfr536hiGz4sNOCC4UGL70Ahkm8TJRACwbYQCxaKZ6d5Rhs54TakiqVCCbor645EmxnqSZz0jFK2Jac48ofCHMcE/aAvdEnf6Sv9pOs7e8Vpj+QuLdnNDpeD2tDh6MpVV5Yru8bbv6x7GKbcTEnMwrj7bUmd2yXfrYIlH8itGHtdVNSo402qni1qBimS6Gp1XtQ8eH+5MrM8EU/SR/olin6gCzoXTb3mH+vTEi8fIScjUfp/AG47a1PF0qsiLb0ulGez4ejHGJ7jhUzANMpYwCIqcu4+jvEZJ8Y745txapx1So2ejDOCf8w4vwFIZLA3</latexit>

Figure 1 shows a one-dimensional system ofparticles under thermal conduction that interact with each other at the Leonard-Jones potential, and each particle is subjected to a constant external force (gravity). Various local
thermodynamic quantities such as temperature, density, pressure, and chemical potential
were determined by molecular dynamics simulations, and from these quantities the free energy and entropy of the whole system were estimated. The thermodynamic relation is foundto give the free energy change of the whole system when the temperature equivalent to thetotal kinetic energy of the system and the spatial average pressure are adopted as temperature and pressure for the system, respectively.
The equation of state has an additional effectfrom the second order of the external force,
implying a limit to the global view of the thermodynamic system under the external force.
<latexit

Tsha1_base64="469r/arnv6Hx41C0Ys0aEcx8564=">AAADCXichVJNSxxBEH2OidFNohu9CF4kixIILLUmYPC0qAePfq0Kriw9Y60ZnC9nehfXYf+A+QE5CEICIYacchOvXvIHcvAnhFwCBrx4sGZ2ICSiW8N0V72qV939KDNw7EgTXfQYvQ8e9j3qH8g9fvJ0cCj/bHgt8huhxRXLd/xww1QRO7bHFW1rhzeCkJVrOrxu7s4l+fUmh5Hte6u6FfCWq3Y8u25bSgtUyw+vvqzWQ2XF1Xl2tGrHU+1avkBFSm38tlPKnAIyW/Tzv1HFNnxYaMAFw4MW34FCJN8mSiAEgm0hFiwUz07zjDZywm1IFUuFEnRX1h2JNjPUkzjpGaVsS05x5A+FOY4J+kFf6JK+01f6Sdd39orTHsldWrKbHS4HtaHD0ZWrrixXdo23f1n3MEy5mZKYhXH325I6t0u+WwVLPpBbMfa6qKhRx5tUPVvUDFIk0dXqvKh58P5yZWZ5Ip6kj/RLFP1AF3QumnrNP9anJV4+Qk5GovT/ANx21qaKpVdFWnpdKM9mw9GPMTzHC5mAaZSxgEVU5Nx9HOMzTox3xjfj1DjrlBo9GWcE/5hxfgNDILA1</latexit>

2<latexit

sha1_base64="vswUlqY0nPTdcJkB/sWm+IVFr14=">AAAC9nichVLLLgRRED3aa4zXYCOxERNiNalBQqyEjeUMBgki3a3Q0S/dPRNM/ABbYiGRkIiIjX+w8QMWPkFsJCQ2Fqp7OhGEqU7fW3WqTt17T0pzTcMPiB5rlNq6+obGRFOyuaW1rT3V0TnvO0VP54LumI63qKk+m4bNhcAITF50PVYtzeQFbWsqzC+U2PMNx54Ldl1esdQN21g3dDUQKD+3mkpThiLr/e1kYyeN2HJO6gXLWIMDHUVYYNgIxDehwpdvCVkQXMFWUBbME8+I8ox9JIVblCqWClXQLVk3JFqKUVvisKcfsXU5xZTfE2Yv+umBrumV7umGnujjz17lqEd4l13ZtQqX3dX2g+7Z96osS/YAm1+sfxia3EyVmIXx99vCOqtKvloFS96VWzG2q6gYYB1jkXqGqOlGSKirXnlRae/kdXZ8pr88QBf0LIqe0yPdiaZ26U2/zPPMKZIyEtmfA/DbmR/KZIczlB9JT0zGw5FAD/owKBMwiglMI4eCnMs4xBGOlR3lTLlUriqlSk3M6cI3U24/ARoZqM8=</latexit>

latexit

中間温度<

Figure 1: Heat conduction system under an
本題：これに定常熱流を流す
external force.
Tsha1_base64="8k6rR3Vgsf1IwFoGWlzCGMF0PQ0=">AAAC+3ichVK7ThtRED0s4RFDwECDRGNhOaKyxoAEokJAQQk4NkgGod3NQBb2xe61JbD4ASo6BFRJRMGj5Qto8gMp+AREgwQSDUVm1ytFARHPau+dOTNn7r1HY/i2FSqi2xat9UNbe0fnx1RX96ee3nRffzn0qoHJJdOzvWDF0EO2LZdLylI2r/gB645h87KxPRvll2schJbnflG7Pq85+qZrbVimrgQqr86xrfT1dJbyFFvmrVNInCwSW/DSD1jFV3gwUYUDhgslvg0doXwVFEDwBVtDXbBAPCvOM/aREm5VqlgqdEG3Zd2UqJKgrsRRzzBmm3KKLX8gzAxy9JvO6ZF+0SXd0cu7vepxj+guu7IbDS77670Hg8XnpixHdoVvf1n/YRhyM11iFsb7b4vqnCb5ZhUseV9uxdhpoqLCBiZj9SxR04+RSFez8aLa3tFjcWopV/9MP+heFP1Ot3Qjmrq1J/NskZdOkZKRKLwegLdOeTRfGMvT4nh2eiYZjk4MYRgjMgETmMY8FlCSc7dwiGOcaPvaT+1Cu2qUai0JZwD/mHb9B92NqvM=</latexit>

熱浴温度差<latexit

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

Macroscopic properties characterized by an extended
thermodynamic functions to nonequilibrium
Naoko NAKAGAWA
Department of Physics, Ibaraki University
2-1-1 Bunkyo, Mito, Ibaraki 310-8512

Thermodynamic functions such as entropy
and free energy play an essential role in equilibrium systems. In order to extend these powerful functions to various non-equilibrium systems, such as thermal conduction systems andexternal force-driven systems, we have performed numerical experiments on the systems
shown in Figs. 1 and 2.

T+2

定常熱流 Jq

Figure 2: Two component fluid separated by asemi-permeable membrane.
局所平衡が成り立つことを仮定して系の定常状態を調べる
11

chemical potentials are not spatially uniform in
the thermal conduction state. In order to explore the principle of determining the osmoticpressure under thermal conduction, molecular
dynamics simulations are performed and thedependence of the osmotic pressure on various parameters is investigated. We used the
LAMMPS in order to achieve efficient results.

References[1] N. Nakagawa and S.-i. Sasa: J. Stat. Phys.
177 (2019) 825.

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

Development of data-driven science tomaterials science
Koji Hukushima
Department of Basic Science, University of Tokyo
3-8-1 Komaba, Meguro-ku, Tokyo 153-8902

In recent years, data-driven science methods
have been attracting attention and have been
applied to condensed matter physics and materials sciences. There are two main purposes toanalyze the data obtained as a result of experiments and numerical experiments. One would
be to extract the mathematical structure behind the data, and the other would be to makepredictions to get the next data.
In this project, we carried out the tasks corresponding to the above-mentioned two purposes as actual problems. As the former problem, the model Hamiltonian was estimatedfrom experimental data with a magnetic compound. As an example of the latter problem, we have conducted research on datadriven scientific methodologies in crystal structure search.

interactions or physical quantities are also discussed. These are the advantages of samplingfrom the posterior distribution. The obtained
eﬀective model is useful to predict hard-tomeasure properties such as spin gap, spin configuration at the ground state, magnetic specific heat, and magnetic entropy.
Next, we studied the eﬃciency of Bayesian
optimization as a method of selecting the next
candidate when searching for a crystal structure. Bayesian optimization is a global optimization machine learning technique for a
black-box function that requires a large amountof computation to be evaluated. A crystal
structure is represented by a numerical vector called a descriptor that is used to quantify how similar structures are. This descriptor is an input of the crystal structure searchwith Bayesian optimization. The eﬃciency ofthe optimization significantly depends on the
choice of the descriptor. We evaluated its efficiency with a varying parameter value in thedescriptor[1]. Applying the crystal structure
search to crystalline silicon shows that the efficiency of the search depends heavily on theparameter value. We find that the eﬃciency
is linked to the distribution of the descriptor.
Therefore, we introduce an information measure of the distribution to estimate an appropriate parameter value for performing the crystal structure search eﬃciently. The measure
can also be used to predetermine an appropriate parameter value. The validity of the measure is confirmed with its applications to siliconoxide and yttrium cobalt alloy.

First, we proposed a data-driven methods toextract the spin Hamiltonian with uncertainty
from multiple experiments. Using our method,an eﬀective model of KCu4 P3 O12 is determined from the experimentally observed magnetic susceptibility and magnetization curves
with various temperatures under high magnetic fields. An eﬀective model is the quantum Heisenberg model on a zigzag chain witheight spins in which four independent parameters are contained. Based on Bayesian statistics, these four parameters were estimated using MCMC sampling from the posterior distribution with experimental data as input.
Our method yields that J1 = 8.54 ± 0.51 meV,J2 = −2.67 ± 1.13 meV, J3 = −3.90 ± 0.15
meV, and J4 = 6.24±0.95 meV, describes thesemeasured results well. Evaluation of these error bars is possible only by evaluating noise. [1] N. Sato, et al: Phys. Rev. Materials 4,
The relations among the estimated magnetic (2020) 033801.

287

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

Uniﬁed understanding of early-time nonequilibrium
relaxation in Monte Carlo simulations
Yoshihiko NONOMURA
MANA, National Institute for Materials Science
Namiki 1-1, Tsukuba, Ibaraki 305-0044

χ(t, T ) ∼ (T −Tc )−γ fsc [ctσ +ln(T −Tc )γ ], (1)which means that all the oﬀ-critical relaxation
data near Tc are scaled on a single curve.
Here we assume that the relaxation process
along the Kibble-Zurek mechanism [3], namely
the controlled annealing process in accordancewith critical relaxation, is expressed by theabove scaling curve, and try to reconstruct theannealing process from the NER data in thetwo-dimensional (2D) Ising model. Althoughvarious attempts within the standard powerlaw annealing are not successful, the following
stretched-exponential modiﬁed KZ process,T (t) = T0 + Td exp[−C(t/τ )σ ],(2)
Tf − aTiTi − Tf
T0 =, Td =; a = e−C , (3)1−a1−agives consistent results for L = 4096 with the106

σ = 1/5
σ = 1/3
σ = 1/2
L = 4096

5

10

104

χ(t)

Recently we numerically revealed that thecritical nonequilibrium relaxation (NER) behaviors in cluster algorithms are described by
the stretched-exponential simulation-time dependence of physical quantities in various spin
systems [1], and we derived such a relaxation
formula phenomenologically [2].
Quite recently, we found a general derivationof a scaling for oﬀ-critical relaxation in localand cluster-update algorithms. As an example, we consider the magnetic susceptibility incluster-update algorithms. From the initialtime critical relaxation χ(t; Tc ) ∼ exp(ctσ )and the temperature dependence in equilibrium χ(t = ∞, T ) ∼ (T − Tc )−γ , we haveχ(t, T )(T − Tc )γ ∼ exp[ctσ + ln(T − Tc )γ ], orreconstructedfrom NER data

C=2

KZ simulation

103
102
101
100

0

20

40

60

80

100

t [MCS]
Figure 1: Simulation-time dependence of themagnetic susceptibility in the 2D Ising modelin the SW algorithm for the modiﬁed KZ process explained in the main text. The dashedand solid lines are the simulated and reconstructed data based on Eq. (1), respectively.
relaxation process in the Swendsen-Wang algorithm with a set of parameters for a rather
slow relaxation: Ti = 1.5Tc , Tf = Tc , τ = 100,
C = 2, and σ = 1/5, 1/3, 1/2 with the exact
critical temperature Tc as shown in Fig. 1.

References[1] Y. Nonomura, JPSJ 83, 113001 (2014);
YN & Y. Tomita, PRE 92, 062121 (2015),
93, 012101 (2016), 101, 032105 (2020).
[2] YT & YN, PRE 98, 052110 (2018).
[3] C.-W. Liu et al., PRB 89, 054307 (2014);
Refs. therein (including originals by KZ).

288

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

Generalization of temperature scaling in early-timenonequilibrium relaxation
Yoshihiko NONOMURA
MANA, National Institute for Materials Science
Namiki 1-1, Tsukuba, Ibaraki 305-0044

χ(t, T ) ∼ (T −Tc )−γ fsc [ctσ +ln(T −Tc )γ ]. (1)We conﬁrm this temperature scaling by analyzing the NER process of the three-dimensional
Heisenberg model in the Swendsen-Wang algorithm from the perfectly-disordered states [3].
Here we take the L × L × L system with theperiodic boundary condition for L = 560, thelargest size used in Ref. [3]. The data up to 225
Monte Carlo steps are scaled with Eq. (1) in
Fig. 1. Minimizing the mutual residuals using
our previous estimates with the NE-ES, Tc =1.442987(2)J/kB and σ = 0.47(1) [3], we have
γ = 1.3945 ± 0.0019, c = 1.2593 ± 0.0043, (2)which is consistent with the previous estimate
from detailed simulations, γ = 1.3960(9) [6].

100γT = 1.4436 J/kB
T = 1.4438 J/kB
T = 1.4440 J/kB
T = 1.4442 J/kB
T = 1.4445 J/kB
T = 1.4450 J/kB
T = 1.4460 J/kB
T = 1.4480 J/kB
T = 1.4500 J/kB
T = 1.4520 J/kB

χ(t,T)(T-Tc)

Recently we numerically revealed that thecritical nonequilibrium relaxation (NER) behaviors in cluster algorithms are described by
the stretched-exponential simulation-time dependence of physical quantities in various classical spin systems [1, 2, 3] and in a quantum
phase transition [4], and derived such a relaxation formula phenomenologically [5]. There
we utilized the nonequilibrium-to-equilibrium
scaling (NE-ES) which connects initial-timeand equilibrium relaxation behaviors.
Here we generalize this derivation to oﬀcritical behaviors. As an example, we consider
the magnetic susceptibility. From the initialtime critical relaxation χ(t; Tc ) ∼ exp(ctσ )and the temperature dependence in equilibrium χ(t = ∞, T ) ∼ (T − Tc )−γ , we haveχ(t, T )(T − Tc )γ ∼ exp[ctσ + ln(T − Tc )γ ], or

10-1

10-2

10-3

10-4
-10

-8

-6

-4

-2

0

2

4

ctσ+ ln(T-Tc)γ

6

8

10

Figure 1: Temperature scaling plot of the magnetic susceptibility in the 3D Heisenberg model
for various temperatures in a semi-log scale using the estimates of Tc and σ in Ref. [3].

References[1] Y. Nonomura, J. Phys. Soc. Jpn. 83,
113001 (2014).
[2] Y. Nonomura and Y. Tomita, Phys. Rev.
E 92, 062121 (2015).
[3] Y. Nonomura and Y. Tomita, Phys. Rev.
E 93, 012101 (2016).
[4] Y. Nonomura and Y. Tomita, Phys. Rev.
E 101, 032105 (2020).
[5] Y. Tomita and Y. Nonomura, Phys. Rev.
E 98, 052110 (2018).
[6] M. Campostrini, M. Hasenbusch, A. Pelissetto, P. Rossi, and E. Vicari, Phys. Rev
B 65, 144520 (2002).

289

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

Finite-temperature properties of theKitaev-Heisenberg models on kagome and triangularlattices
Takami TOHYAMA
Department of Applied Physics, Tokyo University of Science, Tokyo 125-8585

Frustrated quantum spin systems such as
the Heisenberg and Kitaev models on various
lattices have been known to exhibit various
exotic properties not only at zero temperature but also for ﬁnite temperatures. Inspiredby the remarkable development of the quantum frustrated spin systems in recent years, we
decided to investigate the ﬁnite-temperature
properties of the S = 1/2 Kitaev-Heisenberg
models on kagome and triangular lattices.
For calculating accurate low-temperature
properties, we developed two types of improved version for ﬁnite-temperature Lanczos
method (FTLM) [1]: replaced FTLM and orthogonalized FTLM. In standard FTLM [2], avery large number of samplings of initial vector is required to obtain good accuracy at low
temperatures. The low-temperature FTLM [2]is known as one of the solutions to this problem. However, this method has a diﬃculty
for large-scale calculations because it requires
huge random access memory to keep all vectors in the Krylov subspace. In RFTLM, lowlying eigenstates obtained by the ground-state
Lancoz method are used in the patrition function and other thermodynamic quantities. In
OFTLM, initial randome vectors are prepared
to be orthogonal to the low-lying eigenstates.
We also developed a Lanczos code that equips
a technique to save the memory for the Hamiltonian by dividing H into two subsystems.
These techniques lead us to calculate ﬁnitetemperature proprieties up to 36 lattice sites
with good accuracy even at low temperatures.

The newly developed RFTLM was also usedin the calculation of speciﬁc heat and magnetic susceptibility for S = 1/2 frustrated J1 J2 Heisenberg model on kagome, triangular,
and square lattice up to 36 sites [3].
In both kagome and triangular lattices for
the Ktaev-Heisenberg model, multiple peaks
are conﬁrmed in the speciﬁc heat. To ﬁnd
the origin of the multiple peaks, we calculate
the static spin structure factor. The origin of
the high-temperature peak of the speciﬁc heatis attributed to a crossover from the paramagnetic state to a short-range ordered state
whose static spin structure factor has zigzag
or linear intensity distributions in momentum
space. In the triangular Kitaev model, the order by disorder due to quantum ﬂuctuation occurs. On the other hand, in the kagome Kitaev
model it does not occur even with both quantum and thermal ﬂuctuations [1].

References[1] K. Morita and T. Tohyama: Phys. Rev.
Research 2 (2020) 013205.
[2] P. Prelovsek and J. Bonca: Ground stateand ﬁnite-temperature Lanczos methods,
in Strongly Correlated Systems, Numerical Methods, Vol. 176 of Springer Series in
Solid-State Sciences (Springer, 2013).
[3] P. Prelovsek, K. Morita, T. Tohyama,
and J. Herbrych: Phys. Rev. Research 2(2020) 023024.

290

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

Constructing the electronic structure database for themolecular design of pigments utilizing near infrared light
Yu KOMATSU
National Institutes of Natural Sciences Astrobiology Center
National Astronomical Observatory of Japan
2-21-1 Osawa, Mitaka, Tokyo 181-8588
NIR

are shown with 12 kinds of metals in fourradiation are being anticipated for dye-sensitizedsolvents. The data will be combined with ansolar

existing dataset like QM9 for further analysis to

Novel

dye

cells

molecules
andabsorbingimaging.

Moreover,search for functional molecules effectively.

cyanobacteria, oxygenic photosynthetic bacteria,
turn to utilize chlorophyll f, which is the reddest

The computation in related work has beenchlorophyll, for photochemistry after growing

done using supercomputers of Tokyo Institute ofunder the far-red light [1]. Thus, the red limit ofTechnology

photosynthesis has been updated. In this study,

Observatory of Japan, as well.

andNational

Astronomicalin order to explore redder molecules withpossibly working as photosynthetic pigments,
the database has been constructed using ab initio
calculations for ring molecules like porphyrins
with certain properties, combined with highthroughput methods.
Thephysical-chemicalquantitiesofmolecules, were estimated at DFT level asfollows: [total energies in ground and excited

Fig. 1: Estimated ionization potentials andstates, orbital energies, ionization potentials

electronic affinities for (a) bacteriochlorophylls(IPs), electronic affinities (EAs), bindingb,energies for central metals, electronic structures

dibenzoporphycenes with different central

for each orbital, …]. Dye molecules aremetals in 4 kinds of solvent conditions (CAM-

calculated with different central metals (Mg, Ca,

B3LYP/Def2tzvp//B3LYP-D3/Def2tzvp).

Ni, Zn, Sr, Pd, Cd, Ba, Pt, Hg, Pb and H2） in

References

vacuum and water, methanol and benzene with

[1] D. J. Nürnberg et al.,: Science 360 (2018)

PCM. In Fig.1, IPs and EAs for three pigments

1210.

291

(b)

phthalocyanines

and

(c)

meso-

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

Application of numerical renormalization groupmethod in condensed- matter physics
Kenji HARADA
Graduate School of Informatics, Kyoto University, Kyoto 606-8501, JapantensorIn the case of the (2+1)-dim. tensor

renormalization group (TRG) methods [1-4]network, the performance of the TRG

to various condensed matter systems. We

method is not stable. The average density

tried to use them for the non-equilibriumof active objects is not consistent with thesystems.

Monte Carlo results, not near a critical

We

have

appliedthepoint.

The directed percolation (DP) show aTheaccuracyoftheobliquecontinuous phase transition in the phaseprojection in the TRG method is bad.

diagram in the long-time limit. We have

However, due to the higher complexity ofstudied the phase diagram through thethe TRG method, we could not increase theRenyi entropy by using a TEBD method on

bond-dimension of tensors even using thethe matrix product states representation ofsystem B. The reduction of the complexity

a state-probability distribution [5]. In theof the TRG method is a future work tocase of one-dimensional DP, we have foundimprove the accuracy of the TRG method.

a critical relaxation at the critical point

Reference

and a new phase transition point in the[1] M. Levin and C. P. Nave, Phys. Rev.

active phase of the one-dim. DP.

Lett. 99, 120601 (2007).

The state-probability distribution of thed-dim. DP can be written as a (d+1)-dim.

[2] Z. Y. Xie, J. Chen, M. P. Qin, J. W.

tensor network. In the case of one-dim. DP,

Zhu, L. P. Yang, and T. Xiang, Phys. Rev.

the TRG shows better accuracy than the

B 86, 045139 (2012).
[3] G. Evenbly and G. Vidal, Phys. Rev.

TEBD method.

Lett. 115, 180405 (2015).

To check the existence of the new[4] S. Yang, Z.-C. Gu, and X.-G. Wen,transition point in the activity phase of

Phys. Rev. Lett. 118, 110504 (2017).

higher dim. DP, we tried to apply the TRG[5] K. Harada and N. Kawashima,method on the two-dim. DP by a massively

Phys. Rev. Lett., 123, 090601, (2019).

parallel computer on the ISSP.

292

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

Clarification of thermal transport spectra in hierarchicalstructured organic bulk materials
Takuma Shiga
Department of Mechanical Engineering,
The University of Tokyo, Hongo, Bunkyo, Tokyo 113-8656
Material design of cellulose nanostructureThis force field is composed of bond stretching,has been recently paid much attention becausebond angle, dihedral angle, van der Waals, andcellulose structure was identified by X-ray andCoulombic interaction. For non-bonding, we

neutron scattering spectroscopies. Cellulose

used the combination of particle-particle particle

nanomaterial is expected as candidates such as

meth and cutoff methods [2]. In the calculation,transparent

thermalwe applied periodic boundary conditions along

insulator since it has high degree of self-

all directions. After equilibration with T=300 K,

organization and high mechanical property.

hot and cold thermostats with T=330 K and

Turning

thermalT=270 K were attached to the simulation cell as

conductivity of cellulose thermal insulators suchillustrated in Fig. 1. Heat flux and temperature

as aerogel and foam is around 15-18 mWm-1K-1,

gradient were evaluated after the system reaches

and it is lower than widely-used glass wool andnon-equilibrium steady-state, and then thermal

urethane foam. However investigation of

conductivity was evaluated through the Fourier's

microscopic heat conduction mechanism oflaw. All calculations were done by LAMMPS

cellulose nanomaterial is very limited. In order

package.

totofurther

conductive

thermalimprovefilm

andinsulation,

thermalinsulation

performance, it is necessary to get insight into
the knowledge on cellulose nanocrystal andfibrils.
We have evaluated heat conduction and its
size effect of I-b type cellulose nanocrystal andFig. 1 Schematics for (a)cellulose nanocrystal

single chain by performing non-equilibrium(CNC) and (b)single chain cellulose.
Size effect of thermal conductivity ofmolecular dynamics (NEMD) simulation. In thethecellulose nanocrystal is shown in Fig. 2. As

GLYCAM06 [1] force field which reproduces

thermal conductivity increases with long-axis

thermal expansion coefficient of cellulose well.

length increasing, we found that phonon

NEMD

calculation,

we

employed293

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

transport inside cellulose nanocrystal is quasiballistic. This feature is attributed to strong
covalent bond of intra-interaction. From theobtained results and linear-relationship between
an inverse of thermal conductivity and an
inverse of long-axis length (Lc), we estimatedthe thermal conductivity of bulk cellulose
nanocrystal as 4.50±0.60 Wm-1K-1, which isFig. 2 Thermal conductivity of CNC (solid line)

similar with reported literatures.

and cellulose single chain (dashed line) as a

As for single cellulose chain, it was expectedfunction of the conduction region length. The

that thermal conductivity of single chain is larger

inset shows extrapolation of inverse of thermal

than that of nanocrystal however these two

conductivity in infinite conduction region length.

thermal conductivities are actually quite similar
with each other, which is different with thereported works on polyethylene [3]. Similar to
the results of poly-dimethylpolysiloxane [4],this result can be understood because non-linear
vibration of side structure arising from thermal
fluctuation impedes phonon transport along
long-axis direction. The non-linear vibration ofside structure is strongly dependent on dihedral

Fig. 3 Thermal conductivity of original single

angle interaction and then we calculated thermal

cellulose chain (solid line), stiff chain (dashedconductivity by strengthening or softeningline) and soft chain (dash-dotted).

spring constant of dihedral angle, which are

References

shown in Fig. 3 denoted as stiff and soft chains.

[1] S. Hashmi, et al., Int. J. Clin. Exp. Pathol., 8

While strength of dihedral interaction increases(2015), 8786.

thermal conductivity, heat conduction along[2] R.W.Hockney, et al., “Computer simulation

long-axis direction is found to be insensitive tousing particles”, Adam Hilger, NY,(1989).

the softening of dihedral interaction.

[3] A. Henry, et al., Phys. Rev. Lett., 10123(2008), 235502.
[4] T. Luo, et al., J. Appl. Phys., 109-7(2011),
074321.

294

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

Heat conduction analysis for phononic crystal
Takuma Shiga
Department of Mechanical Engineering,
The University of Tokyo, Hongo, Bunkyo, Tokyo 113-8656
It has been known that phononic crystal

Figure 1(b) shows thickness-dependent in-plane(PnC) such as superlattice and thin-film with air-thermal conductivity of silicon thin-films atholes leads to significant modulation of phonon

T=300 K. We found that the thickness

dispersion relation, which finally results in large

dependence is largely different with the

reduction

convectional theory considering size effect.

ofthermalconductivity

whilecrystallinity remains. In order to qualitatively
evaluate a thermal conductivity reduction, we
need to directly calculate thermal conductivityof PnC by means of microscopic simulation.
Before considering actual PnCs, we haveconsidered heat conduction of silicon thin films
with a thickness of an order of nanometer byperforming anharmonic lattice dynamics (ALD)
simulation. For interatomic force constants
which are needed in ALD calculation, we
employedtheoptimized

Stillinger-Weberpotential [1]. Surface orientation of silicon film
is set to (001) and 20´20 sampling mech for twodimensional Brillouin zone is chosen for the
calculation of scattering rate [2].

Fig. 1 (a) phonon dispersion relation of silicon

Figure 1(a) shows phonon dispersion relation

thin film with 10 nm thickness. (b) thickness-of silicon thin-film with the thickness of 10 nm.

dependent thermal conductivity at T=300 K.

Because of the presence of surface, there are

References

surface phonon modes labeled by S1-S5. From

[1] Y. Lee, et al., Phys. Rev. B 85, 125204 (2012).

the analysis of eigenvectors, we identified S1 and[2] T. Tadano, et al., J. Phys. Condens. Mater 26,

S2 are Love and Rayleigh waves, respectively.

225402 (2014).

295

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

Comprehensive study on transition between heat conduction
and radiative heat transfer
Takuma Shiga
Department of Mechanical Engineering,
The University of Tokyo, Hongo, Bunkyo, Tokyo 113-8656
It has been so far revealed that thermal

we realized the framework for investigatingtransport across two substrates with nanogap istransition between heat conduction and near-

known to be attributed by phonon-polariton

field radiation.

resonance and tunneling of acoustic phonons
and therefore such a near-field heat transfer canexceed black-body limit. In order to explore
thermally-functionalized devices utilizing nearfield heat transfer such as thermal switch,
thermal rectifier, and so on, microscopic

Fig. 1 Spectral thermal boundary conductancemechanism of transport regime changing fromacross the interface between Ar and heavy-Ar

near-field heat transfer to heat conduction issolids at different temperatures.

needed to be investigated.
In order to realize such investigation, we
chose a combination of spectral heat current andnon-equilibrium molecular dynamics (NEMD)
simulations.

Figuredependent

spectral1

shows

frequency-

thermalboundary

conductance at the interface across solid argon
and heavy-mass argon crystals for different

Fig. 2 Frequency-dependent dielectric function

temperatures.
inelastic

As

terms

of

temperature

increases,of b-cristobalite SiO2 at T=300 K.

interfacial

scattering

References

contribute to overall heat conduction. In addition,[1] K. Sääskilahti, et al., Phys. Rev. B 90,

we134312 (2014).

calculated

dielectric

functionofb-cristobalite SiO2 on the basis of fluctuation-[2] Y. Chalopin, et al., Appl. Phys. Lee. 104,

dissipation theorem. These two results are in

011905 (2014).

good agreement with previous works [1,2], and296

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo                                                                             


3.0
Jr

F

                



2.0

!"

1.0
0.0

XY1

−1.0
−2.0

 
 
       
 

 
  

  
# 
 
!"
         $        ! $ "         
!   "%           % #                       %       $      & '                              # &   &       (    &     )     *     $     +     &  , &               -    $ &$              .  $       # /  $          0 $ +  1   $        $                         #      .  $  1 1  (  1       $            0 , 2    ,             $   * &  $      &   &          $                 !  
& "            
¼

¼



¼

¼



¼

¼

¼

¼

¼

¼

RS
sXY1

  

H

sF

sN
RT

−3.0
−3.0 −2.0 −1.0

¼

N

0.0

1.0

2.0

Δ

3.0

.  3 4        #&          )
            !   & "                  (  1 $         1  1    !    $        "'                         5      $ )  &          &  !(0"   $     ,  !6+4"   $   $               ,  !57"            .   $        )    8   
   $  &  '         $                 )   %     !."$      !."$   !9:"$  !9:"$    !+ "$     !+0"$   !"$ ;< !;"$
    ;< !;"         .$ ;$  ;      &  .  2
  =  1  >$              !?     $   1 1       " '¼297

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo     .         ;            9: 
9:       (  *    $           !"      /  $         !" !"      /   >(                  .     $  & $   !   "   !   "$       &                      $       ,  ! /$ $ $  "       !    "    &                 $  $   $ 
 ! "    &     &              &           8  ,            
57          8     &   !  "   .     !  /"      '     $      8  ,  &    ! #"  &  &   #!  "$      $  
                      8  ,   &          +0          +   4    '   $  @ 
,1 (0   A  &              $     &  * 
 !  / "   !  / "!2"      8  ,   &            $         +0 $ 
8   &       & &   #!//  #  /"  *  //$    
0         +0 
9:       +  9:
    B, ) @  ,   C$        &           (0   &  ;  @ , D $  8  ,  &       
 !  "   !  / " !>"        +0 ; $      ; $    +  .    7 '    $   )            6+4
¼

¼

¼

¼

¼

¼

¼

  E $      8  ,   & $   &  6+4 * $# !  /"  !

" #  !

 /"  !A"#  !  /"   !  /"   !  /" !C"'     8  ,   &         .  9:     .  .      
 !  "   !  /" !D"     $       .  9: $  8   &  
#    & &    !// "  *  /$    . $       .  . $  8  ,  &       & &   # !#  /"  * //$    ( $           ;  ;     8 $      9:  9:     
  &         = 1  >'&        @ F) $ @ ;
$   0)      ) 
 0 ( $ GG 0 $ H  B  $  GI3 ; J 6  !/D" 2/D 0 0 +$  0 7$ 7 +) $ ? 0!7"$ - 03 J H H H   !/D"D
2  J 0 ,3 6 +& B  !KEC" C2D
> @  = $ ? @ ,$  @ ;3 6+& B  !//A" />>>K
A ? @ ,3 J 6 ?  !KKD" (EA
C I ( B, ) 3 0& 6J5 6  !KD"C/% J H @  ,  7 J  3 J 6
L3 0  0   6  !KD2" ED @ ;  ? @ ,3 J 6 ?3 H 
4  !KKE" D2>
E H 6 ;   3 6  ?  !KDC" AC298

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

Spin multipole dynamics and spin transport infrustrated magnets
Hiroaki ONISHI
Advanced Science Research Center, Japan Atomic Energy Agency
Tokai, Ibaraki 319-1195

O−−− (q, ω) =1
1
− Im ⟨0|Oq−−−†
O−−− |0⟩,π
ω + E0 − H + iη qwhere Oq−−− is the Fourier transform of−−
Oi−−− = Si− Si+1
Si+2
. Note that we perform
one DMRG run to obtain the spectral weight
for a set of q and ω, so that we need to do
many runs to scan a wide range of the q-ωspace. The computations are accelerated byparallel simulations using the system B of the
ISSP supercomputer.
We present intensity plots of O−−− (q, ω) in
Fig. 1. We find a finite gap in the quadrupole
regime [Fig. 1(a)], while we observe a gaplessmode at q = π in the octupole regime in accordance with quasi-long-range antiferro-octupole
correlations [Fig. 1(b)]. Since a bound threemagnon cluster is excited with zero energy in22

ω / J2

(a)

1

1

0
2

0
2

(b)

ω / J2

In a spin-1/2 J1 -J2 Heisenberg chain withferromagnetic J1 and antiferromagnetic J2
in a magnetic field, there appear a seriesof spin multipole liquid ground states, i.e.,quadrupole, octupole, hexadecapole, etc. To
clarify magnetic and transport properties in
the quadrupole state, we have studied spinand quadrupole excitation spectra [1] and spin
Drude weight [2] by numerical methods. We
have argued that low-energy excitations aregoverned by bound two-magnon pairs, so thatmagnon pairs would carry spin current.
Here, we extend our focus to the octupole
regime. By using a dynamical DMRG method,
we examine the dynamical octupole structure
factor at zero temperature, given by

1

1

0

0
0

0.5

1

q/π

Figure 1: O−−− (q, ω) at (a) J1 /J2 = −2 andM = 0.35 and (b) J1 /J2 = −3 and M = 0.35,
where M is the magnetization. DMRG resultswith 40 sites.
the octupole state, we envisage that magnon
clusters would contribute to the spin transport.
The analyses of current correlation functions
and wavepacket dynamics are also on going toclarify the relationship between magnetic andtransport properties.

References[1] H. Onishi, J. Phys. Soc. Jpn. 84, 083702(2015); J. Phys.: Conf. Ser. 592, 012109(2015); Physica B 536, 346 (2018).
[2] H. Onishi, J. Magn. Magn. Mater. 479,
88 (2019).

299

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

Molecular Dynamics Simulation of Ferroelectrics
Using a Shell Model V
T. Hashimoto
Research Center for Computational Design of Advanced Functional Materials (CD-FMat),
National Institute of Advanced Industrial Science and Technology (AIST),
Tsukuba Central 2, 1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan

Amorphous BaTiO3 can be formed by sputtering. Its structure is assumed to be a random
network of TiOn local bonding units (LBUs).
By X-ray absorption fine-structure (XAFS)
spectroscopy, n was found to be 6[1]. But, due
to the experimental diﬃculty, the LBU linkage for amorphous BaTiO3 is not clear. We
have developed a classical molecular dynamics(MD) simulation code that can handle shellmodels that reproduce physical properties of
BaTiO3 [2], and applied it to the study of amorphous BaTiO3 .
In the shell model, each atom is composed ofa core and a shell. The intra-atomic core-shellinteraction is expressed by V (r) = c2 r2 /2 +c4 r4 /24, where r is the core-shell distance andc2 and c4 are parameters. The inter-atomic interaction is through the Coulomb interaction
and the Buckingham type shell-shell interaction V (r) = A exp(−r/ρ) − C/r6 , where r is
the inter-atomic shell-shell distance and A, ρ,
and C are parameters. The computations werecarried out in constant temperature and constant pressure (NPT) ensembles. The pressure
and the temperature were controlled by the
Parrinello-Rahman method and the massive
Nose-Hoover chain method, respectively. The
externally applied pressure was set to 0 Pa.
Because we cooled the system slow enough,the densities were very close to each other, andthe susceptibilities were also very close to eachother, both in agreement with experiments[3].
As shown in Fig. 1, the probability of TiO6

1
nTiO

Probability

0.8

1
2
3
4
5
6

0.6
0.4
0.2
0

n[n]Ti[m]Ti

5

6
7
8
5
4.5
4
3.5
3
2.5
2 apex
1.5
1 edge0.5 face044 45 46 47 48 54 55 56 57 58 64 65 66 67 68 74 75 76 77 78 84 85 86 87 88

11
10
9
8
7
6
5
4
3

nBa[m]Tim
4mn

Figure 1: The probability of [m] Ti and thenumber of Ti around O (nTiO ) (upper panel),and n[n] Ti[m] Ti and nBa[m] Ti (lower panel).
was the largest, and the probability of NBOs
(purple) was larger for smaller m. The numberof Ti (Ba) around [m] Ti was larger (smaller)
for larger m. [7] Ti and [8] Ti were in an O rich,Ti rich, and Ba poor environment with fewer
NBOs and many edge or face-sharings, while[5] Ti was in an O poor, Ti poor, and Ba richenvironment with many NBOs.
[1] A. I. Frenkel et al.: Phys. Rev. B 71
(2005) 024116.
[2] S. Tinte et al.: J. Phys.: Condens. Matter
16 (2004) 3495.
[3] T. Hashimoto and H. Moriwake: Physica
B: Condens. Matter 579 (2020) 411799.

300

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

Transport properties of the antiferromagneticclassical XXZ model in two dimensions
Kazushi AOYAMA
Department of Earth and Space Science, Graduate School of Science, Osaka University
Machikaneyama-cho, Toyonaka-shi, Osaka 560-0043

Transport phenomena in magnetic systems
reflect dynamical properties of interacting
spins, such as magnetic excitations and fluctuations. In this work, we theoretically investigate transport properties of two-dimensional
antiferromagnetic insulators, putting emphasis
on how the occurrence of a phase transition isreflected in spin and thermal transports. In
the classical nearest-neighbor (NN) antiferromagnetic XXZ model on the square lattice,)∑(
H = −J
Six Sjx + Siy Sjy + ∆Siz Sjz , (1)⟨i,j⟩the anisotropy ∆ plays a role to control theuniversality class of the transition, i.e., either
a second-order transition at TN into a magnetically ordered state or the Kosterlitz-Thouless(KT) transition at TKT , which respectively occur for the Ising-type (∆ > 1) and XY -type(∆ < 1) anisotropies, while for the isotropic
Heisenberg case of ∆ = 1, a phase transition
does not occur at any finite temperature. The
dynamics of the spins Si is determined by thefollowing semiclassical equation of motion:dSi= Si × Heﬀ
i ,
dt
∑
Heﬀ
= J
(Sjx , Sjy , ∆Sjz ),(2)ij∈N (i)where N (i) denotes all the NN sites of i.
The spin current and the associated spins are respectively given by
conductivity σµν∑Jzs (t) = J
rij (Si × Sj )z ,
⟨i,j⟩
s
σµν

=

1
T L2

∫ ∞

z
z
dt ⟨Js,ν(0) Js,µ(t)⟩, (3)0301

where rij and L denote a vector connecting twosites i and j and the linear system size, respectively. We numerically integrate Eq. (2) withinitial equilibrium spin configurations generated by Monte Carlo simulations and calculate
z (0) J z (t)⟩ at eachthe time correlations ⟨Js,ν
s,µtime step. By using the second order symplectic method, we perform long-time integrations
typically up to t = 100 |J|−1 − 800 |J|−1 withthe time step δt = 0.01 |J|−1 until the timez (0) J z (t)⟩ is completely lost.
correlations ⟨Js,νs,µThe thermal conductivity κµν can also be calculated in a similar way.
It is found that the spin current probes thediﬀerence in the ordering properties, while thethermal current does not. Figure 1 shows
the temperature dependence of the longitudis for the Ising-type
nal spin conductivity σxx(∆ = 1.05), XY -type (∆ = 0.95), andHeisenberg-type (∆ = 1) spin systems. Fors exhibits a diverthe XY -type anisotropy, σxxs ∝gence at TKT of the exponential form, σxx√
exp [B/ T /TKT − 1 ] with B = O(1), while
for the Ising-type anisotropy, the temperatures is almost monotonic withdependence of σxx
out showing a clear anomaly at TN and sucha monotonic behavior is also the case in the
Heisenberg-type spin system [see the dashedcurve in Fig. 1 (c)] [1]. Noting that in the
XY -type spin system, the inter-free-vortex distance ξs diverges toward TKT in the form of√ξs ∼ exp [(π/2)/ T /TKT − 1 ], the significant
s at Tenhancement of σxx
KT can be understoodas a manifestation of the topological nature of

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo(a)ing finite at any finite temperature. It is foundthat the longitudinal spin-current conductivity
exhibits a divergence at the Z2 -vortex transition temperature Tv , while the thermal conductivity only shows a monotonic temperaturedependence with no clear anomaly at Tv [2].
Our results suggest the strong association
between the spin transport and the bindingunbinding topological transition of vortex excitations.

L=24
L=48
L=96
L=192

Δ=1.05108σxsx 6

TN

4
2

(b)

35
30
s 25
σx x 20
15
10
5
0

(c)

Δ=0.95

L=24
L=48
L=96
L=192
L=384

TKT

References
L=24
L=48
L=96
L=192
L=384

Δ=1

10
8

σxsx 6[1] K. Aoyama and H. Kawamura, Phys. Rev.
B 100, 144416 (2019).
[2] K. Aoyama and H. Kawamura, Phys. Rev.
Lett. 124, 047202 (2020).

4
2
0

0.2

0.4

0.6

T/ |J |

0.8

1.0

Figure 1: The temperature dependence of thes inlongitudinal spin-current conductivity σµµ
the (a) Ising-type (∆ = 1.05), (b) XY -type(∆ = 0.95), and (c) Heisenberg-type (∆ = 1)
spin systems. In (a) and (b), red arrows indicate the magnetic and KT transition temperatures, TN /|J| ≃ 0.75 and TKT /|J| ≃ 0.6,
respectively. In (c), a dashed curve representss (T ) curve extrapolated to the thermothe σxx
dynamic limit of L → ∞.
a vortex whose lifetime τv gets longer toward
TKT due to the diﬀusive vortex motion characterized by the relation τv ∝ ξs2 .
Such an anomalous spin transport is also
found in the NN antiferromagnetic Heisenberg
model on the triangular lattice. In contrast to
the XXZ antiferomagnets on the square lattice, the triangular-lattice antiferromagnet is afrustrated system. In the latter frustrated system, a KT-type binding-unbinding topologicaltransition of Z2 vortices is predicted to occur
at Tv , although spins do not order except at
T = 0 with the spin correlation length ξs be-302

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

Study on Dynamic Mechanical Properties of Slide-Ring Gels
using Coarse-Grained MD Simulations
Koichi Mayumi

Material Innovation Research Center (MIRC) and Department of Advanced Materials
Science, Graduate School of Frontier Sciences, The University of Tokyo, 5-1-5
Kashiwanoha, Kashiwa, Chiba 277-8561, Japan
Polyrotaxane (PR) is a supramoleculara coarse-grained model of PR with 7-

polymer composed of many ring molecules, α-membered rings which can slide freely along

cyclodextrins (CDs) and one linear polymer,

the axial chain[2]. By cross-linking rings of thepoly (ethylene glycol) (PEG), threading on

PRs, we prepared SR gels. We conducted uni-these rings. By cross-linking rings belonging

axial elongation simulations of the SR gels,to PRs, Okumura and Ito succeeded inand obtained Young’s moduli ESR, which were

fabricatinggels

lower than theoretical values by the classicalcalled slide-ring gels (SR gels)[1]. SR gels do

affine network model. Then, we analyzed the

not have direct cross-links between polymer

partial chain length distribution in SR gelschains, but topological interlock by figure-of-undereight shaped cross-links composed of two ring

distribution was described by a single

molecules. Because SR gels have slidable

exponential function in un-deformed state,cross-links, they show unique mechanicalwhile it became bimodal by stretching: longer

properties. For example, SR gels show lower

and shorter chains appeared.

supramolecularpolymer

uniaxial

deformation.

While

the

Young’s moduli than conventional fixed crosslinking polymer gels. So far, no molecularmodels have succeeded in representing thequantitative relationship between the sliding
motion of the cross-linking points and the low
Young’s moduli of SR gels. In this study, weinvestigate the relationship between the slidingof figure-of-eight cross-links and Young’s
moduli by means of coarse-grained molecularFig. 1. λ dependence of the average legth ofdynamics simulations.

long and short chain.
The average length of short and long chains

For the modeling of SR gels, we fabricated303

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyochange linearly with λ (Fig.1). We defined thisdecay rate of the segment number of theshorter partial chain length as Nslide. By
considering the chain movement through theslidable cross-links, we have succeeded inestablishing a molecular model for SR gelsbased on the classical three chain model [3](Fig. 2.), and obtained the symple equation :𝐸

𝐸

1Fig. 3. The Young’s moduli of SR gels plottedagainst the cross-linking deensity.

𝑁
𝑁

Where EAffine is the theretical prediction by

Referencesclassical threee-chain model, N0 is the average

1. Okumura Y.; Ito K. Adv. Mater., 2001,

partial chain length in undeformed state.

13, 485
2. Y. Yasuda et al., Macromolecules, 2019, 52,3787.
3. James, H. M.; Guth, E. J. Chem. Phys. 1943,
11, 455.
4. Y. Yasuda et al., ACS Macro Letters, underreview.

Fig. 2. Schematic illustration of the proposed
In order to verify this equation, we compare
the Young’s moduli obtained from thesimulation results and theoretical estimation.
We found that this estimation is realistic
enough to reproduce the simulation results.
Thus we succeeded in.

quantatatively

evaluateing the effect of the sliding on theelasticity of the slide-ring gels.

304

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

Study on Microscopic Dynamics of Polyrotaxanes using FullAtomistic MD Simulations
Koichi Mayumi
Material Innovation Research Center (MIRC) and Department of Advanced Materials
Science, Graduate School of Frontier Sciences, The University of Tokyo, 5-1-5 Kashiwanoha,
Kashiwa, Chiba 277-8561, Japan
Polyrotaxane (PR) is a supramolecularpolymer which consists of an axial linear
polymer and multiple ring-shaped molecules
threaded on the chain. A unique feature of PR
is that the rings can slide and rotate along theaxial polymer, and this feature has beenutilized for fabricating biomaterials[1] andtopological gels[2]. However, the internal
dynamics such as the sliding and rotation havenot been disclosed yet. In this study, weinvestigated the internal molecular dynamics
of PR composed of α-cyclodextrins (CDs) and

Figure 1. Simulation snapshot of the full-poly (ethylene glycol) (PEG) chain in solution

atomistic model for (a) PRs, (b) free CD, (c)byfree PEG, and (d) MSD for translational

means

of

full-atomistic

molecular

dynamics (MD) simulations [3].
We

measured

the

mean

diffusion of CD in PR(red solid line), freesquaredCD(red dotted line), and free PEG (blue dotted

displacement and estimated the diffusionline)

coefficient of free PEG, free CD, and CD in
PR. The simulation results agree well with theTable 1. Diffusion coefficients of CD in PR,

diffusion coefficients measured by the quasi-

free CD, and free PEG from the asi-elastic

elastic neutron scattering experiment (Table 1).

neutron scattering (QENS) experiment andFrom this result, we consider that our model is

MD simulation.

realistic enough. Then we tried the analysis of
intramolecular motion: the sliding dynamics.

305

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

In order to evaluate the sliding dynamics in PR,

In conclusion, we succeeded in representing

we tracked the one-dimensional trajectory ofthe diffusion coefficient for the sliding motion

CDof PR by the combination of the jump diffusion

along

PEG.

The

mean

squareddisplacement corresponding to the sliding

term with the Einstein-Stokes equation:

motion (MSDslide) increased linearly with time,𝐷

which suggests that the sliding motion is

∝

𝑘 𝑇
exp
𝜂

𝐸
𝑘 𝑇

diffusive (Fig. 2.). From the time dependenceThis relationship is valuable for controlling theof MSDslide, we calculated from the diffusion

sliding speed in polyrotaxanes.

coefficient for the sliding motion Dslide at 300
K as 1.42Å/ns, which is quite smaller value
than that of the translational diffusion of CD in
PR.

Fig. 3 Temperature dependence of the
diffusion coefficient of free CD, free PEG
and sliding motion.
Fig. 2. MSDslide of CD along the PEG chain.

References
In order to clarify the molecular origin of this(1). Yui, N. et al., Chem. Eur. J. 2006, 12,

slow sliding dynamics, we analyzed the6730

temperature dependence of Dslide (Fig. 3). We(2). Okumura, Y.et al., Adv. Mater. 2001,

found that the diffusion coefficients of the13, 485

sliding dynamics was dominated by the energy(3). Yasuda, Y. et al., J. Am. Chem. Soc.

barrier between CD and PEG. This energy

2019, 141 (24), 9655.

barrier was 8.92 kJ/mol, which corresponds to(4). Liu, P. et al., J. Phys. Chem. C 2012,

the energy potential barrier for CD on PEG,

116, 17913.

2kcal/mol = 8.37 kJ/mol [4].

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

Comparison of X-ray pinhole topographs experimentallyobtained and computer-simulated based on the Ewald-Laue andTakagi dynamical theories
Kouhei OKITSU
Institute of Engineering Innovation, School of Engineering,
The University of Tokyo, 2-11-16 Yayoi Bunkyo-ku Tokyo 113-8656

Fig.1 [Sx(T-T)], [Ex] and [Sx(E-L)] (x∈{h, v}) show T-T simulated, experimentally obtainedand E-L&FFT simulated eight-beam pinhole topographs for horizontally (x = h) andvertically (x = v) polarized incident X-rays [2].

The present author verified the simple relation(x = v) polarized incident X-rays whose photon

between the n-beam Ewald-Laue (E-L) and Takagi

energy was adjusted or assumed to be 18.245 keV.

(T-T) X-ray dynamical theories explicitly [1-4] forThe

the first time. It can be described with Fourier

polarizationtransformation. Then, the behavior of diffracted

synchrotron radiation by using the four-quadrant

X-rays in a perfect crystal can be described by bothphase retarder system [1, 2, 4] at the BL09XU ofthe E-L and T-T dynamical theories.

SPring-8. A [1 -1 1]-oriented floating zone silicon

experiment
state

was
of

performed
the

with

the

monochromatedFigs. 1 [Sx(T-T)], [Ex] and [Sx(E-L)] (x∈{h, v})

crystal with a thickness of 9.6 mm was used as theshow T-T simulated, experimentally obtained andsample. The indices of reflections are as described

E-L&FFT

pinholein Fig. 1 [Sh(T-T)]. The T-T simulated topographs

topographs for horizontally (x = h) and verticallywere computed by integrating the T-T equation

simulated

eight-beam

307

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

Fig. 2 Enlargement of the 0 6 6 transmitted-reflected images in Fig. 1 [2].
with the finite-difference method. However, theconcerning the too large R factors in protein

E-L&FFT topographs were obtained by fast

crystallography by estimating them with n-beam

Fourier transforming the solution of the E-L

approximation in place of the conventional two-

dynamical theory.

beam approximation.

References

Figs. 2 [Sx(T-T)], [Ex] and [Sx(E-L)] areenlargements of the 0 6 6 transmitted-reflected[1] K. Okitsu, Y. Imai, and Y. Yoda: Recent

images in Figs. 1 [Sx(T-T)], [Ex] and [Sx(E-L)].

Advances in Crystallography, (2012) pp. 67-86.

Good agreements among Figs. 2 [Sx(T-T)], [Ex]InTech Open. http://dx.doi.org/10.5772/47846.

and [Sx(E-L)] are found both for x=h and for x=v,[2] K. Okitsu, Y. Imai, Y. Yoda and Y. Ueji: Acta

respectively. However, evident discrepancies are

Cryst. A75 (2019) 474-482;

also found between the cases of x=h and x=v.

https://doi.org/10.1107/S2053273319001499.

The present author has a hypothesis that too[3] K. Okitsu, Y. Imai, and Y. Yoda: Acta Cryst.

large R factors for protein crystals is caused by

A75 (2019) 483-488;

interferences of multiple-reflected wave fields andhttps://doi.org/10.1107/S2053273319002936.

bankruptcy of the two-beam approximation that[4] K. Okitsu: J. Jpn. Soc. Synchrotron Rad. Res.

has been used for over 100 years. The present

33 (2020) 61-80 [Invited].

author has prepared to verify this hypothesis308

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

Numerical verification of the higher-order ETH inquantum XXZ spin ladder
Eiki IYODA
Department of Physics, Tokai University
Kitakaname, Hiratsuka-shi, Kanagawa 259-1292

Recently, isolated quantum many-body systems, which are realized in artificial systems such as cold-atoms or superconducting
qubits, have attracted much attention, because
they are useful to study fundamental problems in statistical mechanics. Especially, it
has been shown that even an isolated quantum many-body system exhibits thermalization. The eigenstate thermalization hypothesis(ETH) [1] is a promising mechanism of thermalization, which is closely related to quantum
chaos. Many numerical studies have shownthat ETH holds in non-integrable systems.
We propose a higher-order generalization of
the ETH, to which we refer as the k-ETH (k =1, 2, · · ·) [3]. The lowest order ETH (k = 1) isthe conventional ETH. The k-ETH gives many
conditions which are obtained by comparing
Hamiltonian dynamics with the Haar-random
unitary of the k-fold channel. As a non-trivial
contribution of the higher-order ETH, we showthat the Page correction of the kth Renyi entanglement entropy of individual energy eigenstates originates from the k-ETH.
In this study, we numerically investigated 2ETH for the one-dimensional spin-1/2 Heisenberg XXZ ladder model, which is composedby two chains whose lengths are p and p + 1
respectively. In order to numerically investigate 2-ETH, we should discuss finite size scaling of indicators of 2-ETH by changing thesystem size. The 2-ETH gives many conditions. For example, 2-ETH states that diagonal and oﬀ-diagonal matrix elements such as

Figure 1: (a) d-dependence of I2 (SX ). (b) λdependence of a.
⟨Ei Ej |O|Ei Ej ⟩ equal some value. There aretotally d2 equalities, where d is the dimensionof the energy shell. We adopt I2 (O) as an indicator of 2-ETH, which is defined as the maximum value of the deviation of the d2 equalities. I2 (O) is a straightforward extension of
the indicator of strong 1-ETH [2]. By numerically exact diagonalization, we calculate I2 (O)by changing parameters about energy shell andthe Hamiltonian.
We numerically calculated I2 (SX ), where SX
is the partial swap operator which swaps theregion X between the original system and a
replicated system. Figure 1(a) shows the ddependence of I2 (SX ) with λ = 1, where
λ is the coupling constant in the rungs andtunes the integrability of the model. In nonintegrable systems, the indicator I2 (SX ) decays
polynomially with d, which implies the 2-ETH
for SX holds in non-integrable systems.
We fit the numerical data of log I2 against
a fitting function f (d) := −a log d + b. Figure1(b) shows the λ-dependence of the exponent

309

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

a. In the non-integrable cases (λ ̸= 0), theexponent a takes a suﬃciently large value. Especially, when the size of the region X equals
p, the exponent a ∼ 1, which is consistent withthe typicality prediction.
In this report, we have numerically shownthat the higher-order ETH holds by using theindicator I2 (O). Numerical studies about finite size scaling with other indicators and theweaker version of higher-order ETH are left as
future issues.

References[1] M. Rigol, V. Dunjko, and M. Olshanii,
Nature 452, 854 (2008).
[2] H. Kim, T. N. Ikeda, and D. A. Huse,
Phys. Rev. E 90, 052105 (2014).
[3] K. Kazuya, E. Iyoda,
arXiv:1911.10755.

T. Sagawa,

310

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

Analysis on Structuring and Dynamics of Ionic Liquid
Forming Electric Double Layer as a Local Interfacial Field
Ken-ichi FUKUIDepartment of Materials Engineering Science, Graduate School of Engineering Science,
Osaka University, Machikaneyama, Toyonaka, Osaka 560-8531
The interfacial ionic-liquid (IL) behavior
causing bias stress in the electric double layer
organic field effect transistor (EDL-OFET) (Fig.
1) has been investigated in the molecular scaleusing MD simulations combined withexperimental techniques [1].
The MD simulations were performed withGROMACS-5.0.7. The rubrene substrate was oflinear dimensions 7.2165 × 7.1930 × 2.7760 nm.
EMIM-FSI ion pairs were sandwiched between
the rubrene substrate (2.7760 nm) and vacuum
layer (~9.5 nm). Simulations were performed forthe uncharged rubrene surface (σ = 0) and for the
positively charged surfaces (σ = 0.6, 1.2, and 1.8
μC cm-2). Single point charges were placed oneach carbon atom consisting of a tetracene
backbone of a rubrene molecule corresponding to
the HOMO of rubrene, and its magnitude wasdetermined by the overall surface charge density
σ. The systems were equilibrated at a constantvolume for 100 ps, followed by a constantvolume simulation for 20 ns.
Fig. 2 shows the plane-averaged charge density
profiles with respect to the amount of injectedcharge into rubrene electrode and number density
profiles of atoms in EMIM cation and FSI anion,
respectively, at the charge density of 1.8 μC cm-2.
At this potential, EMIM cations and FSI anions
form a checkerboard structure with their major
molecular axis perpendicular to the surface. Thus,
negative atoms of -SO2F of FSI anions arepointing the rubrene surface (O×2 and F). Thisconfiguration probably causes the bias stress by
efficiently trap the hole carrier at the interface
with structured interfacial IL as observed by force
curve measurements with EC-FM-AFM [1].

Fig. 1: (a) Molecular structures of rubrene andEMIM-FSI. (b) Schematic of EDL-OFET using
rubrene semiconductor and IL gate dielectrics.

Fig. 2: (a) Plane-averaged charge density profiles
with respect to the amount of injected charge into
rubrene electrode. (b),(c) Number density profilesof atoms in the EMIM cation and the FSI anion,
respectively, at the charge density of 1.8 μC cm-2.
Each line corresponds to the shaded atom withthe same color. Gray-shaded area is the rubreneelectrode region. Dotted line at z = 0.04 nm
suggests the position of negative charge density
next to the rubrene electrode surface.

Reference[1] D. Okaue, I. Tanabe, S. Ono, K. Sakamoto, T.
Sato, A. Imanishi, Y. Morikawa, J. Takeya, K.
Fukui, J. Phys. Chem. C 124 (2020) 2543.

311

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

Calculation of ordered structures and their optical properties
of soft materials
Jun-ichi FUKUDA
Department of Physics, Kyushu University, Motooka 744, Nishi-ku, Fukuoka 819-0395
This year we have studied how a topologicalthe effect of the variation of the numericaldefect in a thin cell of a nematic liquid crystal

aperture. We have found that the core of alook under optical microscope. Such defects are

topological defect appears as a dark spot, andcommonly observed in experimental labs

also that the spot becomes sharper as thestudying

numerical aperture increases. More details will

liquidcrystals.

However,usual

be presented in a manuscript in preparation.

identifications of topological defects in anematic liquid crystal relies on the observation

Another subject is the lattice orientation of aof the orientation profile around the topological

cubic lattice of cholesteric blue phase liquid

defects. How the image should be at the centercrystalof topological defects has given little attention

orienting surfaces. Previous experiments foundso far. It is because the length characterizingthatthe defect core is ~10nm, much smaller than theorientation with respect to the direction of theresolution of optical imaging.

surface orientation. We have simulated the

the

in

contact
cubic

with

latticeunidirectionallyexhibits

specific

We have made use of the techniquestructures of cholesteric blue phases in such adeveloped previously [1] which reproduced

setup using the Landau-de Gennes theory

excellently the microscope images of a thin cell

describing the orientational order of the liquidof a chiral liquid crystal. We solve the fullcrystal by a second-rank tensor. This study is

Maxwell equations for the electromagnetic

now under way.

wave to investigate the response of the liquid
crystal to monochromatic incident light with

References

given wavevector. By doing calculations for[1] A. Nych, J. Fukuda, U. Ognysta, S. Žumer

different

and I Muševič, Nature Phys. 13, 1215 (2017).

incident

wavevectors,

one

can

construct a microscope image, whose details[2] M. Takahashi et al. J. Phys. D: Appl. Phys.

are presented in Ref. [1].

51, 104003 (2018)

Calculations have been done to investigate

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

Ultrafast relaxation of inhomogeneous electron
distribution around the Dirac cone in graphene
Shota ONO
Department of Electrical, Electronic and Computer Engineering, Gifu University
Gifu 501-1193

(β)

!% #$

!"#

!"/
!"#

!"!
$!"#

!"#

%&*"+'()%&!"!'()%&,-'()

!"!
$!"#
$!"#

(ν)

∂nq /∂t = (∂n/∂t)ph-e , where fks is the electron distribution function with the wavevector k = (kx , ky ), the valleys ν =K and K’,(β)
and the band index s = ±1. nq is thephonon distribution function, where β and q
denote the phonon mode and wavevector, respectively. The collision terms include theelectron-electron (e-e) and the electron-phonon(e-ph) scattering contributions. The e-e interaction matrix element is calculated by using
the expression given in Ref. [3]. The expression for the e-ph interaction matrix elements
given in Ref. [4] were used to compute thescattering rates, while diﬀerent e-ph couplingstrengths [5] were used for the coupling withthe Γ and K phonons, respectively. The initialcondition for the electron distribution is imposed to show the inhomogeneous population
around the Dirac cones [6].
Figure 1 (upper) shows the time-evolution(ν)
of fk,+1 : The excited electrons inhomogeneously populated around the Dirac cone initially go to the K point due to the e-e scattering mainly, followed by electron thermalization due to both the e-e and e-ph scattering.
The electron relaxation is slowed down from∼ 50 fs to more longer when the e-e scattering is neglected (lower). Using the electron
distribution, the time-evolution of absorption!".

%&,-'()%&*"+'()%&!"!'()

!"!

!% #$

It has been known that the inhomogeneouselectron distribution can be created around the
Dirac cone of monolayer graphene after photon absorption [1]. Motivated by recent experiment using ultrafast laser pulse [2], we simulate the ultrafast dynamics of the photoexcited Dirac electrons coupled with the optical phonons by solving the Boltzmann equa(ν)
tion: ∂fks /∂t = (∂f /∂t)e-e + (∂f /∂t)e-ph and

!"!
!" #$

!"#

$!"#

!"!
!" #$

!"#

$!"#

!"!
!" #$

!"#

Figure 1: The time-evolution of the electron
distribution around the Dirac cone: In upper
figure, both the e-e and e-ph scattering are considered, while in lower figure the e-ph scattering is considered only.
spectra can also be calculated. We expect thatthe present study will stimulate further experimental work on this system.
The author thank I. Katayama (Yokohama
National University) for useful discussions.

References[1] A. Grüneis et al.: Phys. Rev. B 67 (2003)
165402.
[2] S. Nara et al.: JPS 2018 Autumn, 11aPS20.
[3] T. Ando: J. Phys. Soc. Jpn. 75 (2006)
074716.
[4] S. Piscanec et al.: Phys. Rev. Lett. 93(2004) 185503.
[5] M. Lazzeri et al.: Phys. Rev. B 78 (2008)
081406(R).
[6] M. Trushin et al.: Phys. Rev. B 92 (2015)
165429.

313

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

Boltzmann equation solver for nonequilibrium
electrons and phonons in solids
Shota ONO
Department of Electrical, Electronic and Computer Engineering, Gifu University
Gifu 501-1193

By taking the wavevector average of the
distribution functions, we calculate the timeevolution of the electron distribution f (ε) andthe phonon distribution n(ω), where ε and ω
are the electron energy and the phonon frequency, respectively. The time-evolution off (ε) and n(ω) is regulated by the Boltzmann
equation ∂f (ε)/∂t = (∂f /∂t)e-e +(∂f /∂t)e-ph +(∂f /∂t)laser and ∂n(ω)/∂t = (∂n/∂t)ph-e ,where the laser excitation as well as the e-e
and e-ph collision terms are considered.
We apply our BENEP code to the calculation of the femtosecond infrared photolumi-

0.0 ps0.2 ps
0.5 ps1.0 ps3

15x10

PL (arb.)

It is of importance to establish a theoretical scheme for studying the ultrafast electron dynamics of laser-excited solids, where
the nonequilibrium electron distribution approaches thermal distribution within a picosecond time scale through the electron-electron(e-e) and electron-phonon (e-ph) scattering.
Recently the author has developed the Boltzmann equation solver to calculate the timeevolution of the electron and phonon distribution [1], where the free-electron, the Debye
phonon, and the deformation potential models are used to describe the electron densityof-states (DOS), the phonon DOS, and theEliashberg function (the e-ph coupling function), respectively. In the present study, basedon ab initio approach to calculate these threefunctions, we develop the Boltzmann Equation solver for Nonequilibrium Electrons andPhonons, BENEP, which enables us to study
the ultrafast electron and phonon dynamics
of realistic materials absorbing a femtosecond
laser pulse. It would be useful to analyze
the pump-probe experimental data beyond thewell-known two-temperature model [2].

10

5

0
0.2

0.4

0.6
0.8
Photon energy (eV)

1.0

1.2

Figure 1: The time-evolution of PL in silver:
The excitation density, the eﬀective Coulombinteraction strength, and the electron diﬀusion
time are assumed to be 460 J/cm3 , 0.37 eV,
and 0.5 ps, respectively.
nescence (PL) of metals. Figure 1 shows thetime-evolution of PL spectra in silver: The PL
spectra changes with time because the electron
distribution thermalizes through the e-e scattering and because the electron energy is transferred to phonons. We find that the agreement
between the numerical simulation and the experimental data is good, which will be demonstrated in our future paper [3].
The author thank T. Suemoto (Toyota Physical and Chemical Research Institute) for useful discussions.

References[1] S. Ono: Phys. Rev. B 97 (2018) 054310.
[2] P. B. Allen: Phys. Rev. Lett. 59 (1987)
1460.
[3] S. Ono and T. Suemoto, arXiv:2004.09000.

314

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

Entanglement in Remote Electron-phonon Systems
Created by Photoirradiation
Kunio ISHIDA
School of Engineering and Center for Optics Research and Education
Utsunomiya University, Yoto, Utsunomiya, Tochigi 321-8585
Recent experiments on quantum entanglement

We derived relevant composite modes obtained

between remote systems has revealed that theby the Heisenberg equation of motion for 𝑐𝑖 , andirradiated photons mediates quantum correlation

found that the heterodyne detection of it will

between noninteracting material systems[1,2].

help us distill the information on phonon

Since,

and/orcorrelation which slowly increases after a certain

dynamical behavior has not been understoodnumber of phonons is created. As shown in thetheoretically, we study the dynamics ofprevious

entanglement

remotebetween remote diamond crystals[1], thesystems by irradiation of a quantized light pulse.

entanglement is generated as a result of

We employed a model of coupled electron-

measurement on the scattered light, i.e.,phonon-photon systems described by[3]measurement corresponds to a disentanglement

however,

its

generation

mechanism

between𝑗

𝜎 +1
2Solving

entanglement

quantum correlation before measurement is also𝑗time-dependent

phononthe projection hypothesis, this means that the+ {∑3𝑖=1 𝜈𝑖 (𝑐𝑖 † + 𝑐𝑖 ) + 𝜆}𝜎𝑥 ].

the

on

process between photons and phonons. Underℋ = ∑3𝑖=1 𝛺𝑖 𝑐𝑖 † 𝑐𝑖 + ∑2𝑗=1 [𝜔𝑎𝑗† 𝑎𝑗 + {𝜇(𝑎𝑗† +𝑎𝑗 ) + 𝜀} 𝑧

study

important to find appropriate methods for

Schrödinger

entanglement control. These results will help us

equation by numerical calculation on the Systemdesign a method of coherent control of

B at ISSP, we found that the quantum mutual

entanglement generation.

information for phonons reveals the dynamics ofphonon entanglement generation.

References
entanglement[1] K. C. Lee, et al., Science 334, 1253(2011).

generation dynamics are also studied and we[2] F. Borjans, et al., Nature 577, 195 (2020).

showed that the scattered light, e.g., Stokes light,[3] K. Ishida, Eur. Phys. J. D73, 117 (2019).

Detection

methods

ofcarries information on the phonon entanglement.

315

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

Molecular simulation of patchy particles
Takamichi TERAO
Department of Electrical, Electronic, and Computer Engineering,
Gifu University, Yanagido 1-1, Gifu, 501-1193
In colloidal systems, many studies have been

structures of colloidal crystals have been

performed to clarify the properties of sphericalanalyzed for various crystal structures. It has

particles with isotropic interactions. Recently,

been revealed that a wide photonic bandgap

there has been much interest in the studies ofexists in a colloidal crystal with diamondpatchy particles with anisotropic interactions. It

structure. Melting transition of such crystals

is of interest to create novel functional

have been studied by Monte Carlo simulations,materials using these particles. In this study, theand it was confirmed that the crystals withstructural formation of patchy particles was

diamond structure becomes thermodynamicallystudied by molecular simulations using cluster

stable if the attractive interaction energyMonte Carlo method. To describe the effective

between patchy particles is about 10kBT.

interaction between patchy particles, the KernFrenkel model has been employed [1].
First, computer simulations are performedfor 2-patch particles. Each particle has degrees
of freedom of rotation in three-dimensional
space, and its motion was constrained in aplane. In order to clarify the phase diagram, weanalyzed the structure of the system using thegeneralized local bond order parameter that we
Fig. 1: Snapshot of tetrahedral patchy particles.

have developed in previous studies [2]. As aresult, a fluid phase, a hexagonal phase, aKagome phase, and their coexistence phase

References

have been observed.

[1] F. Romano, E. Sanz, P. Tartaglia, and F.

Next,

we

studiedtetrahedral

Sciortino: J. Phys.: Condens. Matter 24 (2012)

patchyparticles in a three-dimensional system (Figure

064113.

1), and the stability with diamond structure was[2] R. Mizuno, K. Okumura, J. Oguri, and T.

numerically investigated. The photonic band

Terao: Mol. Simul. 45 (2019) 743.

316

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

Eﬀects of Dzyaloshinskii-Moriya interactions andexchange randomness on low-temperature specificheats of spherical-kagome spin-systems
Yoshiyuki Fukumoto
Department of Physics, Faculty of Science and Technology, Tokyo University of Science
2641 Yamazaki, Noda, Chiba 278-8510

For cluster magnets described by Heisenberg models, in general, stepwise structures
are expected in zero-temperature magnetization curves. However, for a spherical-kagome
system {W72 V30 } [1], which contains 30 vanadium ions with S = 1/2, it was reported
that there is no such structure up to 50 T
in a magnetization measurement at 0.5 K [2].
So far, the eﬀects of (i) a distribution in thestrength of the exchange interaction [2] and(ii) the Dzyaloshinskii-Moriya (DM) interaction [3] have been investigated. If in (i) the distribution width is about 30% of the average exchange interaction, or if in (ii) the DM interaction strength is about 10% of the exchange interaction, it has been shown that the low temperature magnetization process of {W72 V30 }can be explained.
In this study, we focus on the eﬀects of (i)and (ii) on the specific heat. To be specific, wefirst notice the fact that the recently reported
experimental results for the low-temperature
specific heat of {W72 V30 } in several magnetic
fields do not match the calculated results in the
Heisenberg model; (a) the peak at about 2K inthe calculated result is not seen in the experimental result, and (b) the observed magneticfield dependence is very smaller than the calculation [4]. In order to resolve these two discrepancies, (a) and (b), we use the method ofthermal quantum pure state [5] to calculate thespecific heat, taking (i) and (ii) into account.
As a result, regarding (a), we conclude that317

the distribution of low-energy singlet states isaﬀected by the distribution of exchange interactions to wipe out the low-temperature peak.
The important finding is that the lowest singlet excitation-energy in the Heisenberg model
is a fraction of the lowest triplet excitationenergy, which leads to the distribution widthof exchange interaction being about 10% of theaverage exchange interaction (a fraction of theprevious study [2]). Regarding (b), we foundthat the DM interaction can have the eﬀect ofsuppressing the magnetic-field dependence ofthe energy eigenvalues, whose mechanism canalso be illustrated by the two-spin model.
In conclusion, the specific heat of the 30-site
system has been quantitatively evaluated up
to suﬃciently low temperatures, and we havequalitatively discussed the diﬀerences and similarities between the eﬀects of two types of perturbations (i) and (ii) on the specific heat inmagnetic fields and the magnetization process.

References[1] N. Kunisada et al., PTEP 2014, 41I01(2014).
[2] J.
Schnacketmat/1304.2603v1.

al.,

arXiv:cond-[3] Y. Fukumoto et al., J. Phys. Soc. Jpn. 87,
124710, 2018.
[4] T. Kihara et al., PRB 99, 064430 (2019).
[5] S. Sugiura and A.Shimizu, Phys. Rev.
Lett. 108, 240401(2012).

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

Data-driven determination of a spin Hamiltonian of
KCu4P3O12
Ryo TAMURAInternational Center for Materials Nanoarchitectonics,
National Institute for Materials Science,
1-1 Namiki, Tsukuba, Ibaraki, 305-0044

The origin of physical properties in materials is often understood through an eﬀectivemodel in materials science. Many methods
have been developed in order to construct theeﬀective model for a target material and they
are mainly divided into two groups. One is abinitio calculations which determine the model
parameters in an assumed eﬀective model only
by giving basic information of the target material. The other is a data-driven approachin which model parameters are determined so
as to fit the experimentally measured data inthe target material. In the latter case, from aview point of eﬀective model estimations, datadriven techniques are considered to be eﬃcientfor the acceleration of automatic search for appropriate model parameters[1] and the extraction of relevant model parameters[2, 3, 4].
Here, by a data-driven approach, aspin Hamiltonian as an eﬀective model ofKCu4 P3 O12 is determined. The crystal structure of KCu4 P3 O12 is defined by a =7.433
Å, b = 7.839 Å, c = 9.464 Å, α = 108.28◦ ,β = 112.68◦ , γ = 92.73◦ , and space group:
P1̄[5]. Cu(II) ions have S = 1/2 isotropic
Heisenberg spin, but magnetic properties ofthis compound have not been reported yet.
The lattice structure of Cu ions can be regarded as a zigzag chain consisting of eight Cuions. Thus, a quantum Heisenberg model ona zigzag chain is targeted for a spin Hamiltonian of KCu4 P3 O12 to be estimated. We estimate superexchange interactions between Cu

ions in this target model by a data-driven approach with susceptibility and magnetization
curves experimentally measured. The estimated model parameters are J1 = −8.54 meV,J2 = −2.67 meV, J3 = −3.90 meV, andJ4 = 6.24 meV, which describe the measured
results well.
Once an estimated spin Hamiltonian is established, we predict, through theoretical analysis of the Hamiltonian, various magnetic
properties which cannot be or have not beenmeasured, such as magnetic specific heat, magnetic entropy, spin configuration, and spin gap.
These predictions are helpful to propose a further experimental plan and design.
This work is the collaboration work withKoji Hukushima, Akira Matsuo, Koichi Kindo,and Masashi Hase.

References[1] R. Tamura and K. Hukushima, PLoS
ONE 13, e0193785 (2018).
[2] H. Takenaka, K. Nagata, T. Mizokawa,and M. Okada, J. Phys. Soc. Jpn. 83,
124706 (2014).
[3] R. Tamura and K. Hukushima, Phys. Rev.
B 95, 064407 (2017).
[4] H. Fujita, Y. O. Nakagawa, S Sugiura, and M. Oshikawa, Phys. Rev. B 97,
075114 (2018).

318

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[5] H. Eﬀenberger, Zeitschrift für Kristallographie 180, 43 (1987).

319

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

Establishment of new analysis method
for extend X-ray absorption fine structurewith sparse modeling
Hiroyuki KUMAZOE
Institute of Industrial Nanomaterials, Kumamoto University
Kurokami, Chuo-ku, Kumamoto, 860-8555

Measurement of extended X-ray absorption
fine structure (EXAFS) is one of the primarymethods to obtain local structure information
around a specific element with atomic scale.
We are building a new method for the analysis
of EXAFS by the sparse modeling. This year,
we focused on Bayes free energy as an information criteria for choosing the regularization
parameter.
To perform the sparse modeling on the EXAFS signals, we employed a simplified modelof the EXAFS formalism based on the singlescattering approximation [1]:y(k) =χ(k) · k 3
=

N
X
k2
j=1

Rj2

exp



−2k 2 σj2Rjexp −2Λ× (aj sin 2kRj + bj cos 2kRj ),(1)matrix, X incorporates other terms in Eq. (1).
In this method, some elements of the coefficient, w is are suppressed to exactly zero withthe moderate value of λ.
Let us consider the sparse regression
problem for the LASSO estimations within
Bayesian inference to choose the regression parameter, λ. In Bayesian inference, this linear
regression problem is optimized by maximizing the posterior probability. Here, we introduce the Bayesian free energy [2, 3] as an information criterion to optimize a regularization
parameter and to extract the physical model
appropriately. As a result, we estimated theradial distribution function with fewer bases
than in previous studies [1].
This work was supported by JST CREST
Grant Number JPMJCR 1861, Japan.

where, y(k) is the k 3 weighted EXAFS signal, χ(k) as a function of the photoelectronwavenumber, k. Rj is the distance from thefocused atom, j. σj is the Debye–Waller factor for atom, j. Λ is the mean free path. To
obtain the local structure information, we convert Eq. (1) to a linear regression problem andsolve it with L1 regularization.
12ŵ = arg min ky − Xwk2 + λkwk1 , (2)2wReferenceswhere the regression coefficient is w =(. . . , aj , bj , . . . )T .
The response vector, ycomes from the target data and the predictor[3] Y. Mototake, Y. Igarashi, H. Takenaka,
K. Nagata, and M. Okada, J. Phys. Soc.
Jpn. 87, 114004 (2018).

320[1] I. Akai, K. Iwamitsu, Y. Igarashi,
M. Okada, H. Setoyama, T. Okajima, and
Y. Hirai, J. Phys. Soc. Jpn. 87, 074003(2018).
[2] Y. Igarashi, H. Takenaka, Y. NakanishiOhno, M. Uemura, S. Ikeda, andM. Okada, J. Phys. Soc. Jpn. 87, 044802(2018).

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

Competiton between dipolar and short-range
interactions on the magnetic thin films
Hisato KOMATSU, Yoshihiko NONOMURA, and Masamichi NISHINONational Institute for Materials Science, Namiki, Tsukujba, Ibaraki 305-0044

Magnetic thin films are known to exhibit avariety of orderings when the exchange andthe dipolar interactions compete with eachother [1]. These systems are described by thefollowing Hamiltonian:∑H = −J⟨i,j ⟩+D′−ηSi · Sj∑  Si · Sj
i<j

∑(

i

(a) η=0,H=0


z

Si

r3
)2

ij

−H

−

3(S i · r ij )(S j · r ij


)

5
rij



∑

z

Si .

i

321

1.0

0.5

0.5

0.0

0.0

-0.5

-0.5

-1.0

-1.0

(c) η=0,H=2(1)

Here, r ij = t (xij , yij , 0) ≡ j − i means thediﬀerence between the position i and j. Previous studies revealed that this system has the
spin-reorientation(SR) transition between the
stripe and the planar ferromagnetic phase [2, 3]
driven by the anisotropy η. However, thesestudies are limited to the cases when the magnetic field H does not exist.
We investigate the relation between the
structure of the ordered phase and the parameters η and H using the Monte-Carlo(MC)simulation. In the simulation, naive MC algorithm require O(N 2 ) simulation time, andthus we adopt the O(N ) SCO algorithm [4].
The coupling constant of interactions are fixed
as J = D′ = 1. Figures 1 and 2 show the snapshots and the η − T phase diagrams given by
our investigation, respectively. These figures
show that the SR transition exists no matter
whether H is zero or nonzero.
The numerical calculations of this study
were mainly performed on the Numerical Materials Simulator at National Institute for Ma-(b) η=5,H=0
1.0(d) η=5,H=2
1.0

1.0

0.5

0.5

0.0

0.0

-0.5

-0.5

-1.0

-1.0

Figure 1: Snapshots of the system at several
cases. The vectors and color mean the values of(S x , S y ) and S z , respectively. The parametersη and H are fixed at (a)η = 0, H = 0, (b)η =5, H = 0, (c)η = 0, H = 2, and (d)η = 5, H =2, and the temperature T = 0.1. Note that thewidths of the stripes depend on H.
terials Science. We applied for the uage of
the ISSP Supercomputer in preparation for the
maintenance period of this simulator. However, we have ended this period without needfor the new large-scale simulation, and do nothave used the ISSP Supercmputer.

References[1] K. De’Bell, A. B. MacIsaac, and J. P.
Whitehead, Rev. Mod. Phys. 72, 225(2000).

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[2] J. P. Whitehead, A. B. MacIsaac, and K.
De’Bell, Phys. Rev. B 77, 174415 (2008)[3] M. Carubelli, O. V. Billoni, S. A. Pighin,
S. A. Cannas, D. A. Stariolo, andF. A. Tamarit, Phys. Rev. B 77, 134417(2008).
[4] M. Sasaki and F. Matsubara, J. Phys. Soc.
Jpn. 77, 024004 (2008).

(a)10
8

stripe

η

6

para
4
2 planar ferro
0
0.0

0.5

(b)10
8

1.0

1.5

T
stripe

η

6

para
4
2 planar ferro
0
0.0

0.5

1.0

1.5

T

Figure 2: The η −T phase diagrams at (a)H =0 and (b)H = 2. The closed and open circles mean the first- and the second-order phasetransition.

322

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

Competiton between dipolar andDzyaloshinskii-Moriya interactions on the magnetic
thin films
Hisato KOMATSU, Yoshihiko NONOMURA, and Masamichi NISHINONational Institute for Materials Science, Namiki, Tsukujba, Ibaraki 305-0044}

on the value of D′ and η. It makes the triangular lattice under certain conditions, and thesquare lattice under the others. Furthurmore,
the chiral magnetic phase also shows thechange of the direction of the structure depending on the parameters. We calculate theorientational order parameter of these phases,
and investigated the condition under which
each structure appear.
The numerical calculations of this study
were mainly performed on the Numerical Materials Simulator at National Institute for Materials Science. We applied for the uage of
the ISSP Supercomputer in preparation for the
maintenance period of this simulator. However, we have ended this period without needfor the new large-scale simulation, and do nothave used the ISSP Supercmputer.



References

[1] S. D. Yi, S. Onoda, N. Nishikawa, andJ. H. Han, Phys. Rev. B 80, 054416(2009).

Studies in recent years revealed that the spin
systems with the Dzyaloshinskii-Moriya(DM)interaction have various types of magnetic
structures such as the skyrmion lattice or thechiral magnetism [1]. Most of these studies
does not consider the eﬀect of the dipolar interaction between the magnetic moments ofthe spins, although it can not be removed.
We investigate the two-dimensional square lattice where the DM, the nearest-neighbor ferromagnetic, and the dipolar interactions compete with each other, by the Monte-Carlo(MC)simulation.
The Hamiltonian we consider in this study
is given as follows:∑

H = −J⟨i,j ⟩−D

Si · Sj∑{(S i × S i+ex )x + (S i × S i+ey )y

i 

+D′

∑  Si · Sj
i<j

−η

∑(

i



Siz

3
rij

)2

−H

−

3(S i · r ij )(S j · r ij ) 

∑

5
rij

Siz .

(1)

i

Here, r ij = t (xij , yij , 0) ≡ j − i means thediﬀerence between the position i and j. We
fixed the coupling constant of short range interaction as J = D = 1. In the simulation,
naive MC algorithm require O(N 2 ) simulation
time, and thus we adopt the O(N ) SCO algorithm [3].
As we can see in the figure 1, the skyrmion
lattice changes its structure depending mainly

323[2] M. Sasaki and F. Matsubara, J. Phys. Soc.
Jpn. 77, 024004 (2008).

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

(a)

(b)

(c)

Figure 1: Snapshots of the skyrmion latticephase at (a)D′ = 0.0, η = −1.2, H = 0.2,(b)D′ = 0.0, η = 0.0, H = 0.5, and (c)D′ =
0.6, η = 2.0, H = 1.5. The temperature ofeach one is fixed at T = 0.1. The skyrmions
make the square lattice at the cases (a) and(c), whereas they make the triangular latticeat (b). Note that direction of the square lattice
is diﬀerent in (a) and (c).

324

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

Investigation of Dynamics of a Supercooled Liquid andProteins using the Free-Energy Landscape
Takashi YOSHIDOME
Department of Applied Physics, Graduate School of Engineering, Tohoku University, Sendai
980-8579
Free-energy landscape (FEL) picture has

the two proteins at 400 K to observe unfolding

been used for understanding the properties ofprocess. Then the trajectories were analyzed. It

glass transition and protein folding [1,2]. The

was found that the experimental result was

free-energy landscape is defined as a function

successfully reproduced.

of a particle configuration [2]. In the present

The calculations were performed using thestudy, we have investigated the dynamics of a

L4cpu and F18acc in the system B, and L4cpu

Lennard-Jones system in a supercooled statein the system C. We used our custom-madeand the dynamics of two proteins upon

program [3] for computing the free-energy

unfolding process from the viewpoint of thelandscape of a Lennard-Jones system and for

FEL.

MC simulations. As for the study of proteinaunfolding, MD simulations were performed

crystallization process in a supercooled state.

using the program GROMACS [5]. The MD

The computation of free-energy landscape wastrajectories were analyzed using our custom-

conducted using the protocol we have proposed

made programs.

We

first

computed

the

FEL

for

[3]. The initial configuration was preparedusing a Monte-Carlo (MC) simulation at a

References

temperature above the freezing temperature. We[1]have successfully obtained a FEL upon

Applications to Clusters, Biomolecules andcrystallization.

Glasses (Cambridge University Press, 2004).

D.

Wales:

EnergyLandscapes:[2] T. Odagaki: J. Phys. Soc. Jpn. 86 (2017)

As for the protein dynamics upon unfolding
process, we focused on the two proteins (goat082001.

alpha-lactalbumin and human lysozyme). It has[3] T. Yoshidome, A. Yoshimori, and T.

experimentally shown that they have quite

Odagaki: Phys. Rev. E 76 (2007) 021506.

similar native structures but different structures[4] H.J.C. Berendsen, D. Van der Spoel, and R.

at the transition state. We first performed all-van Drunen, Comput. Phys. Commu. 91 (1995)

atom molecular dynamics (MD) simulations of

43.

325

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

Non-monotonic Behavior of the Binder Parameter inPotts Models
Hiroshi WATANABE
Department of Applied Physics and Physico-Informatics, Keio University
Yokohama, Kanagawa 223-8522, Japan
Satoshi MORITA, Yuichi MOTOYAMA, and Naoki KAWASHIMA
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277–8581

4.5
4

L = 32
L = 64
L =128(a)Binder Ratio

3.5
3
2.5
2
1.5
1
0.76

0.78

0.8

0.82

0.84

Temperature
1

(b)
0.8
Binder Ratio

The Binder parameter, a.k.a, the forth-order
cumulant, is a useful tool for analyzing the critical phenomena [1]. Since the scaling dimension of the Binder parameter is zero, its value
does not depend on the system size at the criticality. While the Binder parameter usually amonotonic function for the systems involving
continuous phase transition, it sometimes exhibits non-monotonic behavior. For example,
a hump appears in the Binder parameter near
the criticality of the Potts model. It exhibit
system size dependence, and therefore, it is
not trivial whether the universal scaling function itself has a hump or not. Additionally, it
is known that the amplitude of the hump depends on the definition of the order parameter.
However, the reason has not been clarified yet.

0.6
0.4
L = 32
L = 64
L =1280.2
00.76

To study the peculiar behavior of the
Binder parameter, we investigate the FortuinKastelyen graph representation. We find thatthe graph representation of the Binder parameter consists of two parts, namely, the low- andhigh-temperature components, and the peculiar behavior of the Binder parameter is causedby the low-temperature one. We also find thespin representation of the high-temperaturecomponents of the Binder parameter, which isa monotonic function. We study the J1 − −J2
frustrated Ising model [2] and find that the
non-monotonic behavior of the Binder parameter can be removed.

0.78

0.8

0.82

0.84

Temperature

Figure 1: The Binder parameters of the J1 −−J2 frustrated Ising model with (a) the conventional and (b) the proposed definitions.
The former exhibits humps while the latterdoes not.
[1] K. Binder, Z. Phys. B Cond. Matt. , vol.
43, pp. 119-140 (1981).
[2] S. Jin, A. Sen, and A. W. Sandvik, Phys.
Rev. Lett. vol. 108, pp. 2-6 (2012).

326

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

Simulating phonon cross section using density
functional theory
Naoki Murai
Materials and Life Science Division, J-PARC Center, Japan Atomic Energy Agency
Tokai, Ibaraki 319-1195, JapanThe studies of lattice dynamics, phonons,in crystalline materials are of fundamental importance as they provide key insights into theunderstanding of some relevant physical properties such as superconductivity and thermal
conductivity. Historically, the direct observation of phonon dispersion using inelastic neutron scattering (INS) by Brockhouse opened
the ﬁeld of phonon spectroscopy as an energyand momentum-resolved technique[1]. Since
then, especially with the development of modern time-of-ﬂight spectrometers at spallation
neutron sources[2], INS has emerged as a powerful tool that allows eﬃcient determination ofphonon dispersion. Meanwhile, as an alternative method, inelastic x-ray scattering (IXS)
has demonstrated its potential to precisely
measure phonon dispersion[3]. Thus, the ﬁeldof phonon spectroscopy has grown quickly overthe past the few decades, and the developmentof new instruments oﬀers a unique opportunity
to explore lattice dynamics in novel materials.
So far, we have been working on phonon dispersion measurements using the state-of-theart INS spectrometers at J-PARC. The present
project was carried out as a joint experimental
and numerical study of lattice dynamics properties in novel materials. In particular, thanks
to recent advances in density functional theory(DFT) calculations, it is now possible to perform ab initio phonon calculations whose only
input information is the chemical compositionof the materials.
In the present research project, we per-formed ab initio phonon calculations, with aparticular focus on unconventional superconductors and thermoelectric materials. Phonon
calculations were performed using the density
functional perturbation theory (DFPT)[4], as
implemented in quantum espresso code[5,6].
Figure.
1 shows representative results obtained for SnS, which is predicted
to be an attractive thermoelectric candidate.
Here, all calculations were performed using the plane-wave pseudopotential method
and the generalized gradient approximation(GGA) for the exchange correlation functional in the Perdew-Burke-Ernzerhof (PBE)
parametrization[7]. Pseudopotentials and energy cutoﬀs of plane-wave basis were chosen
according to the results of convergence tests
provided in the standard solid-state pseudopotential (SSSP) library[8], in which the precisionand performance of publicly available pseudopotential libraries are extensively tested tofacilitate the optimal choice of pseudopotentials. In the present work, we used ultrasoft
pseudopotentials from the GBRV library[9]with cutoﬀs of 70 Ry and 560 Ry for the expansion of the wave functions and charge densities, respectively. The Brillouin zone integration was performed over a 12 × 12 × 12 k mesh.
Dynamical matrices were computed on 6×6×6
meshes in q space, which were then interpolated to determine the full phonon dispersion.
Comparing the experimental phonon dispersion of SnS with DFPT calculations, we foundthat DFPT calculations are quite successful in327

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

predicting the experimental phonon frequencies. From the calculated phonon eigenvalues
and eigenvectors, we also computed the dynamical structure factor S(Q, ω), which provides an estimation of the phonon intensity. A
comparison with the neutron scattering data
revealed that the experimental phonon intensities of SnS were perfectly reproduced by thissimulation.
The present work thus serves to highlightthe potential to use DFPT phonon calculations
to simulate phonon cross sections as well as
phonon frequencies. A more detailed report is
now in course of preparation and will be published elsewhere.

References[1] B. N. Brockhouse and A. T. Stewart Phys.
Rev. 100, 756 (1995).
[2] K. Nakajima et al., Quantum Beam Sci.
1, 9 (2017).
[3] A. Q. R. Baron, High-Resolution Inelastic
X-Ray Scattering, in Synchrotron Light
Sources and Free-Electron Lasers: Accelerator Physics, Instrumentation andScience Applications, edited by E. J.
Jaeschke, S. Khan, J. R. Schneider, and J.
B. Hastings (Springer International Publishing, Cham, 2016) pp. 1643–1719, (seealso arXiv:1504.01098).
[4] S. Baroni et al., Rev. Mod. Phys. 73, 515(2001).
[5] P. Giannozzi et al., J. Phys. Condens.
Matter 21, 395502 (2009).
[6] P. Giannozzi et al., J. Phys. Condens.
Matter 29, 465901 (2017).
[7] J. P. Perdew, K. Burke, and M. Ernzerhof,
Phys. Rev. Lett. 77, 3865 (1996).
[8] G. Prandini et al., npj Computational
Materials 4, 72 (2018).
[9] K. F. Garrity et al., Computational Materials Science 81, 446 (2014).

Figure 1: Calculated phonon dispersion of SnS.
To facilitate the comparison between theory
and experiment, the dynamical structure factor, S(Q, ω), is weighted on the calculated
phonon dispersion curves.

328

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

Molecular dynamics simulations of water confined incarbon structures with the reactive force field
Tomoko MIZUGUCHI
Faculty of Materials Science and Engineering, Kyoto Institute of Technology
Matsugasaki, Sakyo-ku, Kyoto 606-8585

Many studies have been done on water in
confinement by both theoretical and experimental techniques and revealed that the properties of confined water may strongly diﬀer
from that of bulk water. It is clearly importantto properly evaluate the influence of confinement on water for various applications such as
ion channels or drug delivery. Until last year,
we had performed molecular dynamics (MD)simulations of confined water in nanoporous
silica with the reactive force field (ReaxFF)using the supercomputer at ISSP. We examined the influence of hydrophilic surface on
the physical properties of water, specifically,
the hydrogen bond structure and dynamics [1].
This year, I embarked on MD simulations using ReaxFF for water in a carbon nanotube asa model of a hydrophobic pore to clarify the
influence of interface on water.
ReaxFF allows us to incorporate bondbreakage and formation, which are particularly important on interface, in MD simulations with a lower computational cost thanquantum-mechanics-level calculations. I usedfour ReaxFF parameter sets: (i) developedby Monti et al. for biomolecules (protein2013) [2], (ii) by Zhang and van Duin to improve the weak interaction between hydrocarbon and water from protein-2013 force field(CHON-2017 weak) [3], (iii) by Yeon and van
Duin to simulate hydrolysis reactions at silica/water interface (Yeon2015) [4], and (iv)by Pitman and van Duin for clayzeolite com-posites (Pitman2012) [5]. I prepared a singlewalled carbon nanotube solvated in water andobserved the inflow and outflow of water into
the nanotube. The diameter of the carbon
nanotube (CNT) was set 2.7 nm for comparison with our previous results in the silica/water system. All simulations were conducted in the NPT ensemble (T = 300 K and
P = 1 atm) using the Nosé-Hoover thermostat
and barostat with LAMMPS pre-installed in
the supercomputer at ISSP. The ReaxFF MD
simulations were carried out with a time stepof 0.25 fs through the USER-REAXC packageof LAMMPS.
When the protein-2013 and Pitman2012
force fields were used, water flowed into the
CNT and filled it in 250 ps. In the 250-ps simulation using the Yeon2015 force field, water
flowed into the CNT somewhat. It seems toneed more time to fill in the CNT with water.
In fact, when I started the simulation with astructure where the CNT was filled with water,
the CNT remained filled with water. On theother hand, water went out of the CNT using
the CHON-2017 weak force field. Some experiments show that water spontaneously flows
inside the CNT, thus the CHON-2017 weakforce field seems not to be suitable for the
CNT/water system. But we should carefully
verify the result. I am now planning to carryout DFT calculations to compare with the simulation results using ReaxFF.

329

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

References[1] T. Mizuguchi, K. Hagita, S. Fujiwara andT. Yamada: Mo. Sim. 45 (2019) 1437.
[2] S. Monti, A. Corozzi, P. Fristrup, K. L.
Joshi, Y. K. Shin, P. Oelschlarger, A. C.
T. van Duin, and V. Barone: Phys. Chem.
Chem. Phys. 15 (2013) 15062.
[3] W. Zhang and A. C. T. van Duin: J. Phys.
Chem. B 122 (2018) 4083.
[4] J. Yeon and A. C. T. van Duin: J. Phys.
Chem. C 120 (2016) 305.
[5] M. C. Pitman and A. C. T. van Duin: J.
Am. Chem. Soc. 134 (2012) 3042.

330

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

Manipulation of electronic states with nanostructures: model
construction and its application
Toshikaze KARIYADOInternational Center for Materials Nanoarchitectonics,
National Institute for Materials Science, Namiki, Tsukuba, Ibaraki 305-0044
We focus on computing orbital magnetic

The improved program has been used tosusceptibility, which is relatively easy to

calculate the orbital magnetic susceptibility ofmeasure in experiment, to characterize novel

Sr3PbO, a Dirac material with antiperovskiteband structure such as Dirac/Weyl fermions.

crystalline structure. Temperature and carrier

There are some situations in which computation

number dependence of the susceptibility isof orbital magnetic susceptibility requires large

obtained using a tight-binding model derived

computing resources. One example of such afrom the ab-initio calculation. The improved

situation(i)

program enables us to have better results in low

nanostructured system, which is the main topic

temperature. It turned out that the susceptibilityof this project. Nanostructures lead to large

at low temperature is affected mainly by the

number of sites in a unit cell, which in turnchange in the temperature dependence of the

makes the computation hard. Another situation

chemical potential associated with the change

is when we handle with (ii) real materials.

in the density of states.

is

when

we

deal

with

Generically, taking more orbitals into accountAnother topic in this project is thegives better precision at the cost of computing

investigation on the bilayers of atomically thinpower. Furthermore, in order to characterizematerials with some mismatch showing moire

singular band structures like Dirac/Weyl cones,

pattern. A generic theory on the symmetry-fine sampling of k-points is inevitable.

based constraints on the interlayer coupling is

A useful strategy to improve the efficiency
ofcomputationisapplying

proposed [1], which is inspired by data sets ofmassivethe interlayer coupling in bilayer graphene

parallelization on k-point sampling, which had

obtained by systematic ab-initio calculations.

been done in the preceding project last year.
This year, we have made some corrections inReferences

memory handling to push up the maximum[1] T. Kariyado and A. Vishwanath, Phys. Rev.

system size allowed in the same facility.

Research 1 (2019) 033076.

331

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

Statistical Properties in a Branching Network
Satoshi YUKAWA, Kei TAKAYA, and Takeru WATANABE
Department of Earth and Space Science, Osaka University
Machikaneyama, Toyonaka, Osaka 560-0043

0.025
A = 2, arid
A=4
A=6
A=8
A = 10, humid
θ = 72°

0.02

0.015

N(θ)

In a river network, there are many interesting statistical laws, such as Hack’s law andHorton’s law. Recently another statistical law
was found in the distribution of bifurcation
angle.[1, 2] The distribution has a universal
characteristic angle depending on the aridityof the climate:[3] For the humid area, the angle is 2π/5 = 72◦ and, on the other hand, forthe highly arid area, it takes a smaller angle ofalmost 45◦ . The angle 2π/5 is theoretically explained for the Laplacian growth network, but
the smaller angle is not.
In this supercomputer project, we developa theoretical model to understand the aridity
dependence of the network based on the coupled map lattice model.[4] The model is constructed on a square lattice with water and soil
variables. Dynamical rules are too complicatedto be expressed here. The resultant networks
generated by computer simulation are shownin Figs. 1.

0.01

0.005

0
0

20

40

60

80

100

θ(°)

120

140

160

180

Figure 2: Bifurcation angle distribution forseveral aridity parameters A.
of the distribution decreases. There is, however, discontinuous behavior around the characteristic angle in the distribution. The origin
of such behavior is unclear. We need additional
computer simulations.

References[1] O. Devauchelle, A. P. Petroﬀ, H. F. Seybold, and D. H. Rothman: Proceedingsof the National Academy of Sciences 109(2012) 20832.

(a) humid case(b) arid case

Figure 1: Simulated network structures
Distribution of the bifurcation angle can beobtained from the simulated network. The result is presented in Fig. 2. We observe slightaridity dependence in the ﬁgure. Increasingaridity index, the peak at the characteristic angle becomes much larger. Besides the lower tail[2] S. Yukawa, T. Watanabe, and K. Hara:
Journal of the Physical Society of Japan88 (2019) 024901.
[3] H. Seybold, D. H. Rothman, and J. W.
Kirchner: Geophysical Research Letters 44(2017) 2272.
[4] B. Murray and C. Paola: Earth Surface
Processes and Landforms 22 (1997) 1001.

332

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

Numerical Study of Cell Guiding Mechanism in
Cellular Tissue
Katsuyoshi Matsushita
Department of Biological Science, Osaka University
Toyonaka, Osaka 563-0043

Cells have function in the organisms. In thedevelopment of organisms, cells are not always
born in the suitable place which the cells function. In these cases, cell migration play a significant role to transport the cells to the suitable place. As a prominent observation, there
is a case that the transport length of the cell
migration reaches mm scales much larger thanthe cell’s typical scale of 10µm [1]. In suchcases, the migration should be properly guided
by surrounding tissues to transport the cell for
its suitable position. The guiding mechanism
is still unclear and is expected to utilize chemical signals.

effective in the case of cell transport on tissue.
We show that the migration direction ofcells on the tissue aligns in the same direction
even when the cells are suspended. The order in cell migration is induced as a polarity
memory effect [6, 7]. The net transport velocity exhibits a non-linear transport and becomes a finite value above a threshold adhesion
strength. Furthermore, the transport distanceis longer than 102 times of the cell scale. This
result show potential of the heterophilic adhesion with polarity as the guiding mechanism ofcell transports.

As a candidate of the chemical signal, thechemical cell-cell adhesion is one of most powerful candidate. The possibility of the cell-celladhesion on the explanation of the long distance transport of cells is not sufficiently examined. In this work, we have investigatedthe cell migration on tissue which is guided by
the cell-cell adhesion. We employed the Cellular Potts model and examine cell-cell adhesion
effect on the cell migration in tissue [2].

References

We consider the model consists of suspended
transported cells on cellular tissue on the
Cellular Potts model[3] and one dimensionaltransport of the cells. The cells and tissue interacts with the heterophilic cell-cell adhesion,
As a guiding mechanism of the cell migration
direction, we consider the polarity of the cellcell l adhesion. The polarity is a hypothetical factor for guiding collective cell migrations
with long distance [4, 5] and is expected to be[1] C. Lois and A. Alvarez-Buylla, Science
264, 1145 (1994).
[2] K. Matsushita: submitted to PRE.
[3] F. Graner and J. Glazier, Phys. Rev. Lett.
69, 2013 (1992).
[4] K. Matsushita: Phys. Rev. E 97, 042413(2018).
[5] K. Matsushita: Phys. Rev. E 95, 032415(2017).
[6] K Matsushita, K Horibe, N Kamamoto,S Yabunaka, K Fujimoto, Proc. Sympo.
Simul. Traffic Flow 25, 21 (2019).
[7] Matsushita, K Horibe, N Kamamoto, K
Fujimoto, J. Phys. Soc. Jpn. 88, 103801(2019).

333

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

Magnetism in the System with the Multiple-Spin
Exchange Interactions
Chitoshi YASUDA
Department of Physics and Earth Sciences, Faculty of Science,
University of the Ryukyus, Okinawa 903-0213, Japan
Helium-3 (3 He) atoms adsorbed ongraphite form the commensurate solid layer
regarded as the spin-1/2 quantum spin
system on a triangular lattice at certain
density of 3 He. In the exchange process of
3 He atoms, there are not only an exchange
between two neighbor atoms but also thatamong three or more neighbor atoms. The
exchange integrals of three- and four-spin
cyclic exchange interactions are larger
than that of two-spin exchange interaction.
Furthermore, the comparison between thetheoretical and experimental studies showsthat eﬀects of the five- and six-spin exchangeinteractions are not negligible. Thus, in
relation to the solid 3 He, the multiple-spin
exchange (MSE) model on the triangularlattice has been extensively studied. Groundstates and thermodynamic properties of the
MSE model on the triangular lattice havebeen investigated using various methods,e.g., the mean-field approximation[1,2], theMonte Carlo simulations, the exact diagonalization, and the spin-wave theory[3,4].
In particular, the model has attracted
attention owing to the existence of novel
states such as the quantum spin liquid andthe spin-nematic state and the eﬀects offluctuations on ordered phases with chirality.
However, the study becomes diﬃcult for thecomplexity of the system with the various
interactions on the triangular lattice.
Recently, a new quantum spin liquid is experimentally observed in the solid 3 He layer.
This layer is a monolayer 3 He adsorbed ongraphite preplated with atomic layers of deuterium hydride. At low density of 3 He, it
was suggested that there exists the quantum spin liquid with novel dependences of334

the heat capacity and the magnetic susceptibility on temperature. In particular, the lowtemperature heat capacity is proportional to
2/3 power of temperature, which is an unprecedented behavior in the quantum spin
system. The details of this experimental
work are not yet clear, however, it is suggested that the solid 3 He layer is a honeycomb lattice.
We investigate the MSE model on thehoneycomb lattice using numerical methods,
which were applied in works of the MSE
model on the triangular lattice. In thepresent project, we concentrate the classical MSE model with the two-spin interaction and the six-body ring exchange interaction. The ground state is investigated using the conjugate gradient method. In thismethod, we need to prepare many random
spin states as the initial states because of
avoidance of the trap to a metastable state.
For large six-body ring exchange interactions,
we found the classical ground state with themulti-sublattice structure. We also advanceresearch on thermodynamic properties using
the exchange Monte Carlo simulations.

References[1] C. Yasuda, Y. Uchihira, K. Kubo, J.
Mag. Mag. Mater. 310, 1285 (2007).
[2] C. Yasuda, Y. Uchihira, S. Taira, K.
Kubo, J. Phys. Soc. Jpn. 87 104704(2018).
[3] C. Yasuda, D. Kinouchi, K. Kubo, J.
Phys. Soc. Jpn. 75, 104705 (2006).
[4] S. Taira, C. Yasuda, T. Momoi, K.
Kubo, J. Phys. Soc. Jpn. 88, 014701(2019).

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

Simulation Analysis of Collision of Two Filaments Driven by
Protein Molecular Motors
Kazuo SASAKI
Department of Applied Physics, Graduate School of Engineering, Tohoku University,
Sendai, 980-8579
It is known that protein filaments driven bysimulations of the particles for several particle

protein motors distributed on the glass surface

densities and parameters of the model, and theirform vortices when the density of filaments iscollective motion was analyzed from thehigh [1]. We studied this collective motion by

viewpointcomputer simulation. In this year, we focused on

percolation transition was found to occur as thethe following two systems. One is the system

particle density was increased if the interaction

consisting of two filaments and a lot of protein

between particles was week, while a new

motors [2]. The other is the system in which thephenomenon, the character of which is to be

filaments were simplified as point particles.

clarified, was observed in the case of strongof

percolation

theory[4]. AIn the first system, we studied the collisioninteraction. Formation of vortices was not

between two protein filaments, which is closely

observed, probably due to the simplicity of therelated to the emergence of the vortices. Fourmodel.

types of collisions were found experimentally

The calculations were performed using the[1]: parallel, anti-parallel, stop, and crossing; theL4cpu in both the systems B and C. We used our

probability of which type was observed strongly

custom-made program for the calculations.

depended on the incident angle between two
filaments. We simulated the collision of the two

References

filaments for various incident angles. We have[1] Y. Sumino, et al.: Nature 483 (2012) 448.

successfully reproduced the dependences of the[2] Y. Ishigure and T. Nitta: Langmuir 30 (2014)

probabilities on collision angle for parallel, anti-

12089.

parallel, and stop of collision types.

[3] F. Peruani, T. Klauss, A. Deutsch, and A.
Voss-Boehm: Phys. Rev. Lett. 106 (2011)

In the second system, we simplified themodel of protein filaments in accordance with

128101.

the previous study [3]: The filaments were[4] D. Stauffer and A. Aharony: Introduction To

modeled as self-propelled particles; and they

Percolation Theory, 2nd ed. (Taylor & Francis,

moved on a lattice. We performed Monte-Carlo

London, 1992).

335

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

Elucidation of the coercivity mechanism by statisticalphysical approaches
Masamichi Nishino
National Institute for Materials Science,
Tsukuba, Ibaraki 305-0044, JapanThe neodymium (Nd) permanent magnet,
Nd2 Fe14 B, which has strong coercivity, is anindispensable material in modern technologies, applied to information-storage devices,hybrid and electric vehicles, generators, etc. [1]This magnet has also another interesting feature. It exhibits a temperature-induced spinreorientation (SR) transition. The spin dynamics on the SR transition in the magnet has
not been clarified. The temperature dependence of the ferromagnetic resonance (FMR)
frequency is an open question.

quency above TR in connection to the temperature, magnetization along to the easy axis,
and the fluctuation of the magnetization along
to the hard axis (hard plane). We found that
this formula is a good description for overalltemperature region, and the last factor (fluctuation of transverse magnetization) in the formula is important for the qualitative nature of
the FMR frequency.
TR
200

In the present work we studied the temperature dependence of the FMR frequency at zeroexternal field. [2] The FMR spectrum I(f ) is
calculated by the auto-correlation function of
spins (power spectrum) for an atomistic modelof the Nd magnet. The temperature dependence of the FMR frequency fR is plotted byblue diamonds in Fig. 1.
We find that fR exhibits a drastic changearound TR (SR transition temperature) andnon monotonic temperature dependence. In
particular, it is found that fR ∼ 0 below TR .
This is totally diﬀerent from the dependenceof conventional magnets with a single uniaxial
anisotropy energy, in which a monotonic decrease of the FMR frequency is observed. We
clarified the mechanism for fR ∼ 0 below TR .
At zero temperature the state of the tilted spinalignment is stable, and the precession aroundthe ground-state (GS) easy axis (parallel to
the GS magnetic moments) does not occur.
We also derived the formula for the FMR fre-f [GHz]

150
100
50
0

0

200

400

600

800

1000

T [K]

Figure 1: Temperature dependences of fR(blue diamonds) for the atomistic Nd magnet
model..

References[1] S. Hirosawa, M. Nishino,
Miyashita, Adv. Nat. Sci:
Nanotechnol. 8 (2017) 013002.

and S.
Nanosci.

[2] M. Nishino, and S. Miyashita, Phys. Rev.
B 100 (2019) 020403(R).

336

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

Transition in the memory effect of paste
Akio NAKAHARA
College of Science and Technology, Nihon University
7-24-1 Narashino-dai, Funabashi, Chiba 274-8501
A dense colloidal suspension, called a paste,

reproduce the transition from memory of shear

remembers the direction of its motion, such as

to vibration, we find that we have to abandon thevibration and shear, and these memories in paste

assumption of simple shear flow and solve more

dominate the preferential direction for cracks torealistic flow motion under vibration.

propagate. When a paste remembers vibration,
cracks propagate in the direction perpendicular
to the direction of the vibration, while, when apaste remembers shear, cracks propagate along
the shear as is shown in Fig. 1 [1-2].
(a) Setup(b) Velocity field

Fig. 2: Direction of a simple shear applied to apaste in numerical simulations.

Fig. 1: Crack patterns induced by memory effect.
(a) Memory of vibration (b) Memory of shear.
To explain the formation of striped crack
pattern, we performed numerical simulations ofFig. 3: Elongated clusters along shear flow withshear motion of colloidal suspension using

its solid volume fraction of paste 3.3%.

LAMMPS. Colloidal particles attract each other

References

via Lennard-Jones potential and receives[1] L. Goehring, A. Nakahara, T. Dutta, S.

Stokes’s drag force from surrounding fluid.

Kitsunezaki and S. Tarafdar, Desiccation

When we assume a simple shear as shown in Fig.

Cracks and their Patterns (Wiley, 2015).

2, we obtain formation of elongated clusters[2] Japan patent JP 6565067, Material Mixing

along shear, as shown in Fig. 3, which dominates

Method, A. Nakahara, Y. Matsuo, M. Ito andpreferential direction for cracks to propagate. To

R. Yoneyama (2019)337

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

Numerical Study of One Dimensional Frustrated
Quantum Spin Systems
Kazuo HIDA
Professor Emeritus, Division of Material Science,
Graduate School of Science and Engineering,
Saitama University, Saitama, Saitama 338-8570

We investigate the ground-state phases ofspin-1 alternating-bond diamond chains described by the following Hamiltonian [1]:
H=

N [
∑
(1)
(1)
(1 + δ)S l τ l + (1 − δ)τ l S l+1

Hence, each ground state consists of an array ofn successive Tl = 1 pairs separated by singlet
dimers. This phase is called a DCn phase. If
no singlet dimers are present, the ground stateis equivalent to that of a mixed spin chain,

l=1
(2)

+ (1 + δ)S l τ l](1) (2)+ λτ l τ l ,
(1)

(2)

+ (1 − δ)τ l S l+1

H=
(1)

(2)

where S l , τ l , τ l are spin-1 operators. The
lattice structure is depicted in Fig. 1.
(1)

τl
1+δ

1−δ

λ
Sl
1+δ

Sl+1
1−δ
(2)

τl

Figure 1: Structure of a spin-1 alternatingbond diamond chain. S = τ (1) = τ (2) = 1.
The exchange parameters are denoted by λ,1 + δ, and 1 − δ, where λ and δ control thestrengths of frustration and bond alternation,
respectively.
Each ground state consists of an array ofspin clusters separated by singlet dimers owing to an inﬁnite number of local conservation
law,
(1)

∀l [T 2l , H] = 0, (T l ≡ τ l(2)+ τ l ).

(2)

A pair of spins with Tl = 0 is a singlet dimer.
For λ > 2, the state with Tl = 2 is not allowed.

338

N [
∑
(1 + δ)S l T l + (1 − δ)T l S l+1
l=1]+ λ(Tl (Tl + 1) − 4) (Tl = 1 or 2)(3)

with inﬁnite length. For strong frustration, we
ﬁnd a series of quantum phase transitions asin the case of alternating-bond mixed diamondchains with spins 1 and 1/2[2]. For intermediate frustration, we ﬁnd the nonmagnetic Haldane or dimer phases according to whether thebond alternation is weak or strong. For weak
frustration and weak bond alternation, we
ﬁnd the ferrimagnetic states with spontaneous
magnetizations m = 1/6 (F1/6 phase) and 1/3(F1/3 phase) per site. The F1/6 phase is accompanied by a spontaneous translational symmetry breakdown. This phase vanishes for strong
bond alternation[1]. The presence of a narrow
partial ferrimagnetic phase is suggested near
the point (λ, δ) = (λc , δc ) ≃ (1.0832, 0.2598).
The phase diagram determined by the numerical calculation for the Hamiltonian (3) using ﬁnite-size DMRG and inﬁnite-size DMRG
methods is shown in Fig. 2 except for the partial ferrimagnetic phase whose width is estimated to be so narrow that it is invisible inthis scale.

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

References[1] K. Hida, J. Phys. Soc. Jpn. 89, 024709(2020)[2] K. Hida, K. Takano, and H. Suzuki, J.
Phys. Soc. Jpn. 79, 044702 (2010).

1
(a)
δ

DC∞
Dimer
(λ,δ)
=(1.03261,0.42495)

F1/3
(λ,δ)
=(2.5338,0.2598)

0.5

(λ,δ)
=(1.0832,0.2598)
DC0
F1/6
DC∞
Haldane
0

1

DC1

2

3

λ

4

1
(b)
δ
DC∞
Dimer

0.5

(λ,δ)
=(1.03261,0.42495)

F1/3

(λ,δ)
=(1.0832,0.2598)
F1/6
DC∞
Haldane
0
0.9

1

1.1

λ

Figure 2: (a) Overall ground-state phase diagram. Small open circles are the phase boundaries between DCn and DC(n − 1) phases. Forn ≥ 4, these phase boundaries are only shownby solid lines to avoid complications. Other
phase boundaries are shown by the large open
circles and solid lines. The exact DC0-DC1
phase boundary is shown by the thick dotted line. The open square is the accumulation point of inﬁnite series of critical points at
δ = δc (≃ 0.2598). The ﬁlled squares are triple
points. (b) The enlarged phase diagram in thesmall λ region. The thick dotted line is the approximate relation for δ ≃ 1. (from Ref. [1])339

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

Construction of A New Method for Analysis of Crystallization
Mechanism from Amorphous Precursor Using Molecular
Dynamics Simulation Method and Unsupervised Machine
Learning
Hiroki NADA
National Institute of Advanced Industrial Science and Technology (AIST),
16-1 Onogawa, Tsukuba, Ibaraki 305-8569

crystallization from disordered systems in

The control of crystallization fromdisordered systems, such as amorphous

MD simulation studies.

crystalAn MD simulation of a supercooledengineering. However, our understanding

Lennard-Jones melt was performed usingof the mechanism of crystallization from

ISSP supercomputer. The methodology

disordered systems is still poor. This is

was utilized to evaluate the structural

because the structure of disorderedsimilarity between the melt and its

systems can hardly be elucidated by

crystalline form, which was regarded as

experimental means.

the geometrical similarity of the pair-

and

melt,

is

essential

to

Molecular dynamics (MD) are often

distribution function between them, at

helpful to elucidate the structure ofeach period during crystallization. An

disordered systems. A question is how toabrupt change of the structural similarity,

evaluate

which corresponded to crystal nucleation,thetime-evolution

of

the

structure in disordered systems during

wassuccessfully

detected.

Itiscrystallization properly.

anticipated this methodology contributes

Recently, a new methodology forto elucidation of the mechanism ofevaluation of the structural similarity

crystallization, even for multicomponent

betweensystems,

different

phases

using

the

such

as

crystallizationofdimensional reduction, which has beencalcium carbonate from its amorphousdeveloped as a method of unsupervised

precursor.

machine learning, is proposed [1]. The
purpose of this project is to demonstrate

Referencesthat this methodology contributes to[1] H. Nada: ACS Omega 3 (2018) 578.

elucidatingthemechanismof340

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

The analysis of the dynamical response to the opticalbistability
Tatsuhiko SHIRAIGreen Computing Systems Research Organization, Waseda University
Waseda-chou, Shinjuku-ku, Tokyo 162-0042

We have studied a dynamical response ofan optical bistable system to a time-periodic
driving ﬁeld. The optical bistability has beenfound in cavity systems with an external driving ﬁeld. As a function of laser intensity, there
is a ﬁnite interval with bistability of a hightransmission state (HTS) and a low transmission state (LTS). When the intensity of the
input laser is time-periodically modulated beyond the bistable regime, the system shows adynamical phase transition when the period of
the modulation changes. At large period, thesystem stays in a diﬀerent stable state depending on whether the laser intensity is increasingor decreasing. On the other hand, at small period, the system is kept in either the HTS or
the LTS. Here, we show that the Floquet dissipative map is useful to characterize the dynamical phase transition. The Floquet dissipative
map is given by a time-evolution operator ofa quantum master equation, and it characterizes the limit cycle and the decay rate of therelaxation dynamics. We found that the decayrate exhibits qualitatively diﬀerent system-size
dependence before and after the phase transition, and it shows a ﬁnite-size scaling of spinodal phenomena around the transition point.
This work was published to Physical Review
A [1].

operator L as a matrix. Naively, the number of elements in L increases exponentiallywith the number of spins N , which gives astrong restriction on N in the numerical simulation. However, we made use of the fact thatthe present model has a symmetry under exchange of spins. This symmetry reduces thenumber of non-zero elements in L to the orderof N 3 [2, 3, 4]. We have used this property andperformed simulations up to N = 100.
This research was supported by MEXT as
the “Exploratory Challenge on Post-K Computer project (Challenge of Basic ScienceExploring Extremes through Multi-Physics
and Multi-Scale Simulations).

References[1] T. Shirai, S. Todo, and S. Miyashita:
Phys. Rev. A 101 (2020) 013809.
[2] S. Sarkar, and J. S. Satchell: EPL 3(1987) 797.
[3] M. Gegg, and M. Richter: New J. Phys
18 (2016) 043037.
[4] T. Shirai, S. Todo, H. de Raedt, and
S. Miyashita: Phys. Rev. A 98 (2018)
043802.

In this work, we develop a numerical algorithm to eﬃciently solve the dynamics governed by the quantum master equation consisting of many spins. To do it numerically,
it is necessary to express the time evolution341

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

Theoretical study for superconductivity insuperstructured strongly correlated materials with
two-particle self-consistent approach
Kazutaka NISHIGUCHI
Graduate School of Science, Technology, and Innovation, Kobe University
1-1 Rokkodai-cho, Nada-ku, Kobe 657-8501, Japan– High-Tc superconductivity in multi-layer
cuprates –We have theoretically studied superconductivity in multi-layer cuprate, which is a series ofsuperstructured strongly correlated materials.
Here 1–5-layer Hubbard model has been considered as an effective model for multi-layer
cuprate superconductors. Then we have investigated the multi-layer superconductivity bysolving the linearized Eliashberg equation forthe superconducting gap function within the
two-particle self-consistent (TPSC) approach.
We have also developed the TPSC approach so
as to treat the multi-layer systems. Many degrees of freedom arising from the multi-layer
properties have been treated by the massive
parallelization in ISSP Supercomputer Systems.
– Possible superconductivity induced by spinorbit coupling in carrier doped Sr2 IrO4 –Possible superconductivity in in carrie- doped
iridium oxide insulator Sr2 IrO4 have studied
by a weak-coupling approach[1]. We havetreated t2g 3-orbital Hubbard model with alarge spin-orbit coupling (SOC) as an effectivemodel for carrie- doped Sr2 IrO4 . Numericallysolving the linearized Eliashberg equation forthe superconducting (SC) gap function withthe random phase approximation (RPA), we
have systematically examined both singlet andtriplet SC gap functions with possible pairing
symmetries and their parameter dependence.
For the realistic SOC λ and Hund’s coupling

J/U relevant to Sr2 IrO4 , namely, for a large
λ and small J/U region, we find that the
intra-band antiferromagnetic (AF) pseudospin
jeff = −l + s fluctuations favor a dx2 −y2 -wave
pseudospin jeff = 1/2 singlet pairing in theelectron-doping (Figure 1). We also find thatthe dx2 −y2 -wave pairing is more stabilized withincreasing the SOC and decreasing the Hund’s
coupling. Furthermore, we show for a small
λ and large J/U region that an s± -wave singlet pairing is favored in the hole-doped region(Figure 2). The origin of the s± -wave pairing
is due to the interband pair scattering arisingfrom the intra-orbital AF spin s fluctuations.
Although the possibility of a pseudospin triplet
pairing is considered, we find it always unfavorable for all parameters studied here. We
have also developed the multi-orbital and spindependent RPA so as to treat the multi-orbital
systems with spin-orbit coupling. Many degrees of freedom arising from the multi-orbitaland spin properties, namely, spin-orbit coupling, have been treated by the massive parallelization in ISSP Supercomputer Systems.
– Electronic structures of doped SrTiO3 as
H2 -generation semiconductor photocatalyst –We have performed the first-principles calculations based on the density functional theory toobtain electronic structures of doped SrTiO3
as a H2 -generation semiconductor photocatalyst. We have set the dopants to be a Rh
and Sb atom and replaced the Ti sites withthem in the 3 × 3 × 3 SrTiO3 supercell. The

342

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

0.9

0.8 0.7

0.3

0.2

0.4

0.5
0.6

0.9
0.8

0.6

0.7

0.5
0.6

Figure 1: Contour plot of the largest eigenvalue of the linearized Eliashberg equation
for the dx2 −y2 -wave pseudospin singlet pairingin the parameter space of the electron density n and the SOC λ for U = 0.8 eV, andJ/U = 0.05. The red (blue) region represents
the dx2 −y2 -wave pseudospin singlet SC (x-AF)
phase.

Figure 2: Contour plots of largest eigenvalueof the linearized Eliashberg equation for thes± -wave pseudospin singlet pairing in the parameter space of n and λ for U = 0.9 eV and
J/U = 0.3. The red region represents the s± wave pseudospin singlet SC phase.
[2] VASP (The Vienna Ab initio Simulation
Package): https://www.vasp.at

O-vacancies have also been considered as theintrinsic defect in the realistic experimental situations. The first-principles calculations havebeen performed by using VASP (The Vienna
Ab initio Simulation Package)[2], in which thehybrid functionals (HSE06) are adopted to reproduce the accurate band gap of SrTiO3 , andthe massive band and k-point parallelization
are used in ISSP Supercomputer Systems. The
numerical calculations have shown from thedensity of states that the Rh-doping produces
the impurity levels within the band gap as both
donor and acceptor levels. We have also observed that the Sb-doping works as electrondoping and then removes the acceptors in the
Rh+Sb-doping systems. In addition, the Ovacancies also works as electron-doping, whichleads SrTiO3 to the n-type semiconductors.

References[1] K. Nishiguchi, T. Shirakawa, H. Watanabe, R. Arita, and S. Yunoki: J. Phys.
Soc. Jpn. 88 (2019) 094701.

343

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

Study on Ising machines from a viewpoint ofstatistical mechanics
Shu TANAKA
Green Computing Systems Research Organization, Waseda University
27 Wasedacho, Shinjuku, Tokyo 162-0042

We have studied the performance of Isingmachines from a viewpoint of statistical mechanics. The Ising machines have been attracting attention in recent years as a computer thathas the potential to solve combinatorial optimization problems with high speed and accuracy. Very recently, quantum annealing machines, simulated annealing machines, and various other types of Ising machines have beendeveloped. In order to draw out the potential of the Ising machine, it is more important
to develop not only the hardware of the Isingmachine but also the algorithm and softwarein accordance with the established mathematical framework. We conducted research on thefollowing two topics on the basis of this trend.
(A) A guiding principle for embeddingin simulated-annealing based Ising machines
When we use the Ising machine, the costfunction and constraints of the combinatorial
optimization problem are first expressed in the
Hamiltonian of the Ising model or equivalentform. Herein, we refer to the Ising model as thelogical Ising model. Next, we represent the logical Ising model on the limitations of the network of spins in an Ising machine. The operation is called embedding. In embedding algorithms, there are hyperparameters that should
be adjusted to improve the computational performance, and how to adjust them is still unknown. In order to solve this question, we
have succeeded in constructing a method of
embedding algorithm based on the calculationof the correlation function by the transfer matrix. The results suggest an improvement in
the computational performance of simulated
annealing-based Ising machines against typicalexisting embedding algorithms [1].
This work was done in collaboration with
Dr. Tatsuhiko Shirai (Waseda Unviersity).
(B) The performance of black-box optimization using Ising machines from aviewpoint of statistical mechanical informatics
We have previously proposed a new methodof black-box optimization using the Ising machine [2]. The results suggest that this method
is superior to existing methods in the search
for metamaterials with excellent thermal radiation characteristics. The theoretical basis forthe superiority of this method is clarified fromthe perspective of statistical mechanical information.
This work was done in collaboration with
Dr. Ryo Tamura (National Institute for Materials Science, The University of Tokyo).

References[1] T. Shirai, S. Tanaka, and N. Togawa,submitted. preprint version is available via https://www.dropbox.com/s/
o8alkr8ch05ahtr/bare_jrnl_compsoc.
pdf[2] K. Kitai, J. Guo, S. Ju, S. Tanaka, K.
Tsuda, J. Shiomi, and R. Tamura: Phys.
Rev. Research 2 (2020) 013319.

344

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

Nonequilibrium phase transition and slow dynamics
in the dense hard sphere systems
Masaharu ISOBE
Nagoya Institute of Technology
Gokiso-cho, Showa-ku, Nagoya, 466-8555

The hard disk/sphere systems are one of
the simplest models, which have been actively
investigated in both equilibrium and nonequilibrium statistical physics. In this project,
we investigated non-equilibrium phenomenain the hard disk/sphere model system withmodern algorithms, especially for Event-Chain
Monte Carlo(ECMC) [1] and Event-Driven
Molecular Dynamics(EDMD) [2], where we
propose the “Hybrid Scheme”, namely, ECMC
for equilibration and EDMD for calculation of
dynamical properties [3, 4].
Dynamic facilitation and direct evidence ofvoid induced structural relaxations in glass formers:
On approaching glass transition, the viscosity of the supercooled liquid would be increased
rapidly. It often behaves as solid with disorder
structure within observation timescale. One of
the challenging issues is to understand whether
the essential properties of glass-forming materials are fundamentally thermodynamic or dynamic in origin. However, the microscopic origin of structural (so-called α) slow relaxationin deeply supercooled liquids at an atomic
scale has remained elusive due to the limitation of both electron microscopy experiments
and computer simulation. Although many theoretical models have been proposed and debated actively, one of a perspective that favors a dynamic origin is called Dynamic Facilitation (DF) theory [5, 6, 7]. DF theory was
known to explain a wide range of empiricallyobserved dynamical features of thermal super-cooled liquids and glasses. The applicability ofDF theory to athermal systems, i.e., systemsof hard particles where the relevant control parameter is pressure rather than temperature,under “super-compressed” conditions was investigated by employing novel eﬃcient algorithms [1, 2, 3]. The novel methodologies
allowed studying in detail the true equilibrium dynamics of high-density systems, whichis consistent with the predictions of DF theory. In such athermal molecular systems under“super-compressed” conditions — where what
is facilitated is the ability of the constituent
particles to structurally relax — giving rise
to correlated and cooperative dynamics, in amanner predicted by theory [8].
In collaboration with the experimentalgroup in Hong Kong, we resolve the coarsegrained trajectories obtained by colloidal experiments to conﬁrm DF theory. The results
undergo a transition from collective “creeping” to the string-like “hopping” motion on approaching the glass transition, where mobile ofparticles can be rearranged in string-like manners or compact geometries. Still, the compact
ones decomposed into connected sequences ofstring-like hopping motion. Furthermore, sucha string-like motion has an ultra-high returning hopping probability, which is proved to be
the leading contributor to the kinetic arrest forslow relaxation. One of the paradoxical issues
in DF theory is how a starter of string-like
hopping has suﬃcient room for the ﬁrst hop
in highly dense packing fraction. A void as a345

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

defect in the crystal structure can diﬀuse one
after another on the crystal lattice positions incase of dense mono-disperse systems. However,in a dense amorphous structure, there are few
voids with a size of particle diameter. Instead,
we found the quasi-voids distributed withina few neighboring shells. The suﬃcient freevolumes come from accumulated quasi-voids
within a local area as a soft spot might induce
the origin of string motion, where free volumes
transport along to the string-like hopping motions. These hypotheses are now carefully investigated by molecular dynamics simulation,
where the results would show direct evidenceof the scenario based on DF theory [9].

[9] C.-T. Yip, M. Isobe et al.: in preparation.

References[1] E. P. Bernard, W. Krauth, and D. B. Wilson: Phys. Rev. E 80 (2009) 056704.
[2] M. Isobe: Int. J. Mod. Phys. C 10 (1999)
1281.
[3] M. Isobe: Mol. Sim. 42 (2016) 1317.
[4] M. Isobe: “Advances in the Comutational Sciences — Proceedings of the Symposium in Honor of Dr Berni Alder’s
90th Birthday”, edited by Eric Shwegler, Brenda M. Rubenstein, and Stephen
B. Libby, World Scientiﬁc, Chapter 6(2017) 83.
[5] D. Chandler and J. P. Garrahan: Annu.
Rev. Phys. Chem. 61 (2010) 191.
[6] L. O. Hedges, R. L. Jack, J. P. Garrahan and D. Chandler: Science 323 (2009)1309.
[7] A. S. Keys, L. O. Hedges, J. P. Garrahan,
S. C. Glotzer, and D. Chandler, Phys.
Rev. X, 1 (2011) 021013.
[8] M. Isobe, A. S. Keys, D. Chandler, and J.
P. Garrahan: Phys. Rev. Lett. 117 (2016)
145701.

346

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

Glass phases induced by correlations
Chisa Hotta
Department of Basic Science, University of Tokyo
Meguro-ku, Komaba, Tokyo 153-8902

We studied the possibility of having a finite
temperature glass phases driven by the corelation eﬀect. Conventional spin glasses in realistic spatial dimensions, d = 2, 3 are realized in a series of Edwards-Anderson models
with classical interactions between spins, Jij ,which distribute at around Jij = 0. There, thedegree of disorder the only energy scale thatcompetes with the thermal fluctuation. Evenin such case the spin glass phases breaking thetime reversal symmetry is allowed to exist at
finite temperature only when d ≥ 3. We studied two diﬀerent types of models [1,2].
One is the classical vector spin Heisenberg
model on a pyrochlore lattice with additional
Ising-like lattice degrees of freedom, σi [1]; σi =±1 indicate the displacement of the lattice siteat each vertex, in and out of the tetrahedra,and the exchange interaction Jσi ,σj between
Heisenberg spins depends on the location of
displacements of sites that the spins live on.
Although there is no quenched randomness inthe system, the elastic energy term ∝ σi σj andthe spin exchange energy term Jσi ,σj Si · Sj arecorrelated with each other via Jσi ,σj and workas dynamical randomness to each other. The
two degrees of freedom finally undergo a simultaneous glass transition, which is detectedby the scaling analysis on the auto-correlation
functions[1].
The other one is the quantum transverse
Ising model on a triangular lattice with verysmall quenched bond randomness[2]. Thistime, the transverse field first works on thespin degrees of freedom and drives the system
to a KT phase, and successively to the clock-347

phase, where the three-sublattice magnetic ordering appears. Introducing a very small bondrandomness breaks the regular three-sublattice
structure of the clock phase into pieces of domains, cooporatively with the emergent glass
behavior of the spin degrees of freedom. The
spins and lattices(bond order) work to eachother as an emergent random fictitious field,
which is the origin of the finite temperature
spin glass in d = 2 system.
There are two common aspects between
these two cases. First, the model has a finite spin coupling constant J, and energy scale
represnting the quenched disorder is either absent or very small, which is contradictory to
the case of conventional spin glass whose energy is scaled by the quenched disorder. However, since these J’s are defined on frustrated
lattices, the interactions themselves are notenough to stabilize any kind of magnetic orderings. This allows the emergent energy scaleof disorder to dominate the low energy prooerties of the system. Second point is that theinterplay of the two diﬀerent species of degrees
of freedom can stabilize the glass phase much
more easily than had been considered before.
Each of the degrees of freedom have a tendency
to order and are on the verge of transition. The
correation with other degrees of freedom would
not allow them to have a regular order, andboth cooporatively goes into the glass phase.

References[1] K. Mitsumoto, C. Hotta, and H. Yoshino,
Phys. Rev. Lett. 124, 087201, (2020).
[2] C. Hotta, K. Ueda and M. Imada, in preparation.

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

Magnetic structures of multiple-Q states infrustrated itinerant magnets
Takashi UCHIDA
Hokkaido University of Science
4-1, 7-15, Maeda, Teine-ku, Sapporo 006-8585

Recently, the magnetic structures ofmultiple-Q states have attracted attention
since the topologically protected vortexlike magnetic structures such as magnetic
skyrmions emerge as 3Q multiple helical spin
density waves in the Dzyaloshinskii-Moriya(DM) interaction driven systems under magnetic ﬁeld. In addition to these systems, thefrustrated systems without the DM interaction have been found to reveal magnetic
skyrmions [1, 2]. In the present work, wehave applied the molecular spin dynamics(MSD) method [3] to the triangular-lattice
single-band Hubbard model without the
DM interaction to investigate the magnetic
structures of multiple-Q states realized aroundthe antiferromagnetic-ferromagnetic phaseboundary.
The MSD method is based on the functional integral method with the static approximation for the spin ﬂuctuation theories andthe isothermal molecular dynamics method.
The method allows us to ﬁnd automaticallythe magnetic structures of a large system withthousands of atoms in a unit cell at ﬁnite temperatures.
In the numerical calculations the most timeconsuming process is the magnetic force calculation at each time step, where the local electronic structures are calculated in real
space by means of the recursion method. We
have adopted the automatic parallel calculation scheme and found it to be eﬀective in saving both computing time and CPU resources.

We have performed the magnetic structure
calculations on the supercell with 20×20 triangular lattice, which is embedded in a large cluster consisting of 3 × 3 supercells, each of which
are connected by the periodic boundary condition and are subject to the self-consistent uniform eﬀective medium. Under zero magnetic
ﬁeld and the ﬁxed value of the temperature
T /t = 0.0005, we have changed the Coulombinteraction strength U/t = 5 ∼ 8 and the electron number n = 1.29 ∼ 1.36. The various
multiple-Q states such as 2Q, 4Q, 2Q+1Q were
found to be realized near the ferromagnetic
boundary. Of these, the 2Q states are accompanied by small spontaneous magnetization. It
is not clear, however, whether the calculated
multiple-Q states are topologically protectedstructures or not. This is partly because thatthe main Q state is often accompanied by small
satellite components, which make the wholemagnetic structure complex. In order to makeclear the magnetic structures of the multipleQ states, the magnetic structure calculations
with the site-dependent self-consistent eﬀective medium are now in progress.

References[1] T. Okubo, S. Chung, and H. Kawamura:Phys. Rev. Lett. 108 (2012) 017206[2] Y. kakehashi, D. Koja, T. Olonbayar, and H.
Miyagi: J. Phys. Soc. Jpn. 87 (2018) 094712.
[3] Y. Kakehashi, S. Akbar, and N. Kimura:Phys. Rev. B57 (1998) 8354.

348

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

Magnetization reversal process in classical spin systems
Taichi Hinokihara
Institute for Solid State Physics,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
We studied a reveal the magnetization reversal process by using the classical spin Heisenberg modelwith the dipole-dipole interactions (DDI).
Recently, theoretical studies for the coercivity based on the atomistic-scale spin model are increasing
due to the development of the first-principles calculations.
While these studies enable us to obtain material-dependent properties, the computational cost also
increases to simulate such the model.
The problem becomes severe for studying the coercivity of magnet because it strongly depends on
the long-range interaction, DDI [1,2].
In this project, we focused on the following two subjects: (i) develop efficient methods to simulate
this system. (ii) propose a coarse-graining method to construct a macro spin model from theatomistic scale system obtained from the first-principles calculations.
For the subject (i), we recently developed the novel efficient method for the simulation of the
dynamics of spins, called time quantified Monte Carlo with the stochastic cutoff method(TQMC+SCO).
As an example, the computational cost for the spin update is for the three-dimensional DDI systemeven when the system consists of the complicated lattice structure.
We analytically prove that the TQMC+SCO gives the same Fokker-Planck equation with thestochastic Landau-Lifshitz-Gilbert (s-LLG) equation, which is widely used for the simulation of spin
dynamics.
Namely, the ensemble average of the magnetization trajectory calculated by the TQMC+SCO is thesame as that by the s-LLG.
We also evaluated the validity of the present method by demonstrating the magnetization reversal
process and comparing it to the result obtained by the s-LLG.
Figure 1 shows three magnetization-trajectories of the total magnetization along z-axis: s-LLG
without DDI, s-LLG with DDI, and TQMC+SCO with DDI.
This simulation is performed in the 10×10×10 spin system with exchange coupling 𝐽 , the349Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

anisotropy 𝐾 = 0.1𝐽, the DDI 𝐷 = 0.05𝐽, and the temperature 𝑇 = 0.2𝐽.
We performed the magnetic field = 0.085𝐽 to drive the magnetization reversal process.
This figure clearly shows two results: (a) the DDI plays an essential role in the system with DDI, and(b) both s-LLG and TQMC+SCO gives the almost same magnetization trajectory.
We have submitted our paper, and it is now under review [3].
Besides, we recently found the other way of the bond-update process, which is suitable for the GPU
calculation.
We are now preparing this paper.
Each macro spin motion follows this landscape.
Next, the exchange coupling between themacro spins is selected to optimize the physicalproperties such as domain-all width.
We study the validity of our model bycomparing it to other physical properties suchas the total-magnetization and the susceptibility.
We are now preparing this paper.

References
Fig.

1:

Magnetization

reversalprocess[1] Taichi Hinokihara, Masamichi Nishino,

calculated by s-LLG without DDI (black), s-

Yuta Toga, and Seiji Miyashita: Phys. Rev. B

LLG with DDI (red), and TQMC+SCO with

97 104427 (2018).

DDI (blue).

[2] Yuta Toga, Masamichi Nihsino, Seiji
Miyashita, Takashi Miyake, and Akimasa

For the subject (ii), we constructed the macro-

Sakuma: Phys. Rev. B 98 054418 (2018)

pin Hamiltonian from the atomistic-scale spin

[3]

model of Nd2Fe14B in the following procedure.

Munetaka Sasaki, and Seiji Miyashita, arXiv:

First, we construct the free-energy landscape

1811.00237v2 (2019)

for the single macro spin from the atomisticscale spin model.

350

Taichi Hinokihara, Yuta Okuyama,

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

Quantum simulation for quantum many body systems using
noisy-intermediate-scale quantum computers
Keisuke Fujii
Graduate School of Engineering Science,
Osaka University, 1-3 Machikaneyama, Toyonaka, Osaka 560-8531
to describe molecular electronic wavefunctions

Noisy Intermediate-Scale Quantum (NISQ)
devices, or medium-scale quantum computersquantitatively.

Ourmethodmitigatesthewhose qubits are not fault-tolerant, are expected

problem by decomposing VQE calculation for ato be realized within a few years. As Google

large system into a lot of VQE calculations forshowed in 2019 [1], NISQ devices may be able

small systems. We applied this method toto perform some tasks faster than the state-of-various molecules to verify accuracy andthe-art classical computer. Such tasks should be

efficiency by numerically simulating the outputs

sufficiently simple and suitable for a quantumof quantum circuits. The largest system wecomputer, because of its incapability to remove

investigated was the Benzene molecule with theerrors in the computation process. Particularly,

cc-pVDZ basis sets. We note that it is unfeasible

calculations for quantum many-body systems

to perform the simulation of the original VQE

and quantum chemical systems are expected toalgorithm on classical computers as it requires

be promising applications of NISQ devices.

228 qubits for this system.

In this project, we have developed a method

Our numerical simulation is based on theto compute properties of quantum many-bodyfollowing programs: a fast quantum circuit

systems, especially molecules in quantum

simulator Qulacs [3] developed by our group;chemistry, with NISQ devices and performed

OpenFermion library [4] which can handle

numerical simulations of them by using the ISSP

mapping fermionic Hamiltonians to qubit

supercomputer. Concretely, we invented a

Hamiltonian; a quantum chemistry libraryperturbative approach based on the variational

PySCF [5]. Our code was highly parallelized

quantum eigensolver (VQE) algorithm [2],using MPI, as our perturbative VQE method canwhich is a quantum-classical hybrid algorithm tobe performed in an embarrassingly parallel

compute the ground state energy of a system.

fashion owing to its additive scheme.

One of the obstacles to utilize NISQ devices in

K.F. would like to thank Wataru Mizukami

quantum chemistry is that the number offor his help and collaboration throughout thisavailable qubits is limited and often not enoughproject.

351

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

References[1] F. Arute et al., Nature 574, 505-510 (2019).
[2] A. Peruzzo et al., Nat. Commun. 5, 4213
(2014).
[3] Qulacs, https://github.com/qulacs/qulacs[4] J. McClean et al., arXiv:1710.07629[5] Q. Sun et al., WIREs Comput. Mol. Sci. 8,
e1340 (2017).

352

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

Study on Complex Systems by Generalized-Ensemble Algorithms
Takuya HAYASHI, Daiki MATSUBARA, and Yuko OKAMOTO
Department of Physics, Graduate School of Science, Nagoya University
Furo-cho, Chikusa-ku, Nagoya, Aichi 464-8602
We had two projects regarding generalized-

per one water molecule [4] isensemble algorithms applied to complex

𝑆 = 0.815533 ± 0.000066.

systems.

This estimate is in good agreement with the

The density of states (DOS) is one of theresults of several other research groups. The

most important physical values in statisticalresults also imply that we need to develop amechanics.

new simulation model in order to estimate the

Recently,

weproposedasimulation protocol, REWL-MUCAREM [1],residual entropy with more accuracy.

in order to obtain the DOS with high accuracy
even in large systems by combining the
Replica-Exchange

Wang-Landau(REWL)

method [2] and the Multicanonical ReplicaExchange Method (MUCAREM) [3]. The
effectiveness of REWL-MUCAREM wasdemonstrated by using the 2-dimensional
Ising model in [1].

Figure 1: The estimates of residual entropy ofWe also applied the REWL-MUCAREM

ice Ih by several research groups.

protocol to the ice Ih system in order toestimate the residual entropy with highIn the second project, we have proposed

accuracy. The residual entropy of ice has

thebecome one of good examples to test the

Method (GCREM) [5], which is one of theefficiency of sampling algorithms. However,multidimensional replica-exchange methods

it seems that there is slight disagreement[6]. This method achieves a random walk inamong the estimates by various calculation

the temperature and chemical potential space,methods (see Fig.1). In this work, we give ouri.e., variables of grand canonical ensemble, inlatest estimate of the residual entropy by

one simulation. In this study, we applied

REWL-MUCAREM protocol.

GCREM to the LJ potential system andOur latest estimate of residual entropy

Grand

Canonical

Replica-Exchange

succeeded in drawing a phase diagram353

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

includinggas-liquid

first-order

phasetransition, second-order phase transition, andsupercritical fluid (see Fig.2). We also
calculated the critical point as (𝑇c∗ , 𝜇𝑐∗ ) =(1.18 ± 0.01, −3.34 ± 0.01). This is in good
agreement with the results of previous study

(a)

(b)

[7].
GCREM could sample various structuralphases in one simulation run (see Fig.3): It
was possible to sample the states of different
numbers of particles such as gas, liquid,
supercritical fluid, because the number of(c)particles in the grand canonical ensemble

Figure 3: Snapshots of the LJ particle systemvaries. By using the structure that definesobtained from a single GCREM run. Fig.3 (a)each of these states, structural analysis based

corresponds to Gas phase, (b) corresponds toon, e.g., radial distribution functions andLiquid phase, and (c) corresponds to supercluster formation, became possible.

critical phase.

References[1] T. Hayashi and Y. Okamoto, Phys. Rev. E.
100, 043304 (2019).
[2] T. Vogel, Y. W. Li, T. Wüst and D. P.
Landau, Phys. Rev. Lett. 110, 210603 (2013).
[3] Y. Sugita and Y. Okamoto, Chem. Phys.
Lett. 329, 261 (2000).
Figure 2:[4] T. Hayashi, C. Muguruma, and Y.

The phase diagram on reduced

temperature and chemical potential space.

Okamoto, in preparation.

The black line indicate 1st order phase[5] D. Matsubara and Y. Okamoto, J. Chem.

transition points and the edge is the 2nd order

Phys. 152, 194108 (2020).

phase transition point, i.e., critical point.

[6] Y. Sugita, A. Kitao, and Y. Okamoto, J.

Dotted lines indicate the critical temperature

Chem. Phys. 113, 6042–6051 (2000).

and chemical potential.

[7] N. B. Wilding, Am. J. Phys. 69, 11471155 (2001).

354

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

Development of a Coarse-Grained Force Field
Designed with a Polar Water Model
Wataru Shinoda
Department of Materials Chemistry, Nagoya University
Furo-cho, Chikusa-ku, Nagoya 464-8653

Coarse-grained (CG) approaches have been
adopted to investigate biological events relevant to the mesoscale. Many models have been
proposed and used for various target systems,
depending on the level of accuracy required in
the CG description. One of the promising CG
FFs, recently renamed as SPICA FF [1, 2, 3],accurately reproduces the experimental thermodynamic properties and structure (distribution function) obtained from AA-MD simulations using the CHARMM FF [4]. However,
the ability of the SPICA FF to predict the morphology of highly charged molecular complexes
in aqueous solution is still limited due to thelack of the polarity with the CG water modeland to the assumption of a high (background)dielectric constant. To overcome this problem,
we developed a CG FF, pSPICA, that is basedon a polar CG water model in order to precisely describe the molecular events that are
difficult to simulate using the SPICA FF with
non-polar water in lipid membrane systems,
such as water defects, membrane electroporation, and morphological changes in the charged
lipid self-assemblies. We follow the parameteroptimization scheme used in the SPICA FF,
namely, the parameters are optimized to reproduce several thermodynamic quantities measured experimentally [1, 2, 3].
In this study, we used the LAMMPS for CGMD simulations, and the simulations were significantly accelerated by GPU. In the pSPICA
FF, the polar CG water particle is composedof a pair of interacting sites, WO and WH

355

!"#

!$#
!" !#
$! "!

Figure 1: (a) Polar CG water, representing
three water molecules, developed in this study.
The blue bead, WO, has a van der Waals radius shown with a black circle. The red bead,WH, is connected to WO via a rigid bond andis positioned within the van der Waals radius ofWO. (b) The CG mapping of a DMPG lipid.
All-atom representation is given in stick and
CG representation is given in transparent particle.

(Fig. 1 (a)), and was designed to represent
three neighboring water molecules in aqueous
solution. The two sites, WO and WH, areconnected by a rigid bond, and have positive
and negative partial charges, respectively. The
choice of the CG sites for lipid molecules is thesame as that in SPICA FF (Fig. 1). The parameters of the CG sites were systematically
optimized to reproduce surface/interfacial tension, density, solvation or transfer free energy,
as well as distribution functions obtained fromall-atom molecular dynamics trajectory generated with the CHARMM FF, following theActivity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

scheme used in the SPICA FF. Lipid membranes simulated using the present CG FF
demonstrate reasonable membrane area andthickness, elasticity, and line tension, which
ensure that the simulated lipid membranes exhibit proper mesoscopic morphology. Using
the pSPICA FF with the supercomputer at
ISSP, we performed a large scale simulation
to investigate the morphology of aggregatesof the lipids dimyristoylphosphatidylglycerol(DMPG) in dilute solution. Thanks to thereasonable variation of line tension and elastic moduli of the lipid membranes dependingon salt concentrations in the pSPICA FF, themorphological changes in the lipid aggregates
showed good agreement with those observed inexperiments (Fig. 2 (a)). CG-MD simulationsof the membrane electroporation processes offlat lipid bilayers and vesicles were successfully
performed by the pSPICA FF. A water channel bridging the lipid bilayer was formed spontaneously in the presence of an electric field
generated by an imbalanced NaCl salt distribution, allowing ion transport across the lipid
bilayer (Fig. 2 (b)). The results were similar to the pore opening process observed using
AA-MD. Thus, the new CG FF enabled us toquantitatively investigate biological events, especially those involving water strings and defects in lipid membranes and ionic molecules.

(a)Salt-free

500 mM NaCl(b)

References[1] W. Shinoda, R. DeVane, and M. L. Klein:
Mol. Simul. 33 (2007) 27.
[2] W. Shinoda, R. DeVane, and M. L. Klein:
J. Phys. Chem. B 114 (2010) 6836.
[3] S. Seo, W. Shinoda: J. Chem. Theory
Comput. 15 (2019) 762.
[4] J. B. Klauda, R. M. Venable, J. A. Freites, J. W. O Connor, D. J. Tobias, C.
Mondragon-Ramirez, I. Vorobyov, A. D.
MacKerell, R. W. Pastor: J. Chem. Phys.
B 114 (2010) 7830.

356

Figure 2: (a) Snapshots of the initial aggregate composed of 5000 DMPG units and thefinally obtained aggregate in salf-free and 500
mM NaCl solution. (b) Snapshot of the electroporation process found in the simulation fora planar membrane system.

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

Screening of Substitutable Element in Rare-earth
Inter-metallic Compounds
Riki KOBAYASHI
Department of Physics and Earth Sciences, Faculty of Science
University of the Ryukyus, Nakagami, Okinawa 903-0213, Japan300

Total Energy of R5Si3 system200

c-axis/a-axis

Energy (eV)

250
150
100
50
0

Y La Ce Pr NdPmSmEuGd Tb DyHo Er TmYb Lu
Rare earth elements

Figure 1: Total energy of R5 Si3 (R = REEs)
system after structure relaxation.

Recently, we have succeeded in growing single crystals of Ce5 Si3 and Ce5 Ga2 Ge with thetetragonal Cr5 B3 -type crystal structure andfound that the Cr2-sites in the Cr5 B3 -type
crystal structure form the Shastry-SutherlandTetrahedral-Lattice (SSTL). SSTL has four geometrically frustrated two-dimensional layers
in a unit cell, called orthogonal-dimer systems,
and the ground state of Ce2-sites in Ce5 Si3
and Ce5 Ga2 Ge is a spin-dimer due to this geometric frustration. There are no other Cecompounds showing spin-dimers at the groundstate other than Ce5 Si3 and Ce5 Ga2 Ge, which
are important targets for the study of heavy
fermion systems with geometric frustration.
However, few rare-earth compounds with the
Cr5 B3 -type crystal structure have been reported, and this lack of targets is a bottle-

1.90
1.85
1.80
1.75
1.70
1.65
1.60
1.55
1.50

Aspect Ratio of R5Si3 crystal structures

Y La Ce Pr Nd PmSm Eu Gd Tb Dy Ho Er Tm Yb Lu
Rare earth elements

Figure 2: Aspect ratio of R5 Si3 (R = REEs)
crystal structures.

neck for the above research. In this study, thestructural stability of R5 Si3 (R=Rare earth elements) is investigated by comparing the total energy of R5 Si3 obtained by ﬁrst-principles
calculations to ﬁnd new Cr5 B3 -type materials.
The calculations were performed using theVienna Ab-initio Simulation Package (VASP),
which is capable of generalized gradient density approximation (GGA), and the K-points
and energy cut-oﬀ values were based on the
8 × 8 × 4 mesh and VASP recommendations,
respectively[1]. In this study, both quasiNewtonian and conjugate gradient algorithms
are used to calculate the structural relaxation,and it is conﬁrmed that there is no signiﬁcantdiﬀerence in the results between the two algorithms.
Figure 1 shows the total energy of R5 Si3357

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

from ﬁrst-principles calculations. Pm5 Si3 andSm5 Si3 were excluded from the graphs because
Pm is a radioactive isotope with a short halflife, and Sm5 Si3 did not converge in the SCF
calculation. The present calculations revealed
that the total energy is lower in the compounds containing Eu and Yb elements thanin the others, and it may be diﬃcult to obtaina Cr5 B3 -type crystal structure in these compounds.
Figure 2 shows the aspect ratio of the lattice parameter of R5 Si3 . Systematic changes in
the aspect ratios were observed for compounds
containing light REE from Y to Nd. This indicates that a speciﬁc axis of the unit cell may
shrink due to lanthanide contraction.

References[1] G. Kresse and J. Furthmuller, Phys. Rev.
B 54, 11 169 (1996).

358