Molecular dynamics simulation of graphene sinking during chemical vapor deposition growth on semi-molten Cu substrate

Copper foil is the most promising catalyst for the synthesis of large-area, high-quality monolayer graphene. Experimentally, it has been found that the Cu substrate is semi-molten at graphene growth temperatures. In this study, based on a self-developed C–Cu empirical potential and density functional theory (DFT) methods, we performed systematic molecular dynamics simulations to explore the stability of graphene nanostructures, i.e., carbon nanoclusters and graphene nanoribbons, on semi-molten Cu substrates. Many atomic details observed in the classical MD simulations agree well with those seen in DFT-MD simulations, confirming the high accuracy of the C–Cu potential. Depending on the size of the graphene island, two different sunken-modes are observed: (i) graphene island sinks into the first layer of the metal substrate and (ii) many metal atoms surround the graphene island. Further study reveals that the sinking graphene leads to the unidirectional alignment and seamless stitching of the graphene islands, which explains the growth of large single-crystal graphene on Cu foil. This study deepens our physical insights into the CVD growth of graphene on semi-molten Cu substrate with multiple experimental mysteries well explained and provides theoretic references for the controlled synthesis of large-area single-crystalline monolayer graphene.


INTRODUCTION
Following its first preparation by exfoliation from graphite in 2004, graphene has attracted much attention both for fundamental studies and technological applications, owing to its unique morphology and excellent mechanical, electronic and optical properties. It is expected that wafer-scale single-crystalline graphene films will function as ideal platforms for the development of future high-performance graphene-based devices 1 . Consequently, great efforts have been made to synthesize large single-crystalline graphene (SCG) monolayers 2-4 on transition metal catalyst through chemical vapor deposition (CVD) growth -the most promising, cheap and readily accessible synthetic approach so far reported for graphene 5 . The graphene islands, formed during the initial nucleation process [6][7][8][9][10][11] , are crucial in the mechanism of graphene growth and determine the quality of the grown graphene by functioning either as nucleation seeds 6,12,13 or as building blocks for the coalescence of graphene nano-islands 9 . In previous studies, large SCG was realized by the seamless stitching of highly oriented islands on single-crystal substrate, such as Cu(111) 14 and Ge (110) 15 . Most recently, ultra-fast growth of meter-sized SCG has been achieved by producing super-large Cu (111) foil 16 . As an epitaxial growth, the unidirectionally aligned graphene islands on the Cu (111) surface are guided by the lattice symmetry of the single-crystal substrate 17 .
Before graphene islands form, various intermediates such as carbon dimers, atomic carbon nanoarches or networks of carbon clusters, have been observed both in early experiments 18,19 and theoretic calculations [20][21][22] . For instance, based on densityfunctional molecular dynamics simulations, the very intricate pathways of reaction for metal organic chemical vapor deposition (MOCVD) of AlN on graphene can be identified at the atomic level with rich chemistry uncovered 23 . Unexpected defects such as pentagons usually appear during the nucleation stage and degrade the final quality of CVD-grown graphene. To circumvent the formation of these detrimental defects created during the initial stages of growth, defect-free polycyclic aromatic hydrocarbons (PAHs consisting of six-membered carbon rings) have been proposed as ideal precursors for the fabrication of graphene with low defect concentration, on Au (111) 24 , Pt (111) 25 and Cu (111) 26 . Among these different carbonaceous motifs, the coronene-like cluster, C 24 , is widely regarded as the ideal carbon precursor for the growth of high-quality graphene by selfassembly on a Cu (111) surface 27,28 due to its six-fold rotation symmetry and dome-shaped structure, a result of the interaction of the peripheral atoms of the carbon cluster with the substrate metal [29][30][31] . Although it is widely believed that these active carbon species readily move on the metal surface, coalescing into larger islands, and eventually forming the final graphene domain 27,32 , the motion and coalescence dynamics at the atomic level of these carbon islands still remain elusive. For example, based on static DFT calculations, a new growth mode called "embedded" growth, has recently been proposed. In the embedded mode, it was found to be energetically favorable for graphene islands to "sink" into "soft" metals like Cu (111), during the growth of graphene 33 .
Different from that on nickel surface, the growth of graphene on copper surface is a surface-mediated process due to the very low solubility of carbon in bulk Cu [2][3][4]34,35 . Typical experimental temperatures during growth are~1273 K (1000°C), which is very close to the melting point of bulk Cu at~1358 K (1085°C) 4,14,16,36 . Under this temperature, the surface of the Cu substrate is observed to be highly mobile 37 , which presents the characteristics of melting. Thus the Cu surface is thought to be nearly molten or liquid and can never be considered as a strict crystalline lattice as before 38 . How does the graphene island sink into such a "semimolten" Cu surface? Why are the graphene islands aligned on the surface molten Cu surface? Will the graphene sinking create detrimental effects on the stitching mechanism of the graphene grains or affect the final quality of larger merged graphene film? In this article, based on a self-developed C-Cu classical potential and density functional theory, MD simulations were performed to study the sinking of the graphene on a semi-molten Cu substrate. Two different sunken-modes, which depend on graphene island size, were revealed. Further, the embedded graphene tends to be unidirectionally aligned and the seamless stitching of the adjacent graphene grains is preferred.

