Edge Formation Mechanism of Single-Walled Carbon Nanotube Revealed by
  Molecular Dynamics Simulations based on a Neural Network Potential
Despite the high potential for applications utilizing the unique properties
of single-walled carbon nanotubes (SWCNTs), the applications remain challenging
due to the difficulty of synthesizing SWCNTs with a specific chirality. To
elucidate the mechanisms that determine their chirality during growth,
intensive efforts have been devoted to classical molecular dynamics (MD)
simulations. However, the mechanism of chirality determination has not been
fully clarified, which can partly be attributed to the limited accuracy of
empirical interatomic potentials in reproducing the behavior of carbon and
metal atoms. In this work, we develop a neural network potential (NNP) for
carbon-metal system to accurately describe the SWCNT growth, and perform MD
simulations of SWCNT growth using the NNP. The MD simulations demonstrate the
process of defect-free chirality-definable SWCNT growth with dynamic
rearrangement of edge configurations, and the consistency between the
probability of edge configuration appearance and the entropy-driven edge
stability predicted from the nanotube-catalyst interfacial energy. It is also
shown that the edge defect formation is induced by vacancy and suppressed by
vacancy healing through adatom diffusion on the SWCNT edges. These results
provide insights into the edge formation thermodynamics and kinetics of SWCNTs,
an important clue to the chirality-controlled synthesis of SWCNTs.
Ikuma Kohata, Ryo Yoshikawa, Kaoru Hisama, Christophe Bichara, Keigo Otsuka, Shigeo Maruyama
# Catalyzed Single-Walled Carbon Nanotube Growth

###### Abstract

Classical molecular dynamics (MD) simulations have been performed to elucidate the mechanism of single-walled carbon nanotube (SWCNT) growth. However, discussing the chirality has been challenging due to topological defects in simulated SWCNTs. Recently, the application of neural network to interatomic potentials has been actively studied and such interatomic potentials are called neural network potentials (NNPs). NNPs have a better ability for approximate functions and can predict complex systems' energy more accurately than conventional interatomic potentials. In this study, we developed an NNP to describe SWCNT growth more accurately. We successfully demonstrated defect-free and chirality-definable growth of SWCNTs on Fe nanoparticles by MD simulation. Furthermore, it was revealed that edge vacancies of SWCNTs caused defect formation and adatom diffusion contributed to vacancy healing, preventing defect formation.
Introduction

Single-walled carbon nanotubes (SWCNTs) [1] have great potential as a next-generation material for nanoscale devices such as transistors [2], microprocessors [3], transparent conductive films [4] and solar cells [5] because of their excellent electrical, thermal, and optical properties [6, 7]. One of the most important properties of SWCNTs is that the chirality determines whether they are metallic or semiconducting as well as their bandgap [8, 9, 10]. Therefore, chirality control is important for the realization of high-end SWCNT-based devices. However, the method of chirality-controlled synthesis has not been fully established and therefore further elucidation of the growth mechanism is indispensable.

In the past few decades, a number of classical molecular dynamics (MD) simulations were carried out to reveal the growth mechanism of SWCNTs and demonstrated a series of the growth process from the nucleation to the sidewall elongation [11, 12, 13, 14, 15, 16, 17]. Although these studies provided insights into the basic growth mechanism of SWCNTs, the discussion of the chirality was still difficult because the simulated SWCNTs had a number of defects and the chirality could not be assigned. Instead, partial chirality was defined and discussed in several studies. Neyts et al. performed hybrid atomistic simulations of MD and force biased Monte Carlo based on the reactive force field (ReaxFF). They defined local chirality by the diameter and the angle of the hexagon to the axial direction of SWCNTs, and demonstrated the initial growth process of SWCNTs and the chirality change from zigzag to armchair [18, 19]. However, the continuous growth of SWCNTs with a particular chirality was not observed. Although Yoshikawa et al. assigned the chirality to SWCNTs by MD simulation with Tersoff-type interatomic potential [20], but the chiralities of the observed SWCNTs were only zigzag and near-zigzag, which was contrary to the experimental reports of near-armchair preferential growth [21, 22, 23, 24].

