J Mol Model (2014) 20:2260 DOI 10.1007/s00894-014-2260-2

ORIGINAL PAPER

Computational study of graphene growth on copper by first-principles and kinetic Monte Carlo calculations Simone Taioli

Received: 31 December 2013 / Accepted: 23 April 2014 # Springer-Verlag Berlin Heidelberg 2014

Abstract In this work the growth of a graphene monolayer on copper substrate, as typically achieved via chemical vapor deposition of propene (C3H6), was investigated by firstprinciples and kinetic Monte Carlo calculations. A comparison between calculated C1s core-level binding energies and electron spectroscopy measurements showed that graphene nucleates from isolated carbon atoms adsorbed on surface defects or sub-superficial layers upon hydrocarbon fragmentation. In this respect, ab initio nudged elastic band simulations yield the energetic barriers characterizing the diffusion of elemental carbon on the Cu(111) surface and atomic carbon uptake by the growing graphene film. Our calculations highlight a strong interaction between the growing film edges and the copper substrate, indicative of the importance of the grain boundaries in the epitaxy process. Furthermore, we used activation energies to compute the reaction rates for the different mechanisms occurring at the carbon–copper interface via harmonic transition state theory. Finally, we simulated the long-time system growth evolution through a kinetic Monte Carlo approach for different temperatures and coverage. Our ab initio and Monte Carlo simulations of the out-ofequilibrium system point towards a growth model strikingly different from that of standard film growth. Graphene growth

This paper belongs to Topical Collection 9th European Conference on Computational Chemistry (EuCo-CC9) S. Taioli (*) Interdisciplinary Laboratory for Computational Science, FBK-Center for Materials and Microsystems and University of Trento, Trento, Italy e-mail: [email protected] S. Taioli Department of Chemistry, University of Bologna, Bologna, Italy S. Taioli Istituto Nazionale di Fisica Nucleare, Sezione di Perugia, Italy

on copper turns out to be a catalytic, thermally-activated process that nucleates from carbon monomers, proceeds by adsorption of carbon atoms, and is not self-limiting. Furthermore, graphene growth seems to be more effective at carbon supersaturation of the surface—a clear fingerprint of a large activation barrier for C attachment. Our growth model and computational results are in good agreement with recent X-ray photoelectron spectroscopy experimental measurements. Keywords Graphene growth . Electron spectroscopy . Ab initio calculations . Nudged elastic band method . Reaction rates . Kinetic Monte Carlo

Introduction Graphene is the most promising candidate material to replace silicon for the development of new electronic devices with unsurpassed characteristics and for the design of novel carbon-based nanocomposites with tailored functionalities [1]. The success of graphene to achieve these goals within a reasonable time-frame depends critically on many factors. First and foremost, a better understanding of the growth mechanisms to control the epitaxy processes on different metallic and semiconductive substrates is needed in order to achieve high-quality large-grain size of the monolayer [2, 3]; secondly, doping and functionalization of graphene have to be improved in order to tune the electronic and optical properties of pristine layers [4, 5] aimed at applications in microelectronics or to reinforce the mechanical properties of polymer composite by graphene intercalation for a specific use [6]. Among the fundamental aspects mentioned above, this work focuses on the computational investigation of the chemical-physical processes involved in graphene growth on metallic substrates, in particular copper. In fact, while

2260, Page 2 of 13

