Free energy surface of initial cap formation in carbon nanotube growth

Initial cap formation is an important process of carbon nanotube (CNT) growth where a hexagonal carbon network is lifted off from the catalyst surface. In this study, free energy surface (FES) of initial cap formation in the CNT growth is investigated by metadynamics simulation. A two-dimensional collective variable (CV) space is newly developed to examine the complicated formation process of the cap structure, which consists of the formation of a hexagonal carbon network and lift-off of the network from the catalyst surface. States before and after the lift-off of the carbon network are clearly distinguished in the two-dimensional FES. Therefore, free energy difference before and after the lift-off can be directly derived from the two-dimensional FES. It was revealed that the cap structure is stable at a high temperature due to the entropy effect, while the carbon network covering the catalyst surface is energetically stable. The new insight in this study is achieved owing to metadynamics simulation in conjunction with a newly developed two-dimensional CV space since it is impossible to explore FES for such complicated processes in the framework of conventional molecular dynamics simulation.


Introduction
Carbon nanotubes (CNTs) 1 can potentially be used in a variety of applications. 2 It is essential to control the diameter and chirality of CNTs during the synthesis for better application. The catalytic chemical vapor deposition (CCVD) method 2 is widely used to synthesize CNTs, in which CNTs grow from dispersed nanoparticles. 3 Diameter and chirality of CNTs are strongly correlated with the properties of catalytic nanoparticles at the initial stage of the synthesis process. Therefore, it is important to understand the role of catalytic nanoparticles at the initial stage of CNT growth. Molecular dynamics (MD) simulations [4][5][6][7][8][9] have contributed to the understanding of the metal-catalyzed growth mechanism of CNTs from an atomistic viewpoint. In conjunction with in situ transmission electron microscopy (TEM) observations, 10,11 the initial growth process of CNTs by CCVD is considered as follows: 12 aggregation of carbon atoms in/on the catalytic nanoparticle, formation of the hexagonal carbon network on the surface of catalytic nanoparticles and the li-off of the carbon network to form a cap structure.
In general, an MD simulation represents one specic trajectory starting from an initial conguration and does not cover all the ensemble space. It is theoretically possible to cover all the free energy surface of target phenomena if the MD simulation is performed for a sufficiently long time. However, it is not realistic since the time scale of MD simulation is limited to the nanosecond order at most. Therefore, most of the MD simulations of the CNT growth were limited as an enormous part of possible trajectories are cut out 3,4 and it is not clear that the discussion based on one-time MD simulation covers the universal feature of CNT growth. To overcome this problem, metadynamics 13 is proposed in which the efficiency of the sampling of target phenomena is improved by setting collective variables (CVs) and applying a bias to the CVs. Target phenomena can be efficiently sampled in conjunction with appropriate CVs. Up to now, numerous CVs have been developed for various phenomena such as bond dissociation reaction, 14 ligands docking on receptors, 15 atomic diffusion 16 and solid-liquid phase transition. 17 For example, a bond length between atoms is used as a CV for describing bond dissociation and docking reactions. In these simulations, a bond for the target reaction is pre-dened and its length is simply used as CV. However, it is not straightforward to set appropriate CVs for complex and multi-step phenomena. For example, the initial formation process of CNTs consists of the formation of a carbon dimer, junction and hexagonal network, and subsequent li-off of the network to form a cap structure. We previously developed a CV based on the coordination number of carbon atoms to discuss a part of the initial formation process of CNTs and we sampled the potential energy surface of carbon segregation from a nickel nanoparticle by metadynamics simulation. 18 Although the states of isolated atoms, dimer, junction, and triple junction was properly distinguished, it was not successful to cover the subsequent formation process of a hexagonal carbon network and li-off of the network to form a cap structure using one unied CV. To this end, a twodimensional CV space representing the formation of a hexagonal carbon network and li-off of the network was newly developed in this study. Using the newly developed CVs, the free energy surface (FES) of the cap formation process of CNTs is investigated by metadynamics simulation.