Currently, there are two possible reasons why it is difficult to simulate SWCNT growth without any intramolecular chirality change by defects, as reported in the experimental study [25]: first, time scale is too short for defect healing, and second, energetics of SWCNTs including defects is not accurately described by interatomic potentials. In order to solve these problems and to reproduce defect-free growth, it is necessary to perform MD simulations in a long time scale using an accurate interatomic potential.
Recently, the application of artificial neural network to interatomic potentials has been actively studied to overcome the limitations of conventional interatomic potentials and such interatomic potentials are called neural network potentials (NNPs). Since the idea of symmetry functions to represent high-dimensional potential energy surfaces (PESs) with translational and rotational invariance was proposed by Behler and Parrinello [26], a number of approaches have been developed to construct NNPs in the field of materials science. Among the NNPs, the deep potential [27] is one of the most popular methods because of its high accuracy and high computational efficiency. In this study, we developed an NNP based on the deep potential to simulate SWCNT growth more realistically, and reproduced defect-free SWCNT growth by MD simulation.

## 2 Method

### Development of NNP

The growth process of SWCNTs consists of the dynamic formation of carbon rings on the catalyst surface by recombination of inter-carbon bonds. To reproduce the growth of SWCNTs by MD simulations, the inter-atomic potential should predict the energy of various morphologies of carbon including unstable structures. To collect dataset structures covering a wide range of PESs, we sampled three types of atomic structures: disordered structures, defective graphene and CNT growth structures. The disordered structures were obtained by randomly placing atoms in a computational cell of 8-20 A per side, heating them to 10000 K by classical MD simulation using a Tersoff-type potential [28], and then quenching them to 1000-5000 K (Figure S1a). The defective graphene structures were collected by preparing graphene in a rhombic supercell of 17.24\(\times\)17.24\(\times\)10.0 A, randomly displacing carbon atoms by 0-4 A, and relaxing by MD simulations using Tersoff-type potential at 1000 K (Figure S1b). The CNT growth structures were sampled from MD simulations of SWCNT growth on Fe\({}_{60}\), Fe\({}_{80}\), Fe\({}_{100}\) and Fe\({}_{120}\) using Tersoff-type potential and NNP under development (Figure S1c). The cell size of the CNT growth structure was set to 15.0\(\times\)15.0\(\times\)15.0 A so that the length of the vacuum space is \(\sim\)10 A. The number of atomic configurations included in the disordered, defective graphene and CNT growth datasets are 21744, 16871 and 5776, respectively. Each dataset was randomly separated into three subsets, with 10 % as
the validation set, 10 % as the testing set and the remainder as the training set.

We performed spin-polarized density functional theory (DFT) calculations using the projector-augmented wave (PAW) method [29], Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) for exchange-correlation functional [30, 31], with 400 eV plane wave energy cutoff, and the k-point mesh for the Monkhorst-Pack grid Brillouin zone sampling (N\({}_{\text{k}}\)1, N\({}_{\text{k}}\)2, N\({}_{\text{k}}\)3) = \((3,3,3)\) for the disordered structures and \((1,1,1)\) for the defective graphene and CNT growth structures, as implemented in Vienna ab initio simulation package (VASP 5.4.4) [32, 33, 34, 35]. The convergence criterion for the electronic self-consistent loop was set to 10\({}^{-4}\) eV.

Deep potential [27], one of the high dimensional NNPs, was used for training the NNP. The machine learning descriptor was constructed from both angular and radial atomic configurations as implemented in DeePMD-kit [36]. The cutoff radius is \(4.0\) A. The smoothing function was set between 3.5-4.0 A. For the embedding net, the sizes of the hidden layers from the input end to the output end are 25, 50 and 100, respectively. For the fitting net, three hidden layers with 240 neurons were used. Parameter optimization was performed in 200000 steps. The learning rate is 0.001 at the beginning and decays exponentially every 5000 steps to 3.51\(\times 10^{-8}\) at the end of training. Comparisons of the energies and the forces between DFT and NNP are shown in Figure S2. The MAEs of the energies are 16.9 meV/atom and 16.1 meV/atom for the training set and testing set, and the MAEs of the forces are 0.268 eV/A and 0.248 eV/A for the training set and testing set, respectively.

### SWCNT growth simulation

