Quantum-Classical Molecular Dynamics Simulation ofSynthesis of Boron-Nitride Nano-Structures at High
Temperature, High-Pressure Plasmas
Predrag S. Krstića,b and Longtao HanaaInstitute for Advanced Computational Science, Stony Brook University, New York, USA
bTheoretiK

Supported by US DoE, Basic Energy Sciences, through Princeton Plasma Physics Laboratory

4th MoD-PMI, Toki-Gifu, June 20191

Collaboration
This research was funded by BES
for nanosynthesis in volume of ahot, high pressure plasma through
PPPL

Prof. Roberto CarCarbon
Dr. Igor KaganovichPlasma theory

Dr. Yevgeny Raitses
experimentDr Longtao Han
Grad student (at the time)

Dr. Yoshifumi
Nishimura
Waseda U.

Prof. Hiromi Nakai

Computing facilities of IACS,XSEDE, DOE and Nagoya U.

Prof. Stephan Irle,Nagoya U. and ORNL

2

Why Boron-Nitride?
•Boron electronic config.
1s22s22p

B(5) and N(3) are immediate neighbors of C(4)in periodic table, BN(8) structure isisoelectronic with a bonded carbon structure
CC(8)Like graphene

Like carbon

Nitrogen electronic config
1s22s22p3

Like CNTCarbon 1s22s22p2
Boron “borrows” one p orbital fromnitrogen and then both are like carbon:The carbon story repeats

All identical structures
to carbon
Like carbon
Exceptional properties, similar to CNT and graphene; many applications
Unlike graphene: Difference in electronegativites (B 2.04; N 3.04) forms ionic bond; 5.9 eV homo-lumo gap :3
insulator , still good thermal conductor; piezoelectric properties

Why B and N form BN structures (and not B-B, N-N, B-B-N-N,…?)• B-N bond is energetically favorable (~4 eV) than B-B bond (3.1 eV) and N-N
single bond (~1.7 eV).
• Even-numbered rings, made of B and N alterations, are energetically favorable.
Odd numbered rings appear as defects.
Each atom in h-BN neighbored by three other atoms.
e.g. when N in center, if one B replaced by N, only single N-N bond can be formed.
Similarly in c-BN, B (or N) is coordinating with 4 atoms of other type.

c-BN

Self-organization of mixture of Bs and Ns will favor creation of BN structures
4t MoD-PMI, Toki-Gifu, June 2019
Oku T, Narita I and Nishiwaki A Mater. Manuf. Processes,
2004, 19, 1215-39.

4

Then: What were we doing here? How?
• Computationally studied synthesis of boron nitride nanostructures from a hot,high pressure gas containing various combinations of precursors B, N, and H.
(typically expected in plasma synthesis)• Why? To elucidate detailed dynamics of self-assembly mechanisms into clusters
and aggregates at various temperatures. understanding the chemical dynamics andpredicting the optimal conditions concerning temperature and gas content.
• Achieved by using quantum-classical molecular dynamics (QCMD) simulations
based on DFTB quantum chemical potential, implemented in the DC-DFTB-K code
(H, Nakai, Waseda U., Tokyo) to study large systems containing about 1300 atoms(10,000 basis functions) in a canonical NVT ensembles for 1 ns, for T’s 1500-6000 K.

Why in plasma?
Usually used process in nanosynthesis is Chemical Vapor Deposition:
Temperature < 1000K, deposition on a substrate, various pressure.

However, experiments on nanosynthesis in hot plasma volume(torch or arc discharge, no substrate) at high pressure showed
higher yield of production and defect-free synthesis
Kim et al, 2014: BNNT by feeding h-BN powder, along with N2 and H2 in the torch (8000K).
Hypothesis: BN dissociate, B agglomerate, H prevent N atoms to associate to N2: Final reaction some BNH precursors withB clusters.
Zetl et al, Plasma torch in N2, B feedstock, a few atm, BNNT, BN cages, BN flakes.
Hypothesis: N2 dissociate, react with B clusters

We have tried by theory and simulations to understand this nonequilibrium process
and suggest the optimal plasma conditions for the synthesis
4t MoD-PMI, Toki-Gifu, June 20196