Free energy calculation with metadynamics
First, a method to calculate the FES by MD simulation is described. The FES for CV q is obtained by sampling the appearance probability P(q) of q as Here, b is 1/k B T, where k B is the Boltzmann constant and T is the temperature. In an ergodic system, P(q) can be obtained by simply recording a histogram N as function of q as eqn (2) with sufficiently long MD calculation. 19 The total free energy of state A, G A , is calculated as: However, it is not straightforward to achieve sufficient sampling within the time scale of normal MD simulation.
In this study, metadynamics was employed to improve the efficiency of sampling. In metadynamics, a bias potential in the form of Gaussian functions with width s and height W dened as: is added to the CV at every ks, where k is an integer and s is the interval. The height W is modied based on the well-tempered metadynamics 20 scheme as Here, W 0 is the initial height and DT is calculated from a bias factor g as P(q) can be obtained by creating a histogram weighed by exp(bV bias (q) À bc(t)) during the calculation. 21 Here, c(t) is dened as 2.2. Two-dimensional collective variable space for the formation of a carbon nanotube cap As described in the introduction, it is most important to dene appropriate CVs for metadynamics simulation. In this study, a two-dimensional CV space is newly developed to represent the formation process of the cap structure of CNT. This consists of the formation of a hexagonal carbon network on the surface of a catalytic nanoparticle and the li-off of the network from the surface. It is essential to distinguish two states before and aer the li-off of the carbon network from the catalyst surface. The hexagonal carbon network covers the catalyst surface in the former (hereinaer called covering-state) and the carbon network is lied off from the catalyst surface forming a cap structure in the latter (hereinaer called cap-state). The rst CV q 1 was dened to describe the network formation. The denition consists of the following three steps (see Fig. 1(a)). (i) Coordination number c i of carbon atom i, for other carbon atoms j is dened using switching function s as Here, r ij is the distance between atoms i and j. Cutoff parameters d 0 and r 0 represent offset distance and width at half maximum, respectively. Exponents n and m represent the degree of decay of switching function at the cutoff distance. Parameters in eqn (10) are employed from our previous work. 15 Coordination numbers of 0, 1, 2 and 3 represent isolated atom, atom in dimer, atom in chain and atom in the center of triple junction, respectively. (ii) The local average of the coordination number is calculated since the coordination number is not sufficient to distinguish the atom in the center of a triple junction from one in the connected hexagonal network. The local average l i of the coordination number around atom i is calculated using the weighting factor w ij as The local average of the coordination number for the atom in the complete hexagonal network is assigned as 3, whereas that of the imperfect hexagonal network is assigned as less than 3, as shown in Fig. 1(a). Therefore, the local average of the coordination number can distinguish these structures correctly. (iii) The largest value of local averages of coordination number assigned to each atom is selected as the rst CV q 1 . Since the CV must be a continuous value to add bias in metadynamics, eqn (13) is employed for the simulation.
a ¼ 0.01 is used in this study. The second CV q 2 is dened to represent the li-off state of the carbon network structure from the nanoparticles. The number of atoms located above a certain distance from the center of gravity of the catalyst counted using switching function is employed as Here, r i is the distance between atom i and the center of the gravity of the catalyst. The second CV q 2 was used only to distinguish states whether the hexagonal carbon network is lied-off or not. Note that q 2 is not used to accelerate the reaction itself.