graphene epitaxy can be achieved by thermal decomposition of silicon carbide (SiC) at high temperature under vacuum [7], the most successful preparation technique for obtaining highquality low-cost large-area films is still provided by chemical vapor deposition (CVD) on metals [8]. Indeed, epitaxial graphene on SiC, a very promising and easily manufacturable insulating substrate [9, 10], would probably be the best and most efficient platform for microelectronics applications, as no film transfer is required from the growth substrate to a different one. However, growing graphene on insulating materials requires temperatures in excess of 1,500 °C [7], i.e., significantly higher than those needed for CVD growth on metals (900 °C) [11]. In most cases, these high temperatures are needed to obtain good quality graphene layers as lower temperatures prevent the interface between epitaxial graphene and the silicon carbide substrate from forming. Moreover, in case graphene layers should be transferred from the growth substrate to a different one, no exfoliation procedure is presently available for SiC. On the other hand, dissolving metallic foils after CVD growth is commonly employed to transfer graphene sheets from metals to other substrates [12, 2], producing remarkably long, structurally well-defined graphene. In this regard, many features make graphene epitaxy on metallic substrates appealing, the most notable being the growth dynamics dependence on details of crystal edges, the equivalence between the binding energy of in-plane carbon–carbon bonds (≃7.4 eV per carbon atom) and of graphene edge-metal substrate (≃7 eV per carbon atom); finally, the reversibility of the growth dynamics is a very useful feature [11]. For these reasons, the growth of graphene by CVD on metals, notably copper and nickel [3, 7, 12, 2], is currently the most promising way to achieve low-cost scalable graphene synthesis. Thus, in this work we decided to focus our attention on this technique. CVD consists typically of a reactor in which hydrocarbon precursors, such as C2H4, CH4 and C3H6, are deposited at relatively high temperature (above 900 °C) and relatively low pressure (below 10−3 mbar) on metallic surfaces. The high temperature in the chamber stimulates the catalytic decomposition of precursor gases, hydrogen desorption, diffusion on the surface and, depending on the metallic substrate, dissolution and segregation of decomposed carbon species into the bulk. Lower temperatures, particularly for the case of nickel where carbon solubility is high, tend to give rise to amorphous and carbidic forms. Upon cooling of the heated sample, bulkdissolved carbon emerges and gives rise to surface nucleation and graphene growth. However, in some metals, such as copper, the last steps are negligible. In this case, graphene nucleates from superficial carbon, trapped upon hydrocarbon exposure at high temperature and then diffusing over the substrate surface. Furthermore, low carbon solubility in copper prevents additional growth by segregation during cooling from CVD temperatures, which is,

J Mol Model (2014) 20:2260

at variance, a known problem with nickel substrates leading to non-uniform multilayer films. The growth of single layer graphene on copper foils via CVD was realized by Ruoff and coworkers in 2009 [12], while Hofmann and co-workers recently followed and characterized the entire graphene CVD process [2]. The goal of this paper was to shed some light on this specific set of experimental observations in which graphene is grown on a copper substrate by CVD using propene (C3H6) as precursor gas. In fact, while the experimental procedures are rather simple and well established, atomic-scale mechanisms driving the growth, e.g., nucleation, metal–carbon interface interactions and attachment of carbon atoms to the growing flake, are not well understood [13]. Growth processes are critically dependent on the type of bonding that graphene layer and metal form, and thus each procedure needs a specific study. Nevertheless, a deep understanding of these basic chemical-physical mechanisms at the graphene–substrate interface is crucial to increase the quality and the grain size of graphene layers and the scalability of the overall process. Furthermore, a better comprehension of graphene epitaxy on metals would also help to understand the nucleation steps of carbon nanotubes, which are still poorly understood two decades after their discovery [14]. Indeed, while the growth of nanotubes is unordered and perpendicular to the metallic surface rather than epitaxial, as result of the fact that nanotubes nucleate from nanometer sized catalyst particles determining the tube diameter and electronic character [15, 16], the growth dynamics in these two cases is similar. To describe the growth process our analysis focused on three different directions: first, we performed C1s core-level theoretical analysis for different configurations using a recently implemented method [17–19] for calculating electron spectra in solids; second, we used the nudged elastic band method (NEB), based on density functional theory (DFT), to find activation energies of the fundamental chemical and physical mechanisms on which graphene epitaxy is based; finally, we used these ab initio data to perform kinetic Monte Carlo simulations to investigate the long-time dynamics of the growth that are not accessible by first-principles. By comparing our theoretical data on graphene growth on Cu with electron spectroscopy experiments, we put forward a model showing how graphene growth on metals differs from standard heteroepitaxy and occurs by a catalytic capture of carbon atoms decomposed from precursors and diffusing on the surface. This paper is structured in the following way. “Computational details” briefly outlines the theoretical approaches and ab initio techniques used to modeling the observed phenomena. Electron spectroscopy of the forming film C1s corelevels and NEB simulations to find the most relevant diffusion and capture energy barriers in graphene growth are discussed in “Electron spectroscopy and growth model of graphene on