Classical or quantum-mechanical approach?
Using Classical MD, varying various potentials and other parameters,
enabling EEM to include dynamics of charging has not resulted inorganized but rather to amorphous BN structures
Inspiration came from Y. Ohta, J. Comp. Chemistry, 2016Showed that can form small nanocages using QCMD with DFTB

10 ns
1 fs step
Varying TOnce we applied QM: The solution for thesynthesis in the plasma volume opened
4t MoD-PMI, Toki-Gifu, June 20197

Why quantum-mechanics works here?
Preliminary and ad hoc Hypothesis: Because of the Quantum Entanglement
What is quantum entanglement?
Quantum state of one atom cannot be described independently of allother particles in the system (that is being self-organized in our case),even if they are at large mutual distances. All particle correlated.
Quantum measurement: Whenever dump the coordinates of the particles
to the disk (each 4 ps)
The entanglement hypothesis of our self-organized systems is in an initial stageof study!! M. Zack did initial research more than 15 years ago.
4t MoD-PMI, Toki-Gifu, June 20198Thermodynamics of the process of self-organization
The change in Gibbs free energy

ΔG = ΔH-TΔS<0 fo the spontaneous process of synthesis.
Depends on temperature:
When ΔS < 0 and ΔH < 0, the process will be spontaneous at low temperatures and nonspontaneous at high temperatures.
A decrease in system entropy can only occur spontaneously if the entropy changeof the surroundings is both positive in sign and has a larger magnitude than theentropy change of the system.
The increase in entropy of the surroundings is accomplished via heat transfer from the system to thesurroundings
Big problem is change of the number of particles due to both implantation and loss
Non-equilibrium thermodynamics
IF, Beograd, June 2017

Methods
Quantum-Classical Molecular Dynamics
QCMD

4t MoD-PMI, Toki-Gifu, June 201910

About 2 orders ofmagnitude faster than
DFT

DFTB+

Many groups. Most versatile
version of DFTB.
Developed and maintainedby the Bremen group for
Computational Material
Science
4t MoD-PMI, Toki-Gifu, June 201911

B. Aradi, B. Hourahine, and Th. Frauenheim. DFTB+, a sparse matrix-based implementation of the DFTB
method, J. Phys. Chem. A, 111 5678 (2007).

4t MoD-PMI, Toki-Gifu, June 201912

Significant step forward for QCMD is DC-DFTB-K of the group ofprof. Hiromi Nakai in Waseda U, Tokyo!
Divide-and-Conquer DFTB: This software is capable, without a significant loss of accuracy!!,
to treat larger systems (thousands of atoms) in an acceptable real time.
DC-DFTB-K
Basic References:[1] H. Nishizawa, Y. Nishimura, M. Kobayashi, S. Irle, and H. Nakai, Three Pillars for Achieving Quantum MechanicalMolecular Dynamics Simulations of Huge Systems: Divide-and-Conquer, Density-Functional Tight-Binding, and MassivelyParallel Computation, J. Comp. Chem 37, 1983 (2016).
[2] H. Nakai, A. Wibawa Skti, and Y. Nishimura, Divide-and-Conquer-Type Density-Functional Tight-Binding Molecular
Dynamics Simulations of Proton Diffusion in a Bulk Water System, J. Phys. Chem, B 120, 217 (2016).

1) DC method => linear-scaling computation with respect to the sys size;2) DFTB theory, which avoids the explicit two-electron integral evaluation while keeping
reasonable accuracy,
3) Massively parallel computation.
The computational efficiency assessed including systems up to one million atoms.
4t MoD-PMI, Toki-Gifu, June 201913
This method is reaching N-scaling

RESULTS

4t MoD-PMI, Toki-Gifu, June 201914

Nanocage formation

Nanoflake formation

Formation of a BN nanocage by selforganization of BN diatomic molecules.
Snapshots are taken at t=0, t=264 ps and
t=480 ps, respectively. The implantation
rate is derived by linear regression. (a)Number of atoms and (b) number ofrings are plotted as a function of time.

Evolution of B36 boron cluster
irradiated by NH3 molecules into theflat nanostructure. Three snapshots aretaken at t=0, t=800 ps and t=2260 ps,
respectively. (a) Number of atoms and(b) numberIF,ofringsJune
areplotted as aBeograd,
2017
function of time.

