GaN-MOVPEの表面反応機構：現在第一原理計算でできること
西永頌先生

押山さん、エピ成長の計算、できませんか？

いやあ無理無理！こんな感じでしょ。

1995年

大野英男さん（一緒に札幌すすきので飲んだくれてた）

え、押山さんって、できることしかしないんだ！

日本学術振興会未来開拓学術研究推進事業 (1996年－2000年）
押山（筑波大）、常行（東大）、塚田（東大）、中山（千葉大）、白石(NTT)、…
1

富岳semicon研究会 10/18-19, 2022

筑波大学計算科学研究センター： 計算機の方々とのおつきあい
宇川彰

 素粒子物理と計算機– QCDPAXによる格子量子色力学（Lattice-QCD）– 岩崎洋一、星野力、小柳義夫の諸先生方

 計算科学研究センター （1992年～）
物性科学にグランドチャレンジってなじまないでしょ！

前例がないくらいの大規模（1万原子）の第一原理（密度汎関数法）計算でもやってみるか！

RSDFT(Real Space Density Functional Theory)
プロジェクトExaFLOPS 202x --K
2011 T2K 2008 PACS-CS 2005
- 2007
CP-PACS 1996 -QCDPAX 1989
岩田潤一

2高橋大介

朴泰佑

佐藤三久

月一回のミーティング、宿題出し合い→ 2011年のゴードンベル賞

Lattice QCD

押山さん、物性でのグランドチャレンジって何ですか？
何かやってくれませんか。

1995 –U Tsukuba

Days of vectorizationand computations
1985 – 1995 NEC
1983 – 1984 IBM
myself
富岳semicon研究会 10/18-19, 2022

Clarification of Microscopic Mechanisms of GaN Epitaxial Growth

Atomic reactions in the step-flow growth of GaN–Decomposition and Diffusion of source-gas molecules on theterrace?
–What is the step structure?
–Incorporation process at the step edges?
–Unexpected dynamical features–K. M. Bui, J.-I. Iwata, Y. Kangawa, K. Shiraishi, Y. Shigeta, and A. Oshiyama, J. Phys. Chem. C 122, 24665 (2018)
K. M. Bui, M. Boero, K. Shiraishi and A. Oshiyama, Jpn. J. Appl. Phys. 59, SGGK04 (2020)
K. M. Bui, K. Shiraishi and A. Oshiyama, Appl. Surf. Sci. 557, 149542 (2021)
M. Boero, K. M. Bui, K. Shiraishi, K. Ishisone, Y. Kangawaa, and A. Oshiyama, Appl. Surf. Sci. 599, 153935 (2022)
K. M. Bui, K. Shiraishi and A. Oshiyama, Appl. Surf. Sci. submitted (2022)

3

富岳semicon研究会 10/18-19, 2022

MetalOrganic Vapor Phase Epitaxy (MOVPE) of GaN

Mechanism ofstep-flow epitaxial growth?
 Fates of arriving species on the terrace?
 Behavior near step edges

1. Reactionsin the gas phaseof GaN-MOVPE?

 TMGa is decomposedin the gas phase,becoming GaH3 or Ga←DFT cal by Shiraishi group
Sakakibara et al., JJAP 2021

 NH3 is not decomposed←Time-of-Flight mass-spectroscopy
more than 99% is NH3(Nitta-Amano)⇒ Growth surface is Ga-rich and NH3 is a key player!
4

富岳semicon研究会 10/18-19, 2022

Epitaxial Growth of Power Semiconductor GaN

DFT calculations has revealed:
Gas Source NH3 decomposed on the growth surface
Resulting NH diffusing to the surface step-edges

Car-Parrinello Molecular Dynamics
CPMD
CPMD @@ 13001300 KK

Ga adatomGa

NHs at the edge incorporated to the GaN film

Subsurface
NHowever, real growth surface seems to be aliquidat growth temperature 1300 K.
5Gadiffuses
adatom anddiffuses andGa adatom
a dimerforms a dimer forms/ is condensed
富岳semicon研究会 10/18-19, 2022

MetalOrganic Vapor Phase Epitaxy (MOVPE) of GaN
Mechanism of step-flow epitaxial growth?
 Fates of arriving species on the terrace?
 Behavior near step edges
2. Structural Identification of GaN growth surfaces?
ad gasad 1 gasad 1 gasFree energy of each surface structure G =Esurface − nGaµGa − nN µ N − nH µH

222Kusaba et al., JJAP 56, 070304 (2017)using 2x2 lateral cellGa-adatom

23Ga-H Ga-adatom surface with the coverage of 0.25
 H-adsorbed (3Ga-H) surface with the coverage of0.75 , both satisfying electron counting rule6

富岳semicon研究会 10/18-19, 2022

Mixed Surface of the Ga-adatom and the 3Ga-H Structural Motifs
Ga-adatom motif

DFT + Bayesian Optimization (AI)Kusaba Kangawa, Kuboyama, Oshiyama,
APL 120, 021602 (2022)

3Ga-H motif

lowest energy before

For 2 nm x 2nm area,
27C8 = 4, 686,825 caseseven with some constraints…….
7⇒ Growth Surface: Ga-adatom and 3Ga-H Motifs
The Key Player is first NH3

富岳semicon研究会 10/18-19, 2022

Surface Steps on Vicinal GaN (0001)
Mono-Bilayer Atomic Steps

Higher atomic steps are energetically unfavorable due to the lack of the edge reconstruction
Sawada, Iwata & Oshiyama: APL 104, 051605 (2014); PRB 93, 235421 (2016)

GaN

Step Formation
Energy8( µGa + µ N =EGaN )富岳semicon研究会 10/18-19, 2022

NH3 Adsorption on GaN (0001) Surface
Bistable Adsorption Structures for NH3(Meta)Stable Geometries

Energy Landscape between the Two Geometries

H

By Car-Parrinello MD (CPMD) with BlueMoon
N 2.2 A

GaGa adatom

1.9 A

Electron Density
weak Ga-Ga bond

DFT total energy

Free energy
Ga adatom

Ga adatomon the adatom site
9on the top-atom site
富岳semicon研究会 10/18-19, 2022

Decomposition of NH3 (NH3 → NH/GaN +H2 or NH/GaN + 2H on GaN )When assuming a reaction pathway, NH intervenes into

Energy profile during the reaction
from NH3 to NH

from here(Red Free) Energy

the Ga-Ga bondto there
Strongerbonds

 Small Barrier (energy cost to make H2 in gasphase ) 0.63 eV for the bond breaking reaction（free-energy gain of H2 in the gas phase is 2.1 eV)