MD simulations were performed using LAMMPS [37]. The time step for numerical integration was 0.5 fs. The temperature of Fe nanoparticles was controlled at \(T\) using the Nose-Hoover thermostat [38, 39]. The simulation box was a cubic cell with the size of 20\(\times\)20\(\times\)20 nm. Periodic boundary condition was applied for each direction of the axis. At the beginning of the simulation, we placed Fe\({}_{60-120}\) nanoparticles at the center of the cell. After relaxation at constant temperature for 2 ns, we supplied C atoms at random positions in the cell controlling the density of C atoms in the gas phase at 8 atoms/cell, which determines the rate at which carbon is supplied to the catalyst nanoparticles. The initial velocity vectors of C atoms \(\mathbf{v}\) were randomly
determined, assuming that the temperature of the system is \(T\) and the magnitude of the initial vector \(|\mathbf{v}|\) satisfies \(|\mathbf{v}|=\sqrt{3k_{b}T/m}\). To describe the energy barrier of the carbon feedstock decomposition, a Lennard-Jones potential was applied between non-covalent C atoms, as in previous SWCNT growth simulations [20, 28]. The visualization of the atomic structures and common neighbor analysis were performed using OVITO [40].

The developed views are made by expressing the coordinates of the SWCNT by the cylindrical coordinates \((r,\theta,z)\) and projecting them onto the x-z plane using the equation \(x=r\times\theta\), \(y=0\), \(z=z\). The chirality \((n,m)\) of each SWCNT was determined by the positive combination of the basis vectors of the developed SWCNT, \(a_{1}\) and \(a_{2}\), according to the original definition of the chirality.

## 3 Validation of NNP for Fe catalyzed SWCNT growth
overestimated the defect formation energy. This indicates that the stability of graphene was overestimated, which may have prevented cap lift-off by pentagon formation during the initial stages of SWCNT growth.

### Edge energy of graphene

The energetics of graphene edge has been considered to be an important factor of SWCNT chirality selectivity [43, 44]. We validated chiral angle dependence of the edge energy of graphene (Figure 2). To calculate

Figure 1: Comparison of defect formation energy of graphene between DFT and interatomic potentials.
the edge energy of graphene by the NNP, we modeled graphene nanoribbons (GNRs) in the cell with one periodic direction and vacuum space of 15 A in other two non-periodic directions, and relaxed them with conjugate gradient method. The edge energy \(\gamma\) is defined as

\[\gamma=(E_{total}-N_{C}E_{C}-N_{e}\mu)/2L, \tag{2}\]

where \(E_{total}\) is the total energy of the supercell, \(N_{C}\) is the total number of carbon atoms in the supercell, \(E_{C}\) is the energy of carbon atom in graphene lattice, \(N_{e}\) is the number of atoms of the terminating element, \(\mu\) is the chemical potential of the terminating element and \(L\) is the periodic length of graphene nanoribbon as defined in a previous theoretical study [43]. We computed the edge energy of (2,1), (3,1), (3,2) and (4,1) GNRs plotted as a function of chiral angle, as shown in Figure 2. For comparison, we also plotted the theoretical curve expressed as

\[\gamma^{\prime}(\chi)=2\gamma^{\prime}_{A}\sin{(\chi)}+2\gamma^{\prime}_{Z}\sin {(30^{\circ}-\chi)}, \tag{3}\]

where \(\chi\) is the chiral angle, \(\gamma^{\prime}_{A}\) and \(\gamma^{\prime}_{Z}\) are the edge energies of armchair and zigzag, respectively. \(\gamma^{\prime}_{A}\) and \(\gamma^{\prime}_{Z}\) were directly computed using the GGA as 1.03 eV/A and 1.19 eV/A, respectively.

On pristine graphene edges, the energies computed with the NNP showed higher stability of large chiral angles, which was in agreement with the GGA. On the other hand, those computed with the Tersoff-type potential indicated higher stability of small chiral angles, which was considered to result in the zigzag preferential growth by previous MD simulations [20].
### Thermodynamic properties of Fe bulk and nanoparticles