Considered systems and conditions
Metadynamics simulation of the cap formation process of CNTs is performed using above dened two CVs. The simulation methodology basically follows our previous work 18 except for the denition of CVs and the calculation system. A Ni 55 C 110 nanoparticle placed at the center of a 35.24 Â 35.24 Â 35.24Å 3 cubic cell is employed as the initial conguration. The nanoparticle is prepared by annealing Wulff-shaped Ni 55 and randomly placed 110 carbon atoms at 2000 K for 100 ps rst and then at 1500 K for 200 ps. Ni was chosen as an element of the catalytic metals since it is one of the most popular catalysts used in the CNT synthesis and is a target of numerous previous computational works. 3,4,8,9,18 The number of Ni atoms is set as 55 to make the particle diameter around 1 nm, which is consistent with the diameter of a small single-walled CNT. It was conrmed from our previous simulation that the number of carbon atoms must be at least larger than that of Ni atom for spontaneous precipitation. 22 Therefore, the number of carbon atoms was set as 110, which is large enough to form the initial cap structure. Repulsive Gaussians of height W 0 ¼ 0.043 eV and width s ¼ 0.05 were employed as biasing parameters in eqn (5). Bias potential is added at every 200 steps. A bias factor of 200 was used for the well-tempered metadynamics scheme. A repulsive wall potential (spring constant of 830 eV) was set to prevent carbon atoms from getting too close to each other and it also acts on atoms with the largest value of coordination number c i when it exceeds 4.2. A ReaxFF reactive force eld 23 was employed to describe the interaction between nickel and carbon atoms, which has already been successfully applied to simulations of CNT nucleation and growth. 8,9 The NVT (the number of atoms, cell size and temperature constant) ensemble was sampled using a Nosé-Hoover thermostat. 24,25 Simulations were carried out for seven different temperatures, viz. 1000, 1250, 1500, 1750, 2000, 2250 and 2500 K. 3.0 Â 10 7 simulation steps were carried out for all cases with a time step of 0.25 fs. All calculations were performed using a Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) 26 with PLUMED plugin. 27   Nanoscale Advances carbon network during the metadynamics simulation. In the normal MD simulation, reactions usually proceed in one direction and the reverse reaction is rarely observed. That is, the hexagonal carbon network does not disappear once it forms 4,8 in the normal MD simulation. Therefore, this oscillation veries that the metadynamics using q 1 successfully accelerated the sampling of possible states including the hexagonal carbon network in the metadynamics simulation. Note that the value of q 1 can exceed three since carbon atoms further than the bond length can have a nite contribution to the coordination number due to the denition of q 1 as having a continuous value. Fig. 3(a) shows the FES along q 1 obtained from the metadynamics simulation at 2000 K. The FES indicates that the states between q 1 ¼ 3.6 and q 1 ¼ 3.8 are the most stable. The insets in Fig. 3 show representative structures for the covering-and capstates. The values of q 1 for both structures are close to 3.8 and therefore it is not able to distinguish these two structures by q 1 only. Therefore, FES is expanded to a two-dimensional space using q 2 . Fig. 3(b) shows the two-dimensional FES along q 1 and q 2 obtained from the metadynamics simulation at 2000 K. Note that the bias potential is added to q 1 only during the metadynamics simulation and q 2 is just recorded during the simulation. Therefore, the simulation itself is the same even if two CVs are recorded. The states successfully separated into different spaces in the two-dimensional FES. That is, covering-and capstates around q 1 ¼ 3.8 spread out horizontally in the twodimensional FES. Fig. 4 shows two-dimensional FES obtained from metadynamics simulations at other temperature conditions. The insets in Fig. 4 show representative structures for the covering-state and cap-state. Both covering-and cap-states are observed in the separated region of two-dimensional FES at all temperature conditions. As the temperature increases, the region of stable state in FES distributes upward at around q 1 ¼ 4. The value of q 2 does not exceed 20 at 1000 K, whereas the values of q 2 exceed 40 at 2250 K and 60 at 2500 K. Moreover, the area of cap-state in FES becomes large with an increase in temperature. This means that the large cap structure is more stable at high temperatures, which agrees with the direct observation of snapshots from the simulation. For example, the cap structure in the snapshot at 1000 K (q 2 ¼ 15.4) has a diameter of approximately 4Å, while that at 2500 K (q 2 ¼ 43.7) corresponds to a diameter of approximately 10Å. Regarding the growth mode, perpendicular growth (i.e., the line contact for the catalyst-tube contact) appears at the Ni catalyst with a high carbon concentration in our simulation. It agrees with the previous report 28 that a high carbon concentration inside the metal particle favors perpendicular growth mode, whereas low carbon fractions in the catalyst lead to tangential growth mode (i.e., the surface contact). Note that the diameter of the cap structure is also affected by the diameter of the catalyst, which is however out of the scope of this study.