Without postulation,
CPMD-with-BlueMoon simulations with theconstraint of reducing the hydrogen coordination
number with N

On-adatom position
0.5 eV12

3On-top-Gaposition

Coordination numberNH2NHGa adatom
HN

NH3 adsorbed ⇒ decomposed into the NH unit intervening into weak Ga-Ga bond10

富岳semicon研究会 10/18-19, 2022

Diffusion of the NH unit on the Terrace of GaN(0001) Surface
Stable NH Unit
Ga-adatom surface
Ga adatom
NSearch for the diffusion pathway byHyperplane Constraint Method

Also on the Mixed Surface
Ga-adatom Assisted Complex Diffusion

HNGa

Mixed (Ga-adatom+ Ga-H) surface
Ga adatom
HHNGa
N
From here
11To thereDiffusion barrier = 0.59 eVDiffusion barrier = 0.66 eV
富岳semicon研究会 10/18-19, 2022

Electron Density at Step Edges: Ga-Ga Weak Bonds[Ga Atoms near Step Edges]

GN Step
Electron Density near GN Step Edge

Ga-Ga bondsideview

Ga-Ga weak bondtopview

Ga2 Step
Ga-Ga bondsideview
topview

12

Ga-Ga weak bond:
Hot spot for the growth

富岳semicon研究会 10/18-19, 2022

NH Incorporation at GN Step Edges,,, and then Step-Flow Growth
NH Intervening in Ga-Ga Weak Bonds at GN Step Edges

 2 isolated NH in Ga-Ga bond on terrace -> 2 NH at step edge + 0.45 eV Exothermic Attraction at Step Edges
New Ga on Ga-rich Surface -> Incorporation of the new Ga with H2 DesorbedµH2 = −k BT ln[ gk BT / pς transς rotς vibr ]= - 2.1 eV

top view