Influence of catalyst phase on SWCNT growth has been reported in experimental studies [45, 46]. In order to accurately simulate SWCNT growth, the thermodynamic behavior of catalyst nanoparticles must be well considered. In this section, the thermodynamic behavior of Fe bulk and nanoparticles is discussed. First, we estimated the melting point of bulk bcc Fe by MD simulation of solid-liquid biphasic system. For MD simulation, we obtained the initial structure of bulk bcc Fe with solid-liquid interface by preparing a supercell of \(28.52\times 28.52\times 74.17\) A, optimizing the lattice constant, and melting one of the halves divided along bcc (100) plane at 5000 K under NVT ensemble (Figure 3a). Starting from the initial structure, we carried out MD simulation for 200 ps under NPT ensemble at 1 atm. The final structures are shown in Figure 3b-d and the potential energies of each structure were well converged as shown in Figure 3e. As the temperature at which the solid and liquid phases reach equilibrium, the melting point was estimated to be 1887.5\(\pm\)12.5 K, which is in agreement with the experimental value 1811 K [47]. In addition, latent heat of melting was estimated as

Figure 2: Edge energy of graphene as a function of chiral angle. The values were directly computed (dots) and obtained from Eq. 3 (lines).
0.15 eV/atom from the difference of the averaged potential energies between solid and liquid, which is also in agreement with the experimental value 0.143 eV/atom [48] (Figure 3f).

Next, we verified the thermodynamic behavior of Fe nanoparticles by monitoring the Lindemann index representing the root-mean-square relative bond-length fluctuation

\[\delta=\frac{2}{N(N-1)}\sum_{i<j}\frac{\sqrt{\langle r_{ij}^{2}\rangle-\langle r _{ij}\rangle^{2}}}{\langle r_{ij}\rangle}, \tag{4}\]

where \(r_{ij}\) is the distance between atom \(i\) and \(j\), N is the number of particles and the bracket is the ensemble average calculated over MD simulation at constant temperature \(T\). The melting point was determined as the temperature at which the Lindemann index changed discontinuously. In previous studies, a threshold value in the range of 0.08-0.15 was used as the melting point criteria [49, 50, 51, 52, 53, 54]. However, in the case of small nanoparticles, the solid and liquid phases coexist over a wide temperature range and the criteria for the melting point has not been well established. In this study, we fitted the Lindemann index by a hyperbolic function expressed as

\[f(x)=\frac{(ax+b)\exp{(\frac{x-v}{s})}+(cx+d)\exp{(-\frac{x-v}{s})}}{\exp{( \frac{x-v}{s})}+\exp{(-\frac{x-v}{s})}}+w, \tag{5}\]

where a, b, c, d, e, s, v, and w are fitting parameters, and took an inflection point at \((v,w)\) as the melting point criteria (Figure S3).

For MD simulation to obtain the Lindemann index, the initial structures were prepared by melting Fe nanoparticles at 2000 K and then annealing them to 0 K at 1 K/ps. The annealed structures of Fe\({}_{200}\), Fe\({}_{600}\) and Fe\({}_{2000}\) are shown in Figure 4a-c. The stable phase of Fe nanoparticles changed depending on the size. For large nanoparticles (N\(>\)700), bcc structure was the stable phase as well as bulk Fe. However, for medium size nanoparticles (300\(<\)N\(<\)700), icosahedral fcc (Ih) structure was stable. For smaller nanoparticles (N\(<\)300), any regular structures were not observed. This phase transition from bcc to Ih by size, which was considered to be due to surface effect, is consistent with high-resolution transmission electron microscopy (HRTEM) observation of Ih Fe nanoparticles [55]. Starting from the annealed structure, MD simulations were performed for 2 ns at constant temperature \(T\) and then the Lindemann index was obtained as the average over the latter half 1 ns of MD simulations. According to Gibbs-Thomson equation, the decrease in the melting point is inversely proportional
to the diameter. As shown in Figure 4b, for large nanoparticles of bcc Fe, the melting point followed the equation. However, the melting point discontinuously changed by 300 K between Fe\({}_{600}\) and Fe\({}_{800}\). This discontinuity of the melting point corresponded to the phase transition from bcc to Ih.