BNNT growth

Growth of BNNT from a zig-zag template.
(a) Increased number of atoms in BNNT. (b)Grown BNNT after 5.8 ns. (c) Evolution ofring types. Growth of BNNT is completed by
self-organization of a BN side chains, which
migrated and accumulated at the top endof the tube. T and flux critical qunatities.

Synthesis of BN fullerene from BN dimers• Stockfeed by of BN dimers at high temperature andconcentration (pressure).

BN chains:

1. branch
2. aggregate
T = 2000K
4ps/BN moleculeP ~ 327 torr

4t MoD-PMI, Toki-Gifu, June 201916

Growth of a BNNT from a zig-zag template• Side chains formed by BNmolecules.
• BNNT grows by self-organizationsof the BN side chains• migrate and accumulate at the topof the tube by proper T• Lowering flux of feedstock provides
T=2500K
extra time for migration of side
24 ps/moleculechainsp~65 torr
4t MoD-PMI, Toki-Gifu, June 2019

Template

Product17

Need to break the symmetry to build a BNNT

4t MoD-PMI, Toki-Gifu, June 201918

Hydrogen rich feedstock leads to the 2D structures• NH3 and HBNH
• on boron cluster

Cage opens withedge saturated byhydrogen.

4t MoD-PMI, Toki-Gifu, June 2019

NH319

Growth of c-BN by particle irradiation at ambientpressure

4t MoD-PMI, Toki-Gifu, June 201920

Surface tension induces the high internal pressure• Surface effects become dominant
when cluster gets smaller.
• Young-Laplace equation:

Pi = 2γ / r• 𝛾𝛾 ~ 3𝐽𝐽/𝑚𝑚2 , 𝑟𝑟 ~ 1 𝑛𝑛𝑛𝑛 ⟶ P ~ 6 GPa• At 2000K, 6GPa is sufficient to bringc-BN to stable phase.
• H changes the hybridization of B
and N atoms, catalyze the transition
to c-BN.
4t MoD-PMI, Toki-Gifu, June 201921

How do we achieve synthesis of the BN
nanostructures form the “soup” of precursors?
Systems: Typically 1300 atoms, enclosed in a box of lateral dimension ~10 nm, with 3D
periodicity. This resulted in particle density of n~1021 cm-3 , i.e. initial pressure p=nkT~
1.4x104xT (Pa), which means 270 bars at 2000K. This large initial pressure promotes reactive
collisions frequency, accelerating the self-assembly process during the short time of 1 ns. The
pressure decreases since n decreases by agglomeration. All system ran for 1 ns, with time-stepof 1 fs.
Basic Toll: Quantum mechanics, Quantum-classical molecular dynamics, using DFTB, utilized
through the new code, DC-DFTB-K, developed at Waseda University, Tokyo, by Prof. Nakai andcollaborators. This code accelerated standard DFTB calculation by close to an order ofmagnitude (for our systems). Still, each point (defined by a precursor and a temperature)required about one week of calculation

Finding optimal precursors and temperature for the BNsynthesis from a hot, high pressure gas

4t MoD-PMI, Toki-Gifu, June 201923

Precursor systems and major products
11 configurations of precursors containing B and N, at 1500K-6000K
Simulations done by DC-DFTB-K
11 precursor systems
GroupA

Subgroup
A1
A2
B1

B
B2

Major productsContentSymbol

1/2B 1/2N
1/3B 2/3N
2/3B 1/3N
1/4B 3/4N
BN
1/3B 1/3N 1/3H
1/4B 1/2N 1/4H
1/2B 1/4N 1/4H
4/5BN 1/5H
HBNH
B3N3H6

B-N
B-2N
2B-N
B-3N
BN
BNH/111
BNH/121
BNH/211
BN-H
HBNH
Borazine
4t MoD-PMI, Toki-Gifu, June 2019

Initial pressure ~100 atm: increase the collision probability

24

Ratings of precursor systems for optimal synthesis of sp2 BNstructures

