Finite Temperature Eﬀect on Structural Changes of Single
Crystalline Graphite under Plasma Irradiation
Seiki Saito1 , Atsushi M. Ito2 , Arimichi Takayama2 , and Hiroaki Nakamura1,2
1

Department of Energy Engineering and Science, Graduate school of Engineering, Nagoya University, Furo-cho,Chikusa-ku, Nagoya 464-01, Japan
2
Fundamental Physics Simulation Research Division, National Institute for Fusion Science, 322-6 Oroshi-cho, Toki
509-5292, Japan

1

Introduction

Carbon materials are promising candidates for divertor
plates in nuclear fusion reactors. Divertor plates of a nuclear fusion reactor are bombarded with hydrogen plasmas.
The hydrogen plasmas erode the carbon divertor plates,
yielding H2 and other hydrocarbon molecules such as CHx
and C2 Hx , which are undesirable impurities in plasma conﬁnement experiments. To understand the nature of chemical and physical interactions between hydrogen plasmas
and the divertor plates, it is important to clarify the elementary processes of reaction.
Under bombardment of plasmas whose thermal energies
are higher than several eV, it is expected that the structures of target materials dynamically changes. For example, when single crystalline graphite is bombarded by hydrogen plasmas, the surface of the graphite becomes hydrogenated amorphous carbon by the chemical and physicalreaction of hydrogen atoms. This structural change aﬀects
the processes of retention, reﬂection and sputtering.
In our previous work [1], hydrogen injection into single crystalline graphite was performed by binary-collisionapproximation-based (BCA) simulation and a hybrid simulation of BCA and molecular dynamics (MD) simulation.
As a result, an interesting phenomenon called ”channeling” was observed in the case of crystalline structure. Because of the channeling eﬀect, incident atoms with speciﬁcincident angle corresponding to crystalline structure movestraightly without causing large angle scattering in the target material. This eﬀect is not negligible in retention process because it makes incident atoms to reach deeper position in the target material. From the simulation, it wasalso found that the channeling eﬀect gradually disappears
because of the amorphousization of the target material.
The mean depth of incident atoms in the case where thechanneling eﬀect exists may strongly depend on the temperature of the target materials because the vibration ofatoms of target material causes dechanneling by hitting
and kicking the projectile away from the channeling path.
The time evolution of the bonding states of the target
material under plasma irradiation was also investigated in
our previous work. In that result, almost 30% of sp content
ratio was obtained at the maximum. In our other previous
work [2], we formed amorphous carbon in a manner ofdeposition on a diamond substrate by MD simulation. In
the simulation of deposition, sp2 rich amorphous carbon
was obtained. In that case, the percentage of sp carbons isless than 10%. We infer the reason of the diﬀerence of sp
content ratio in the amorphous carbon materials formedby plasma irradiation and deposition process is that thesimulation of plasma irradiation was executed under thecondition where the temperature of target material was
set to 0 K. It is expected that more sp2 or sp3 carbons aregenerated if the simulation is executed in the case where
the thermal vibration exists, because the sp carbons may
immediately make bonds with other carbons by annealing
process.
In this paper, therefore, to investigate the eﬀect of theﬁnite temperature of the target material, we perform hydrogen injection into a single crystalline graphite whichhas ﬁnite temperature by taking the thermal vibration ofatoms on Debye model in BCA simulation and attaching
Langevin thermostat to the MD simulation box in BCAMD hybrid simulation.
2BCA-MD Hybrid Simulation [3]To explain the hybrid simulation of BCA and MD more
concretely, let’s consider a hydrogen atom injection into agraphite material (Fig. 1). An injection kinetic energy ofthe hydrogen atom is set to 1 keV, which is higher thanEth . While the kinetic energy of the hydrogen atom ishigher than Eth , trajectories of the hydrogen atom and thesurrounding carbon atoms are calculated by BCA simualtion. The hydrogen atom dissipates its kinetic energy byinteracting with carbon atoms in graphite and then thekinetic energy becomes lower than Eth . At that moment,
MD simulation starts to calculate the motions of the hydrogen and the carbon atoms instead of BCA simulation.
Some of carbon atoms which collide with the incident hydrogen are kicked out on the way of BCA simulation. Themotions of the kicked carbons in the cascade process aresolved by the same procedure as the incident hydrogen.
Unlike BCA simulation, MD simulation cannot treat alarge scale material. To save the computation time, cubic simulation boxes (MD simulation boxes) whose side
is 34 Å long are picked up for MD simulation. The motion of equation for all atoms in the MD simulation boxes
is solved. After each injection, the structural relaxationof all atoms in the MD simulation boxes are executed tominimize their potential energies.
33.1

Method for Finite Temperature
Thermal Vibration in BCA Simulation

To investigate the temperature dependence of the channeling eﬀect, the thermal vibration of target atoms in BCA