J Mol Model (2014) 20:2260

copper” and “Nudged elastic band calculations of the process activation energies”, respectively. The ensuing calculations of the reaction rates for the above processes and the kinetic Monte Carlo simulations of the long-time growth dynamics at different carbon monomer concentrations are reported in “Kinetic Monte Carlo of the C–C chemical reactions on copper surface”. Concluding remarks are finally sketched in “Conclusions”.

Computational details C1s core-level spectra of gas-phase hydrocarbons and isolated graphene islands on copper surfaces were performed via the in-house developed ab initio code suite SURPRISES [17–19]. This approach provides an extension of the Fano’s resonant multichannel scattering theory [20], and enables the construction of a continuum wave-function of emitted electrons with appropriate boundary conditions, accurately including the main electron–electron correlation effects. In this method, initially one calculates bound and core-hole excited states of the system under investigation by defining a model Hamiltonian projected onto the Hilbert space spanned by Gaussian basis sets centered on the nuclei. Bound and core-hole states were calculated at Hartree-Fock (HF) level of theory as a starting point, and the accuracy of the calculation was increased subsequently by using standard quantum-chemistry methods, notably configuration interaction (CI). Once the CI Hamiltonian for ground and core-hole states is diagonalized, the continuum wavefunction is obtained separately by solving a Lippmann-Schwinger equation at the same energy with a potential given by a long-range Coulomb interaction between the photo-emitted electron and the remaining excited ion. As above, the Coulomb potential is projected onto the previously used Gaussian functions or onto an enlarged Hilbert space for increasing accuracy. Discrete and continuum orbitals for several open and closed channels corresponding to different electronic configurations interact via the Fano’s resonance mechanism [20], and the projected Hamiltonian needs to be diagonalized on the basis of the vectorial products of the discrete and continuum orbitals. Finally, from the knowledge of the orbitals for the global photoionization process, the photoemission intensity can be calculated in dipole approximation by first-order perturbation theory. The reader can find more details on the theory and calculations of this method in Refs. [18, 17, 19]. In particular, for calculating gas-phase hydrocarbon spectra, a cc-pVDZ basis set, including s, p and d-symmetry Gaussian functions was centered on each carbon and hydrogen atom, and single-double excitations starting from the HF ground and C1s core-hole excited states were performed. Basis sets used in these calculations can be obtained from use of the basis set exchange (BSE) software and the EMSL

Page 3 of 13, 2260

Basis Set Library [21, 22]. For calculating the graphene island C1s core-level spectra adsorbed on copper, cc-pVDZ basis sets were adopted for each carbon atom where the core-hole was created, and for the first and second stars of carbon and copper atoms surrounding it. To lower the overall computational cost, a 6–31G* basis set was used for the carbon and copper atoms beyond the second star of neighbors. Single excitations starting from the HF orbitals for both ground and core-hole excited states in the active space generated by the full basis set were included in the CI procedure. In this case, only the binding energy of two carbon atoms in the graphene island was calculated to model tight and weak interactions between copper and carbon in the growing speckle. DFT calculations were performed using the ab initio totalenergy and molecular dynamics program VASP [23–25] based on the finite temperature implementation of DFT. The ion-electron interaction is described using the projector augmented (PAW) technique [26] with single-particle orbitals expanded in plane waves with an energy cut-off of 400 eV, which ensures convergence of the electronic structure and of the total energy within chemical accuracy (0.01 eV). The non spin-polarized local Density Approximation (LDA) was used to describe the exchange-correlation functional, after accuracy tests that will be discussed further below. LDA functionals present in the VASP library are the Perdew-Zunger (PZ) [27] fit to Ceperley and Alder’s (CA) Monte Carlo calculation [28] of the total energy of the homogeneous electron gas. The Cu(111) surface was modeled by a supercell equal to (5×5) surface unit cells in a five-layer slab, separated from its successive images in the z-direction by a vacuum region of 10 Å to avoid spurious interactions among periodic images. The bottom-most layer of the Cu(111) slab was fixed during energy minimization and in all calculations; all other layers and adsorbates were allowed to relax. Due to the large size of the calculation supercell and after converge tests on the number of k-points, the Cu(111)-(5×5) slab was sampled at the Γ point only. The climbing image nudged elastic band (CI-NEB, from now onwards NEB) method [29] was used for accurate identification of the saddle points of the reactions considered in this work and for calculating the activation energy barriers along the reaction coordinate for each elementary step. NEB is a family of methods able to find the minimum energy reaction pathway (MEP) in the configuration space of the system under investigation by assuming that the working temperature is lower in comparison to the energy barrier. MEP defines a path that connects initial reagents and final products, having orthogonal components of the force equal to zero. Thus, it is a constrained minimization method that can be pursued by any minimization algorithm, starting from an initial guess of the reaction path. In this work we used the steepest descent method for minimizing the ionic positions. The path was discretized in a suitable number of points (images) specific