RESULTS AND DISCUSSION
Sinking of carbon nanocluster Figure 1a-f show the cross-sections of a C 24 cluster on the Cu (111) substrate after 120 ps MD simulations at 1200, 1250, 1300, 1350, 1400, and 1450 K, respectively. The 3D views of the corresponding configurations are shown in Supplementary  Fig. 2. In our MD simulations, only the top layer Cu atoms melt while other layers of the substrate remains in a crystalline shape even at 1450 K owning to the approach of semi-infinite surface adopted in our model. Thus the substrate is in a state between liquid and solid. We name the substrate a "semi-molten" substrate. As shown in Fig. 1h, during the MD simulation, the Lindemann index 39 of the substrate oscillates between those of liquid and solid. As shown in Fig. 1a, b, the upmost surface of Cu substrate is stable and retains its original crystallinity at lower temperatures (T Previous studies have reported that during graphene growth, graphene edges bound to a metal step are energetically preferred 6,8,40 . It has also been observed that graphene ribbons tended to stand upright on a flat metal terrace 41 while a graphene island like C 24 prefers to form dome-like structure with its edge bent towards the terrace 20 . In our MD simulation, the C 24 cluster also presents a dome-like geometry (Fig. 1a), similar to that reported in the previous studies 20 . The height of the dome-like C 24 on Cu surface is~0.6 Å, which is~50% of that of C 24 on a Ni (111) surface 42 . At 1300 K, the top layer of Cu atoms start to melt slightly with a small proportion of Cu atoms diffusing out of the surface to attach to the edge of C 24 (Fig. 1c). These diffused Cu atoms can be viewed as a newly-formed atomic layer passivating the edge of the C 24 cluster. Therefore, different from the dome-like flake on the Cu surface shown in Fig. 1a, C 24 can form metal-carbon σ bonds with surrounding Cu atoms and its shape becomes flat (Fig.  1c). When the temperature approaches 1350 K, the substrate becomes highly disordered and the C 24 starts to sink into the substrate (Fig. 1d). The higher the temperature, the deeper the C 24 cluster sinks (Fig. 1e, f). Besides, as it sinks deeper, the domelike structure of C 24 becomes increasingly flatter (Fig. 1g). To quantify the degree of sinking of the C 24 cluster, the average heights (z coordinates relative to z = 0) along the z direction (see the axes in Fig. 2d) of the mass centers of C 24 (h-C 24 ) and the surrounding Cu atoms (h-SCu) during the final 20 ps are compared in Fig. 1i. With the increase of the temperature, the h-C 24 decreases obviously whereas h-SCu remains almost unchanged with only a slight fluctuation at~3.5 Å. The height difference Δh, where Δh ¼ h À C 24 À h À SCu, therefore drops to below 0 Å at temperatures higher than 1350 K, indicating the complete sinking of the C 24 cluster.
In order to investigate the detailed sinking process, the smoothed instantaneous heights of C 24 (h-C 24 ) and surrounding Cu atoms (h-SCu) as a function of MD time are also calculated. As shown in Supplementary Fig. 1, at temperatures <1300 K, the h-C 24 and h-SCu are almost parallel, implying that both the sinking of C 24 and the surface melting cannot take place (Supplementary Figs. 1a, b). At 1300 K, both h-C 24 and h-SCu start to rise after 90 ps, which can be mainly attributed to change in the mass centers of C 24 and the surrounding Cu atoms (Supplementary Fig.  1c). In other words, the original dome-like shape of C 24 , stable at lower temperatures, begins to flatten at 1300 K due to edge passivation by several Cu atoms that diffused upwards. Therefore, the mass center of the C 24 cluster is raised by~0.3-0.5 Å, which is consistent with a previous density function theory (DFT) calculation on the height difference between C 24 and coronene on a metal substrate 20 . In addition, the mass center of the surrounding Cu atoms is also raised by~0.5-0.75 Å because of their upward diffusion. Up to this point, we can see that the C 24 floats on the Cu (111) surface without sinking. However, as shown in Fig. 2a, when the temperature is 1350 K, the h-C 24 begins to drop sharply bỹ 1 Å after~60 ps (red arrow). The height difference between C 24 and Cu surface is greatly reduced from~2.25 to~0.5 Å at 66.5 ps. After 110 ps, the height difference is further reduced to be less than 0.20 Å, indicating that the C 24 has already been fully embedded in the semi-molten Cu surface. The potential energy profile in Fig. 2g shows that the sinking process after 60 ps is energetically preferred because of the downward trend. The snapshots of the sinking process at five different MD times are shown in Fig. 2d. Clearly, at 54.35 ps, the C 24 is still on the Cu surface. At 66.5 ps, the C 24 sinks into the Cu. After 70 ps, the C 24 rises up at 87.2 ps and fall down at 89.6 ps, corresponding to the height fluctuation of the C 24 in Fig. 2a. At 112.35 ps, the C 24 is fully embedded in the Cu substrate with their height difference as small as 0.13 Å. With the increase of temperature, the onset time for C 24 sinking, denoted by red arrows, is advanced from 60 (T = 1350 K) ps to 40 ps (T = 1400 K) and then to 10 ps (T = 1450 K) (Fig. 2b, c). Moreover, the height difference between C 24 and surrounding Cu atoms are further reduced even with higher Cu atoms observed in Fig. 2b, c. In previous studies, Yuan et al. proposed that the sinking of graphene in CVD growth can be achieved either by removing metal atoms under a graphene island, or by the adsorption of fast moving metal atoms around a graphene island to form a new atomic layer 33 . Clearly, the former pathway is the dominant sunk mechanism at high temperatures ≥ 1350 K. To distinguish this sunken-mode with the later one, we define it as "sunken-mode I".
For comparison, Fig. 2f shows the C 24 embedded in the liquid substrate, 8-layers free Cu atoms, at 1350 K. Different with the semi-molten substrate, the crystallinity of the Cu layers completely disappears. As shown in Fig. 2e, the heights of the C 24 cluster drops from 10.5 to 8.75 Å very quickly in only 20 ps with the height difference reduced to~0.15 Å. After 67 ps, the h-SCu begins to exceed h-C 24 with the height gap near 0.5-1 Å. Compared with the semi-molten substrate, it is, therefore, easier for C 24 to sink and be embedded in the liquid surface. In order to verify the observations in the classical MD simulation shown above, we performed an extensive DFT-MD simulation to explore the behavior of a C 24 cluster on a Cu (111) substrate at 1350 K. Owning to the high computational cost of DFT-MD, we use 3-freelayers Cu (111) substrate to accelerate the surface melting, which can reduce the computation time effectively. As shown in Supplementary Fig. 7, at 0.23 ps, the surface has not melted yet and the dome-like C 24 appears with the dome height~0.9 Å ( Supplementary Fig. 7b), which is close to 0.6 Å in our classic MD simulations. After 2 ps, an obvious sinking of the C 24 cluster into the substrate can be observed ( Supplementary Fig. 7c-f) with the height difference~0.5 Å (Supplementary Fig. 7g). Thus, the results of the classic and DFT MD are qualitatively consistent.
After knowing the sunken-mode of C 24 at higher temperatures, we may ask a question as to whether a larger graphene island also shows the same sunken-mode or not. To answer this question, larger carbon clusters of C 54 and C 96 , including 19 and 37 hexagons were selected for further simulations. As shown in Fig.  3a, unlike C 24 , the C 54 did not sink into the substrate but just float on the substrate surface at 1350 K. Only some melted Cu atoms diffuse upwards and form metal step with C 54 edge partially passivated. For this reason, Fig. 3d shows that there is no obvious drop of C 54 height line and a height gap of~1.5 Å between C 54 and the underlying Cu surface always exists (Fig. 3d). Even at higher temperatures like 1400 and 1450 K (Fig. 3b, c), C 54 prefers to float on the Cu surface and the surface Cu atoms eventually melt and diffuse upwards to surround it. Although the height gap between C 54 and surrounding Cu atoms is greatly reduced in Fig. 3 (e, f), this behavior can be attributed only to the lifting-up of the surrounding Cu atoms rather than the sinking down of C 54 . In other words, the C 54 cluster cannot sink into the substrate as easily as C 24 . Thus, we define such a sunken-mode as "sunken-mode II", which is different form "sunken-mode I" aforementioned. The reason for this difference is that, due to the larger size of the C 54 cluster, more underlying Cu atoms need to out-diffuse to form a basin for the cluster to sink, which in turn, requires a greater energy cost. Besides, the C 54 cluster is energetically more stable on the Cu (111) surface owing to its lower edge atoms/inner atoms ratio. Thus, there is a smaller driving force for the sinking of the C 54 cluster when compared to C 24 . In this respect, the sunkenmode II of C 54 at high temperatures is achieved by the adsorption of fast moving metal atoms underlying a graphene island to form a new atomic layer. As for a much larger cluster like C 96 , it always floats on the Cu (111) surface for a wide range of temperatures from 1300 to1450 K (Supplementary Fig. 4). We can thus conclude that the larger carbon nanoclusters prefer the "sunken-mode II" with only the carbon edges embedded in the metal atoms diffused from the underlying semi-molten substrate.
Orientation of carbon nanocluster It is well known that the orientation of the graphene islands formed at the early stage is very crucial to the final formation of larger graphene, such as the morphology, grain boundary, moiré pattern and so on 43,44 . By controlling the orientation of graphene island as well as associated growth conditions, various shapes of graphene grains can be synthesized 12,[45][46][47] . To eliminate the anisotropic effect of the solid surface on the orientation of graphene domain, the liquid Cu surface is proposed to provide isotropic surface for the uniform formation of the graphene grains 35,48,49 .
Recent theoretic study shows that coalescence of randomly rotated graphene islands on the liquid Cu surface only has the 0°( single crystal) and 30°-twinned polycrystal as highly stable grain boundaries. In our simulation, the orientation of the C 24 is proved to depend strongly on surface states of the metal catalyst at different temperatures. At the lower temperature far below the melting point, the C 24 cluster only vibrates rigidly, and sporadically, it rotates by~20°and then returns to 0°very soon ( Supplementary Fig. 5). These rotations only result in many sharp peaks but with no obvious orientation change. Based on our previous DFT calculation 43 , the orientation of graphene island, like C 24 and C 54 , is mainly constrained in a specific alignment at lower temperatures by the interaction between the edge atoms and Cu surface with the rotation barrier as high as 0.96 and 5.18 eV, respectively. For this reason, the rotation freedoms of larger clusters of C 54 and C 96 are expected to be reduced compared with C 24 at the same temperatures due to the higher rotation barriers ( Supplementary Fig. 5). When the temperature rises to 1350 K or higher, the C 24 starts to rotate gradually on the semi-molten surface as a function of MD time at temperature of 1350 K with its orientation changed apparently ( Fig. 4e and Supplementary Fig.  5). It is noticed that the distribution of the rotation angle for C 24 also depends strongly on the surface state. On one hand, as shown in Fig. 4a, c, the rotation angles on the semi-molten surface at 1350 and 1400 K present narrow distributions only at 0°and 60°, respectively. Due to unique sunk-mode of the C 24 aforementioned, its orientation is constrained by the unmolten subsurface with six-fold lattice symmetry. On the other hand, there is a broader distribution of the rotation angles for the C 24 on the liquid Cu surface at the same temperatures due to the isotropic lattice of liquid Cu (Fig. 4b, d).
Stitching of graphene As mentioned before, the seamless stitching of the graphene during the merge of the graphene domains is a significant step for formation of large-area single crystal graphene [14][15][16]50 . When the graphene sinking occurs, will these separated graphene domains coalescence into a larger one without defects? And, what are the atomistic details? To answer these questions, we select periodic graphene nanoribbons with a narrow gap as stereotypical model for mimicking the stitching process. Firstly, three graphene nanoribbons on Cu (111) surface, containing 3, 5, and 8 rows of hexagonal rings along the width direction (GNR-3, GNR-5, and GNR-8), are designed to study their structures on the Cu substrate. Figure 5 shows their structures after MD simulation of 200 ps at temperatures of 1200, 1350, and 1400 K, respectively. Clearly, at the lower temperature 1200 K, the upmost surface has not started to melt and the edges of all GNRs bend toward the Cu surface with arched ribbons formed (Fig. 5a-c). When the temperature is increased to 1350 K, the higher temperature leads to the formation of a metal step (Fig. 5d-f) with one-atom layer at one side of the GNR edges. The arched ribbons are therefore flattened due to the edge passivation. With further increase of the temperature, more metal atoms in the substrate diffused out and form the metal steps for GNR-3 and GNR-5( Fig. 5g-i). Thus, the GNRs on the semi-molten Cu substrate can be viewed as embedded ribbon based on sunken-mode II. Again, the DFT-MD simulations are carried out to verify our classic MD results. As shown in Supplementary Fig. 8, the GNRs bend toward to the surface and form the arched ribbon at the early stage when the surface has not melted. After~6-7 ps, the arched GNRs are flattened with the metal step formed near the edges. Height of the ribbon edges and the metal steps during the final stage indicate that the ribbon edges are embedded in the metal steps, which are highly consistent with our classic MD results in Fig. 5. So far, a clear conclusion is that with increase of the size, the carbon nanoclusters and/or GNRs can be well flattened on the Cu surface by the edge passivation of the metal steps, particularly on the substrate with surface melting. The following question, however, is that will these metal steps affect the final formation of the high-quality graphene, particularly on the semi-molten substrate. Or say, will these metal steps remain there and create defects during the stitching of the GNRs or disappear eventually with perfect graphene left? To solve this problem, we performed a MD simulation by adding carbon atoms randomly into the gap between two edges of GNR-8 on the Cu substrate. Very surprisingly, as shown in Fig. 6a-j, the metal step (yellow atoms), which is clamped by the two GNR edges, disappears with the stitching of graphene (red carbons) and a near perfect graphene sheet is eventually formed with only one 5|7 defect left. The atomistic details show that the atoms of metal steps can not only diffuse into the substrate below but also facilitate the defect healing and the hexagon formation. As shown in white rectangles of Fig. 6c, d, with the assistance of the yellow atoms, the two carbon chains are connected and transformed to a hexagon successfully. Also, a 5|7 defect in Fig. 6e is healed to be two hexagons (Fig. 6f, g). Besides, a carbon vacancy can also be healed by a carbon atom newly added on the upmost Cu atom. As shown in Fig. 6h-j, with the sinking of upmost Cu atom (yellow) into the substrate beneath graphene, the adatom C fills the carbon vacancy perfectly. The stitching of graphene on Ni (111) substrate can also occur between two adjacent graphene ribbons as shown in Supplementary Fig. 9 although the atomic details are different from Cu. Because of the higher melting point of the nickel substrate as well as the large carbon solubility in Ni, our simulations show that a simulated graphene ribbon cannot sink as easily as on a Cu substrate and tends to maintain an arc-bridge shape on the flat Ni substrate. Further, compared with graphene stitching on a Cu substrate at the same temperatures, more defects in the stitching area are observed. Therefore, the sinking of the graphene on a Cu substrate promises a better method for graphene stitching.
Here, we may propose a new explanation to the mysteries that why semi-molten substrate is suitable for the growth of large SCG graphene from another point of view by fully addressing the significance of the unique sinking kinetics discovered in our results. In details, the orientations of the small carbon clusters at the early stage of graphene growth can be well aligned on the semi-molten substrate via the sunken-mode I and the grain boundaries can be therefore greatly reduced. When the graphene islands grow larger, they can be well flattened by metal steps attached to their edges due to the sunken-mode II. Such flat graphene islands can not only facilitate the seamless merge of the graphene islands but also prevent the small interstitial atoms (e.g., carbon precursors) from diffusing into the subsurface beneath the graphene sheets and reduce the formation probability of the defects and new nuclei for the growth of second-layer graphene.
In conclusion, the sinking kinetics of graphene islands on semimolten Cu (111) substrate, including carbon nanoclusters and graphene nanoribbons, were systematically investigated at the atomic scale via classical MD simulations. It was confirmed that C 24 can sink into the semi-molten Cu substrate at temperatures T ≥ 1350 K. However, such a sunken-mode (Sunken-mode I) transformed to the second mode (Sunken-mode II) with the increase of cluster sizes. Thus, large clusters like C 54 and C 96 preferred floating on the upmost surface with its edges passivated by molten Cu atoms. Compared with the liquid substrate, the orientation of the C 24 on semi-molten substrate showed a narrow distribution with the rotation angle concentrating around 0°and 60°owning to the specific sunken-mode and the lattice symmetry of unmolten Cu layers. For graphene nanoribbons, the sunken-mode II was preferred. They also present arched or flattened shapes depending on the metal steps formed. These metal steps, surprisingly, can assist in defects healing during the seamless stitching of graphene nanoribbons instead of creating defects. This study serves to deepen our physical insights into the CVD growth of graphene on semi-molten Cu substrates and explains multiple experimental mysteries. These calculations can also provide theoretic references for the controlled synthesis of large single-crystal graphene. Finally, the modeling approach, including the C-M potential and classic MD developed in this article were confirmed against DFT-MD and can thus be applied to many simulations of graphene growth on similar metals.

METHODS
Four-layer Cu (111) slabs with periodicities of 35.75 × 39.84 × 100 Å and 23.00 × 22.14 × 100 Å were adopted as the metal substrates for the exploration of carbon nanoclusters (i.e., C 24 , C 54 , C 96 ) and graphene nanoribbons, respectively. To mimic the semi-infinite surface, Cu atoms in the bottom layer were fixed during the simulation.
Pure classical molecular dynamics (MD) simulations were performed with a time step of 0.5 fs. Similar to our previous work 51 , the velocity Verlet algorithm and the Berendsen thermostat 52 (see the Supplementary Methods) were adopted in MD simulations. The empirical potential energy surface (PES) employed in this study has already been successfully applied to simulate the growth of carbon nanotubes and graphene on Ni surface 42,51 . More specifically, this new PES is based on the modified second generation of the reactive empirical bond order (REBO2) potential of carbon-carbon interaction 53 , the carbon-metal interaction [54][55][56] and the many-body tight-binding (TB) potential for metal-metal interaction 57 . With 26 parameters, we can further tune the carbon-metal interactions precisely. In particular, the angle-dependent graphene edge-substrate interaction was included in order to simulate the dome-like shape of graphene clusters on a Cu surface 42 .
The benchmarked structures used for the parameters fitting is similar to those used in our previous work 42 . The energies listed in Supplementary Tables 1 and 2 are the calculated formation energies of the benchmark structures. After parameter fitting, three graphene clusters (C 20 , C 21 , and C 24 ) on Cu (111) surface were adopted to test the PES with the results shown in Supplementary Table 3. The low carbon solubility in the Cu substrate has been considered in our potential ( Supplementary Fig. 10). Moreover, to further verify the results of classical MD simulations, DFT-MD simulations were also carried out for comparison with the calculation details presented in the Supplementary Methods.

DATA AVAILABILITY
All data generated and/or analyzed during this study are included in this article and its Supplementary Information file. The raw simulation data and the other datasets are available from the author on reasonable request.

CODE AVAILABILITY
The results were simulated by using a code developed by us. We will consider releasing the code later and the usage of this code requires signing an agreement with the code developer.  6 Stitching process of graphene nanoribbons at 1350 K. a-j Snapshots extracted from the MD trajectory at different times. Cyan and ochre balls represent C and Cu atoms, respectively. The Cu atoms of metal step and the added C atoms are highlighted in yellow and red, respectively.