hydrogenThe position of the target atom is shifted by W (r) at eachbinary collision.

BCA

3.2

Thermostat in MD simulation

Langevin thermostat is attached to the MD simulation
box to keep the material temperature the preset temperature. Langevin thermostat controls the temperature ofatoms by applying the random external force. The equation of motion for the atoms are described as follows:

MDMDMDṙ i =34Fig. 1. Schematic diagram of the BCA-MD hybrid simulation of hydrogen injection into a graphite.[3]∂H∂H, ṗi = −− [γpi + Ri (t)]∂pi∂r i(7)

where, r i and pi are the position and the momentum ofthe i-th atom, respectively. H is the Hamiltonian of thesystem. γ is the viscous damping coeﬃcient. The random
force Ri is described as the white noise as follows:⟨Ri (t)⟩ = 0, ⟨Riα (t1 )Riβ (t2 )⟩ = 2kB T δαβ δ(t1 − t2 ) (8)simulation is taken into account.
According to the Debye model, the internal energy of asolid in a thermal equilibrium state at the temperature Tis written as∫ ∞ϵdϵ,(1)
U=D(ϵ) ϵ/k T
B
e−1−∞
where kB is the Boltzmann constant, ϵ is the energy ofphonon, D(ϵ) is the phonon density of states described asfollows:
23V
0 ≤ ϵ ≤ kB TD
 2π2 (h̄c)3 ϵ ,D(ϵ) =(2) 0,ϵ < 0, kB TD < ϵwhere c is the sound velocity in the solid, V is the volume,h̄ is the Planck constant devided by 2π. TD is the Debyetemperature described asTD =h̄ckB(6π 2 NV)1/3,

11√dr,π (2ϵ/k)2 − r24

Research Approach

Hydrogen atom injection into single a crystalline
graphite material, which has ﬁnite temperature by taking
the thermal vibration of atoms on Debye model in BCA
simulation and attaching Langevin thermostat to the MD
simulation box, is performed by BCA-MD hybrid simulation. Hydrogen atoms are injected in succession, so thatthe structure of the target material gradually changes intime evolution. The temperature dependence of the timeevolution of the density of incident atoms in target material is extracted from the simulation results. The temperature dependence of the time evolution of the bonding
states of target material is also investigated.

(3)
Acknowledgement

where N is the number of atoms which constitute the solid.
The probability F (ϵ)dϵ that the average energy per atomu ≡ U/N in a range from ϵ to ϵ + dϵ is written as3l3ϵ30 ≤ ϵ ≤ kB TD 2π2 (h̄c)3 eϵ/kB T −1 dϵ,F (ϵ)dϵ =0,ϵ < 0, kB TD < ϵ(4)where l ≡ (V /N )1/3 is the mean interatomic distance.
Harmonic oscillation is assumed for the thermal vibration of atoms. In the case of harmonic oscillation with theenergy ϵ in classical mechanics, the probability W (ϵ, r)drthat an atom is observed in a range from r to r + dr is
described as
W (ϵ, r)dr =where, T is the preset temperature, δ is the Dirac delta
function, α and β denote the x-, y-, or z-component. Thebracket ⟨ ⟩ denotes a long-time average.

(5)

where k is the spring constant. The probability distribution function W (r) which is integrated for the energy range
from 0 to ∞ is described as follows:∫ kB TD
ϵ31
√dϵ,
(6)
W (r) = Aϵ/kT −1
Be(2ϵ/k)2 − r20
where A is the normalization constant.
The spring constant k is obtained by MD simulation.

Numerical simulations were carried out using thePlasma Simulator at the National Institute for Fusion
Science, Japan. The present study was supported by
the National Institutes of Natural Sciences as part ofthe Forming Bases for Interdisciplinary and International
Research through Cooperation Across Fields of Studyand Collaborative Research Programs (NIFS09KEIN0091
and NIFS10KEIN0161), a Grant-in-Aid for JSPS Fellows(No. 23-340), and a Grant-in-Aid for Scientiﬁc Research(No. 19055005) from the Ministry of Education, Culture,Sports, Science, and Technology, Japan.
References[1] S. Saito, A. M. Ito, A. Takayama, and H. Nakamura,
Structural Change of Single-Crystalline Graphite underPlasma Irradiation, submitted to Jpn. J. Appl. Phys.
[2] S. Saito, A. M. Ito, A. Takayama, and H. Nakamura,Anisotropic Bond Orientation of Amorphous Carbon byDeposition, Jpn. J. Appl. Phys., 51, 01AC05, 2012.
[3] S. Saito, A. M. Ito, A. Takayama, T. Kenmotsu, and H.
Nakamura, Hybrid Simulation between Molecular Dynamics and Binary Collision Approximation Codes for Hydrogen injection into Carbon Materials, J. Nucl. Mater, 415,
S208-S211, 2011.