2260, Page 4 of 13

to the reaction under investigation, from 5 to 12 as described below. The starting mesh points were obtained by linearly interpolating the initial and final optimized ionic positions. To perform accurate NEB simulations, initial and final structures of the path were fully relaxed until forces acting on the atoms were smaller than 0.001 eV/Å . In general, we noticed that NEB simulations converge more rapidly by keeping the starting number of images reasonably small and by adding more images to further NEB runs. Care was used in choosing the repeated images to keep the center of mass fixed in all images. To avoid images sliding one on the top of the other and to keep an evenly spaced grid, a “virtual spring” was introduced in the NEB approach to connect one point of the mesh to the other. In our NEB calculation spring force between images was set to 5.0 eV/Å 2. Furthermore, transition state theory in the harmonic approximation was adopted to calculate the reaction rates characterizing carbon monomer diffusion on the copper surface and carbon uptake by the growing graphene layer. Within this approach, the reaction rates are assumed to be Boltzmann distributed. Finally, an ad-hoc implementation of the kinetic Monte Carlo method (KMC), based on the acceptance method proposed Bortz, Kalos and Lebowitz (BKL) [30] with further improvements for overlayer dynamics by Voter [31], was used to study the long-time dynamics mechanisms characterizing the graphene film epitaxy. An exhaustive description of kinetic Monte Carlo methods can be found in [32] and a relevant application to the CO oxidation at RuO2 (110) of the two-step philosophy followed in this paper is explained in [33]. A typical KMC simulation implementing this approach proceeds in the following steps: (1) At time t, a list of reaction rates kn is used to assess the transition probability for each of the possible jumps that each carbon atom can do (each atom can hop in the up-down and left-right directions finding 0, 1 and 2 carbon atoms representing the fundamental building blocks of graphene flakes). Due to the complexity of the transition patterns with a large number of atoms, the list of possible reaction rates was reduced to the three classes of jump: diffusion, dimer and trimer formation. For hops that are forbidden k=0. (2) Once all the possible events for each carbon atom are identified, an ordered list of events weighted by its occurrence probability is created. In our approach, a total activity function A is defined as the sum of all possible reaction rates at time t. By dividing the reaction rate of a given jump by A, one obtains the transition probability Pn =kn/A for that event. (3) Finally, an event among those possibly occurring is selected randomly using an ad-hoc implementation of the residence-time algorithm [30] as follows. By choosing a random number u on (0:1), a final state m is selected by comparing u with the m probability of hopping to such state: ∑ m−1 n=1 Pn

Computational study of graphene growth on copper by first-principles and kinetic Monte Carlo calculations.

In this work the growth of a graphene monolayer on copper substrate, as typically achieved via chemical vapor deposition of propene (C3H6), was invest...
4MB Sizes 2 Downloads 4 Views