Step Edge Proceeds
Possible Microscopic Process of Step-Flow Growth13

富岳semicon研究会 10/18-19, 2022

Heating up GaN(0001) by Car-Parrinello Molecular Dynamics
Sudden Increase of Diffusivity ofGa adatoms at growth temperature 1300KGa adatom
GaSubsurface NCPMD@ 1300
CPMD@ 1300K KTemperature (K)1300 K

Heat it up by CPMD
Time (ps)Thermal Motion aroundthe stable positions

14adatomdiffuses
Ga Gaadatom
diffusesandandforms
a dimeris condensed
forms
a dimer/ is/condensed
富岳semicon研究会 10/18-19, 2022

2-Dimensional Liquid Phase of Ga Adatom at Growth Temperature
Diffusion Coefficient D of Ga adatom
on GaN(0001)
1 t →∞ 1 Nv I (t ) v I (0) dt
D= ∫∑03
N I =1 Arrhenius behavior with 0.28 eV barrier
at low temperature Sudden increase in D at 1300 K by an order
15

Pair correlation function of Ga adatomon GaN(0001)

 At low temperature, peak at 6 Å, characteristic
to the adsorption structure At 1300K, peek at 2.6 Å, characteristic
to Ga-Ga bond
富岳semicon研究会 10/18-19, 2022

Epitaxial Growth on Single Monolayer Ga Liquid?
NHx Incorporation

16

2D Ga
liquidlayer

富岳semicon研究会 10/18-19, 2022

Neural-Network-Assisted Order N - OFDFT Scheme
Imoto, Imada, Oshiyama: Phys Rev Research 3, 033198 (2021)

Computational time measured for
SiC 576～4704-atom systems

DFT / CPMD calculations formillions of atoms are feasible
with moderate size ofcomputational resources

17

富岳semicon研究会 10/18-19, 2022

Further Efforts to Large-Scale Calculations: Order N Scheme Current scheme of DFT (Kohn-Sham Scheme) scales as cube of the target size N, O(N3) But solving Euler equation in DFT costs O(N) in itself, if we get rid of the orbitals {ϕi (r )} 1 2(;())−∇+rrvneff
 2 ϕi (r ) =ε i ϕi (r )
δ EXC [n]n(r′)
veff (r ) =vnucl (r ) + ∫dr ′ +| r − r′ |
δ n(r )

δ Ts [n]δ EH [n] δ EXC [n]µ+ vnucl (r ) ++=δ n( r )δ n(r)
δ n( r ) {ϕi (r )} is unnecessary if we know the functional form of the kinetic energy

Ts [n] = ∫ τ TF (r ) F [n] drEnhancement factor

 We rely on the Neural Network to develop Ts [n]Imoto, Imada, Oshiyama:
Phys Rev Research 3, 033198 (2021)
18

富岳semicon研究会 10/18-19, 2022

Ts [n] = ∫ τ

TF(r ) FNN(r ; n) drNt
δ TsKS [n]11  δ TsNN [n]=−()()
Lrr
Neural-Network Training∑ mm Nt m=1 2  δ nδnσ ( x)

Output from j-th neuron in l-th layer:z (jl ) = σ (l ) (a (jl ) )

2x

D(l )

l −1
Input a j is the output(l )a j = ∑ W jk(l ) zk(l −1)
from the previous layer:
k =0Descriptors

z1(0) = s = ∇n  2(3π 2 )1 3 n 4 3 Output

z2(0) = q = ∇ 2 n  4(3π 2 ) 2 3 n5 3 

DN −1
DN −2
NN 2(N )(N )( N −1)( N − 2)z1
F=
( s , q ) ∑ W1i σ  ∑ Wijσ  ∑ W jk σ ( ⋅⋅⋅)  =j 0 =(In the present computation, # hidden layers=3,k0=i 0=
DN# neurons per layer =3)

Prepare n( r ) , s ( r ) and q (r ) from the training data (in carbon diamond).
(l )Construct F as functionals of n, q and s with the NN parameters Wij .
Train the NN (optimize NN parameters Wij(l ) ) so as to minimize the cost function L
19

富岳semicon研究会 10/18-19, 2022

Structural Propertiesof 11 solids

20

富岳semicon研究会 10/18-19, 2022