The melting points of Fe\({}_{120}\)C\({}_{n}\) were also calculated as shown in Figure 4c. The melting point decreased from 814 K to 740 K with \(\sim\)2.5 % carbon content, while it increased with more carbon content. The observed V-shape diagram is consistent with the phase diagram of bulk Fe-C [56] and previous theoretical studies [49, 57].
Figure 3: Melting point measurement of solid–liquid biphasic system. (a) Initial structure. (b-d) The final structure after MD calculation for 200 ps by NPT ensemble under 1 atm at 1875 K, 1877.5 K and 1900 K. (e) Time evolution of potential energy at each temperature. (f) Potential energy distribution and its average per position of (c). The latent heat of melting \(\Delta\)H was estimated as 0.15 eV/atom from the difference of the potential energies between solid and liquid phase.
Figure 4: Thermodynamic behaviors of Fe nanoparticles. (a-c) Snapshots of cross sectional views of annealed structures of Fe\({}_{200}\), Fe\({}_{600}\) and Fe\({}_{2000}\) nanoparticles. (d) Melting point of Fe nanoparticles as a function of inverse diameter. (e) Melting points of Fe\({}_{120}\)C\({}_{n}\) nanoparticles obtained by adding carbon atoms up to 17.5 %. The blue and red balls denote atoms in BCC and FCC structure, respectively. The white balls denote atoms out of BCC and FCC structures. The structures were categorized by the common neighbor analysis [58].
## 4 MD simulation of SWCNT growth

We performed MD simulations of SWCNT growth on Fe\({}_{55}\), Fe\({}_{80}\), Fe\({}_{100}\) and Fe\({}_{120}\) nanoparticles at \(T\)=1100-1500 K in 100 K increments to investigate the growth trends by the temperature (Figure S4). Figure 5 shows the typical growth on Fe\({}_{120}\) nanoparticles at 1100 K and 1500 K. At low temperatures (-1300 K), carbon atom mobility was low and isolated hexagons nucleated simultaneously at several sites on the surface. Subsequently, the multiple graphitic islands were formed from the hexagons and the islands covered the surface of the nanoparticles and deactivated them. At high temperatures (1400-1500 K), a single hexagon nucleated, one large island was formed by the addition of carbon atoms, and the cap formation was observed.

We observed defect-free SWCNT growth at high temperatures, which was considered to be suitable for

Figure 5: Snapshots of the initial stages of SWCNT growth at (a) 1100 K and (b) 1500 K. The blue and green atoms denote iron and carbon, respectively.
SWCNT growth. Figure 6 shows the SWCNT growth process at 1500 K. At this temperature, SWCNTs grew through the typical growth stages; carbon saturation, cap formation and side wall elongation, as reported in computational [20, 28] and experimental [59] studies.

In the saturation stage, supplied carbon atoms in the gas phase were adsorbed on the catalyst surface and preferentially occupied the subsurface sites (Figure 6a). When the subsurface sites were fully occupied, carbon atoms started to form small molecules such as dimer and trimer on the surface (Figure 6b). When a sufficient number of carbon atoms gathered on the surface and the catalyst became supersaturated, the small molecules started bonding with each other and formed short chains. The first hexagon was formed by cyclization of the chain. If the chain acting as a precursor was larger than the hexamer, the excess carbon atoms were detached from the newly formed hexagon (Figure S5).

In the cap formation stage, additional hexagons were formed next to the first hexagon and a small graphitic island was formed (Figure 6d). The island was enlarged as carbon atoms were supplied and lifted off as a pentagon was formed (Figure 6e-f). The sixth pentagon formation finished the cap formation (Figure 6g). The cap pentagons were formed by the addition of a carbon atom at an armchair-like site and left inside the island by the formation of additional hexagons surrounding them (Figure S6). The formed cap was perfect, which had exactly six pentagons and no heptagons.

After the completion of the cap, the side wall of the SWCNT started to elongate. In this stage, defects such as pentagons and heptagons were formed on the open tube edge as well as hexagons, but healed to hexagons by bond recombination, which realized defect-free growth. Figure 6i,j compares the time evolution of the numbers of polygons: hexagons, pentagons and heptagons. The number of hexagons increased at constant rate after the first formation. We should note that in the case of the (8,7) SWCNT, the hexagon formation rate was 0.61 ring/ns and the growth rate in length was \(6\times 10^{5}\)\(\mu\)m/min, which was \(10^{2}\sim 10^{5}\) times higher than the experimental growth rates of SWCNTs. The number of pentagons increased to six during the cap formation and repeatedly increased and decreased but never fell below six during the sidewall elongation, which indicated that the pentagons in the cap were maintained. The number of heptagons was almost always zero, and even if it reached one, it immediately returned to zero.
Figure 7 shows the detailed mechanism of defect healing on the SWCNT edge. A pentagon was formed on the edge and healed into hexagons by the addition of monomer and C-C bond reformation (Figure 7a). Heptagons were also formed but less frequently than pentagons, and quickly healed to hexagons by rejection of monomer