Free energy difference between covering-state and capstate
Next, the free energy difference between the covering-and capstates is investigated. The free energy of each state can be derived from eqn (3) once the region of each state is dened in FES. Here, the region of q 1 $ 3.4 and q 2 < 7.0 in FES is regarded  as the covering-state and that of q 1 $ 3.7 and q 2 $ 7.0 is regarded as the cap-state. The free energy difference DG is calculated by the difference between the free energy of the covering-and cap-states as DG ¼ G cap À G covering . Fig. 5 shows the free energy difference between the covering-and cap-states as a function of temperature. A positive value of DG indicates the free energy of the covering-state is lower (i.e., more stable) than that of cap-state and vice versa. As shown in the gure, the covering-state is more stable than the cap-state at a low temperature and DG becomes small with an increase in temperature. The cap-state becomes more stable than the covering-state at 2500 K. The line tted to the plots corresponds to the equation DG ¼ DH À TDS according to the denition in thermodynamics, where DH and DS represent enthalpy difference and entropy difference between covering-and cap-states, respectively. From the vertical intercept and the slope of the tted line, DH and DS are estimated to be 0.74 eV and 0.33 Â 10 À3 eV/T, respectively. This suggests that the covering-state is enthalpically (i.e., energetically) stable due to the positive value of DH. However, the cap-state can be stable at high temperatures since entropy effect becomes dominant via the term ÀTDS. This result makes sense since the hexagonal carbon network covering the catalyst surface has more catalyst-carbon bonds than the network lied-off from the surface. The catalyst-carbon bonds make the covering-state more stable energetically. In contrast, the cap structure is entropically stable since there is no limit to the number of atoms consisting of the cap structure, whereas there is a geometric constraint in the  arrangement of the carbon network covering the nanoparticle surface. In fact, the region of the stable cap state is spread out along the q 2 axis in FES at 2500 K, which indicates that cap structures consisting of numerous atoms are sampled in the metadynamics simulation. This result provides a reasonable answer to the long-standing issue that the li-off of the cap structure is hard to occur at low temperatures 5 from a thermodynamic viewpoint.

Conclusion
Initial cap formation process of the CNT growth is investigated by metadynamics simulation on a two-dimensional CV space. Newly developed two-dimensional CV successfully distinguished states of the hexagonal carbon network covering on the catalyst surface and that lied-off from the catalyst surface in the FES. The free energy difference between the covering-and cap-states with respect to temperature suggests that the capstate is stable at high temperatures due to the entropy effects, while the covering-state is energetically stable. These ndings suggest appropriate directions for efficient CNT synthesis and the formation of the CNT cap is suitable at a higher temperature although the lower temperature is in general preferred as it prevents the nanoparticles from aggregating. Moreover, our result provides a reasonable answer to the long-standing issue that the li-off of the cap structure is hard to occur at low temperatures. In the practical synthesis process, additional carbons are supplied throughout the reaction process although the number of carbon atoms is xed in this study. It is expected that the degrees of freedom of carbon atoms for the cap structure will increase as the number of carbon atoms increases. The cap structure is expected to be more stable from an entropy viewpoint with increasing degree of freedom even at lower temperatures. The effect of the number of carbon atoms as well as the effect of the size of catalytic nanoparticles will be examined in the next study.

Data availability
The data that support the ndings of this study are available from the corresponding author upon reasonable request.

Conflicts of interest
There are no conicts to declare.