Probability of synthesis of h-BN structure form precursor X at T is product of probability of hexagon formation attime t and probability of formation of B into sp2 ideal hybrid• PSQ: parameter of the synthesis
qualityP( X , T ) = Phex ( X , T ) PspB2 ( X , T )N hex
Phex = max
N hexmaxhex

N=P

Bsp 2 + N sp 22B
sp 2=Bsp 2NB+ 2( Bsp 3 + N sp 3 )• Best choices:• 2000K• BN, BN-H, Borazine and HBNH
4t MoD-PMI, Toki-Gifu, June 201925

Stages of structure evolution for BN precursor system• BN molecules form BN chains (sp);• Chains branch (sp, sp2);• Branched chains join into various sp2 BN structures

0 ns0.03 ns

0.13 ns

1 ns
26

Evolution of B & N atoms hybridization• Sharp increase of sp atoms➔ fast creation of BN chains• Gradual increase of sp2
atoms with decrease of spatoms➔ chain branching andcreation of h-BN structures
4t MoD-PMI, Toki-Gifu, June 201927

Borazine at 2000K

4t MoD-PMI, Toki-Gifu, June 201928

Stages of structure evolution for borazine precursor system• Borazine molecules needs torelease Hs before assembly.
• Higher temperature:flakes• Promotes dissociation ofborazine and release of H
• Reduces ratio of H in final
product products

flakesflakes + cBN

4t MoD-PMI, Toki-Gifu, June 201929

Mechanisms of BNNT self-assembly from gas phaseprecursors
1. Merge of patches (a)• BN chains branch and join
• BN cluster connects withpatches of BN flakes (BNf)• Merge and reorganization

2. Self-folding of BNf (b)3. Coalesce andreorganization of closed
BN cages (c)4t MoD-PMI, Toki-Gifu, June 201930

What novelties have
we reported? Summary
Summary• BN-containing molecules are efficient feedstock for building various BN nanostructures attemperatures and pressures provided by electric arc and torch discharge.
• Hydrogen saturates the dangling bonds at the edges and stimulate creation of nanoflakes, andprevents formation of nano-onions• The atomic hydrogen catalyzes hybridization of B and N atoms from sp2 to sp3, leads to formationof nano-size c-BN at ambient P.
• Compound precursors containing B-N bond (BN, BN-H, Borazine and HBNH) are optimal forsynthesis of sp2 BN structures; optimal temperature is 2000 K.
• Understood some mechanisms of creation of BNNTs in a hot high-P plasma. Self organizationof BN diatomic molecules

Details of this talk in:

Han et al, Nanotechnology Lett. (IOP) 28, 07LT01 (2017)
Krstic et al, J. Phys. Chem. C 122, 936 (2018)
Krstic et al, Chem. Sci. 9, 3803 (2018)31

Thank you!

4t MoD-PMI, Toki-Gifu, June 201932

How to define the synthesis quality?

Hexagons of alternating B -N bonds basic structural “brick” for all ordered BN nanostructures.
1) Need counting of the hexagons (form the ring counting algorithm, adapted form TINKER package)
This number at time t is Nhex(t) – only hexagons with 3B and 3N counted.
2) At any time t, maximum number Nmax of hexagons formed from both trivalent (sp2) and tertravalent (sp3) hybridized
B and N atoms, taking only those B and N that are bound to their complement.
max
hex

N=

Bsp 2 + N sp 223) Probability of hexagon formation at time t is Phex =+ 2( Bsp 3 + N sp 3 )

N hex
max
N hex

Only hexagons in clusters with more than 10 atoms taken into account (in case of borazine 13, to avoid counting of theborazine hexagons.
4) We want h-BN based structures. To obtain probability of h-BN purity in the structures, multiply by probability of formationof B into the sp2 ideal hybrid
Bsp 2BPsp 2 =NB

Finally, probability of synthesis of h-BN structure form precursor X at T isP( X , T ) = Phex ( X , T ) P ( X , T )B
sp 2

PSQ = parameter of the synthesis quality

Toward root growth mechanism for BNNT

• Our simulations indicate that the root growth mechanism may not valid for BNNT:• N atoms at the nanotube root can dissolve by diffusion on boron cluster• Incorporation of feedstock happens at the top end of BNNT

4t MoD-PMI, Toki-Gifu, June 201934