Figure 6: (8,7) SWCNT growth on Fe\({}_{120}\) at 1500 K. (a-h) Snapshots of MD simulation. (i) Time evolution of the numbers of hexagons, pentagons, and heptagons. (j) An enlarged view of Figure (i). The red-colored polygons denote the pentagons.
(Figure 7b). Such defect healing resulted in defect-free SWCNT growth. In this study, we observed (8,7), (9,9) and (10,7) defect-free SWCNTs, which contained exactly six pentagons in the cap and no other defects (Figure S7).

The edge shape of the SWCNTs changed dynamically during the growth depending on the local structure such as armchair, zigzag and Klein edge. On armchair edges, pentagons were easily healed to hexagons by the addition of monomer because such pentagons had dangling bonds and highly reactive, as shown in Figure 7a. However, on zigzag edges, pentagons without dangling bonds were often formed at edge vacancies (Figure 8b). Such pentagons were less reactive, resulting in less susceptibility to monomer addition and chances of being left unhealed inside the sidewalls. They were occasionally left inside the sidewalls without being healed. Heptagons were often formed below the pentagons to compensate for the curvature change by the pentagon formation, resulting in pairs of adjacent pentagons and heptagons (5-7 defects). 5-7 defects caused chirality changes as reported in previous experimental studies [60, 61, 25]. An SWCNT with 5-7 defect mediated chirality change is shown in Figure 8a. The SWCNT grew up to a certain length with the initial chirality (11,7) determined by the cap and changed to (10,7) by a 5-7 defect. We also observed the chirality change from (10,3) to (8,3) by two 5-7 defect.

Figure 7: Edge healing processes on the developed view of the (8,7) SWCNT. (a) Reformation of a pentagon into a hexagon on the (8,7) SWCNT. (b) Reformation of a heptagon into a hexagon on the (8,7) SWCNT. The red and blue polygons denote the pentagon and heptagon, respectively.
defects (Figure S8).

On zigzag edges, adatoms often formed a protruding pentagon, which was thermodynamically unfavorable. The adatoms frequently jumped to adjacent zigzag sites with repeated bonding and dissociation like a 'trapeze'. And when reaching an armchair site, it formed hexagons and stabilized (Figure 9). Such adatom diffusion healed the vacancies, preventing topological defect formation (Figure 10).

In this study, unhealed vacancies occasionally resulted in topological defects as shown in Figure 8, because the growth rate in the simulations was much higher than that in the experiments, and vacancy healing did not occur in time for SWCNT growth. In the experiments, however, SWCNTs are considered to grow while maintaining smooth edges through adatom diffusion, resulting in a low probability of defect formation.

Figure 8: Snapshots of an SWCNT with 5-7 defect mediated chirality change. (a) Snapshot and developed view of a SWCNT with chirality change from (11,7) to (10,7) after 280 ns of carbon supply. (b) 5-7 defect formation process. The pink and orange areas are (11,7) and (10,7) sidewalls, respectively. The light blue areas can be defined as both (11,7) and (10,7) since there are two possible translation vectors. The chirality of the white areas can not be assigned. The red and blue polygons denote the pentagons and heptagons, respectively.
Figure 10: Edge vacancy healing by adatom diffusion on the (8,7) SWCNT.

Figure 9: Adatom diffusion on the zigzag edge of the (10,7) SWCNT.
## Conclusions

We demonstrated defect-free SWCNT growth by MD simulation based on the developed NNP. During the growth process, it was observed that defects were frequently formed on the open edge and healed into hexagons by metal mediated C-C bond recombination, resulting in perfect SWCNTs with (8,7), (9,9) and (10,7) chiralities. Furthermore, it was shown how 5-7 defects were formed at edge vacancies, which remained unhealed in the sidewall and caused the chirality change. Also, it was observed that adatom diffusion contributed to keeping smooth edge without vacancies, which indicated the origin of why a pre-determined chirality was sustained in the growth of SWCNTs, a key mechanism for the chirality selectivity.
