Epitaxial Graphene Growth and Shape Dynamics on Copper: Phase

Oct 22, 2013 - Graphene: One Material, Many Possibilities—Application Difficulties in Biological Systems. Marta Skoda , Ilona Dudek , Anna Jarosz , ...
3 downloads 0 Views 2MB Size
Download PDF
Letter pubs.acs.org/NanoLett

Epitaxial Graphene Growth and Shape Dynamics on Copper: PhaseField Modeling and Experiments Esteban Meca,† John Lowengrub,† Hokwon Kim,‡ Cecilia Mattevi,§ and Vivek B. Shenoy*,∥ †

Department of Mathematics, University of California, Irvine, California 92697-3875, United States CEA, LETI, Minatec Campus, F-38054 Grenoble, France § CASC, Department of Materials, Imperial College London, London SW7 2AZ, United Kingdom ∥ Department of Materials Science and Engineering, University of Pennsylvania, Philadelphia, Pennsylvania 19104-6272, United States ‡

S Supporting Information *

ABSTRACT: The epitaxial growth of graphene on copper foils is a complex process, influenced by thermodynamic, kinetic, and growth parameters, often leading to diverse island shapes including dendrites, squares, stars, hexagons, butterflies, and lobes. Here, we introduce a phase-field model that provides a unified description of these diverse growth morphologies and compare the model results with new experiments. Our model explicitly accounts for the anisotropies in the energies of growing graphene edges, kinetics of attachment of carbon at the edges, and the crystallinity of the underlying copper substrate (through anisotropy in surface diffusion). We show that anisotropic diffusion has a very important, counterintuitive role in the determination of the shape of islands, and we present a “phase diagram” of growth shapes as a function of growth rate for different copper facets. Our results are shown to be in excellent agreement with growth shapes observed for high symmetry facets such as (111) and (001) as well as for high-index surfaces such as (221) and (310). KEYWORDS: Graphene, phase field modeling, epitaxial growth, anisotropic diffusion

A

growth parameters, as such pressure of methane and hydrogen,13−15 has been experimentally investigated, the evolution of the nuclei shape over time has not been fully understood both experimentally and theoretically. As the grain boundary orientations are dictated by the shape of the nuclei before they are stitched together, it is of fundamental importance to understand their time evolution. Models that predict growth shapes as a function of underlying copper crystallinity, gas pressures, and their time evolution up to reaching a continuous film can aid in identifying growth conditions that will lead to high quality graphene with desired properties. From a theoretical perspective, analysis has been restricted to first principles simulations that only consider the very early stages of growth and not the complex nonequilibium evolution of morphology at the mesoscopic level.16−20 Also, complex edge structures and nonintuitive surface reconstructions make predictions based on certain atomistic calculations21 not consistent with recent experimental results.22 In this paper, we consider the growth and shape evolution of micrometer sized graphene islands by developing a continuum phase field model. We focus on tracking the evolution of an initially seeded single growing

mong various transition metals, which can catalyze graphene synthesis, copper1 has attracted great interest owing to its unique capability to selectively grow only one continuous graphene layer1 and the ease with which graphene can be transferred onto other substrates after growth.2 Graphene grown on copper is polycrystalline on a micrometer scale and can be uniform in thickness over more than 98% of the surface. Its polycrystalline nature derives from a growth process consisting of formation of stable nuclei, which grow and merge with each other at the “grain” boundaries. Since the orientation and density of the grain boundaries can significantly alter the electrical, mechanical, and thermal properties of graphene,3−5 engineering the density and orientation of twodimensional grain boundaries can provide a means to engineer polycrystalline layers for different applications. The development of methods to lower the density of grain boundaries is also of considerable interest. In particular, to tailor the grain boundaries, it is essential to achieve control over the nucleation sites and nuclei shapes in relation to the underlying crystal symmetry of copper. Indeed a clear correlation between the shape of the graphene nuclei and the symmetry of the crystalline substrates has started to emerge.6 For instance the 4-fold symmetry of Cu(100)7−10 facets lead to the formation of 4-lobed nuclei shape (4-fold), while the 6-fold symmetry Cu(111) facet lead to 6-lobes (6-fold).11,12 In addition, while the dependence of the nuclei shape from the © XXXX American Chemical Society

Received: September 11, 2013 Revised: October 8, 2013

A

dx.doi.org/10.1021/nl4033928 | Nano Lett. XXXX, XXX, XXX−XXX

Nano Letters

Letter

crystal. Our model includes the anisotropy in the edge energies of the nuclei, the kinetic barriers that the diffusing species encounter as they attach to a growing edge, and the anisotropy in diffusion reflecting the crystallinity of the underlying copper facets. The model allows us to quantitatively study the competition between these competing factors and to derive a phase diagram as a function of material and growth parameters. We compare our simulations with experiments on both highsymmetry and high-index copper facets. Model. The starting point of our modeling is local mass conservation of diffusing species, which leads to the following equation for the concentration, c, of “adatoms”:23 c ∂tc = ∇·(D∇c) + F − (1) τ

zigzag and armchair orientations of the growing nuclei. Denoting ψ to represent the orientation of the high-symmetry direction of the copper facet (that is used to define the diffusion tensor) with the x-axis, the diffusion tensor can be written as

where D is a diffusion tensor, the second term in the right-hand side corresponds to the rate of deposition of adatoms, and the third term accounts for desorption with a time constant τ. The equilibrium concentration of adatoms at the edge of a growing domain is given by

Note that, for facets such as (111) and (100), D|| = D⊥, while D|| ≠ D⊥ for high-index orientations such as (221) and (310). We define two new parameters from D|| and D⊥:

0 ceq = ceq [1 + Γ(ξs(θ ) + ξs″(θ ))κ ]

⎛D D = R(ψ )⎜⎜ ⎝0

(7)

where D|| denotes the largest eigenvalue of the diffusion tensor (without loss of generality) and the rotation matrix R(ψ) is given as: ⎛ cos ψ −sin ψ ⎞ R(ψ ) = ⎜ ⎟ ⎝ sin ψ cos ψ ⎠

D̅ =

(2)

1 (D + D⊥); 2

(8)

δ=

D − D⊥ D + D⊥

(9)

Scaling space with a given characteristic length l and time with l2/D̅ , and defining u = (c − c0eq)Ω, we obtain the following nondimensional system: u ∂tu = ∇·(D̃ ∇u) + f − τv (10)

where c0eq is the equilibrium concentration for a flat interface, κ is the curvature of the domain, and ξs(θ) + ξs″(θ) is proportional to the stiffness of the domain edge that depends on θ, the orientation of the normal to the edge. The constant Γ (that gives a measure of the edge energy) is defined in such a way that, in the absence of anisotropy, ξs(θ) = 1. Note that the anisotropy in the edge energy has been calculated for graphene flakes.21 Under nonequilibrium growth conditions, the concentration of adatoms at the edge of the growing domain will not be equal to equilibrium concentration. If we define c+ (c−) as the concentration directly ahead (behind) the edge of the domain, mass conservation at the edge leads to: v n ·(D∇c)+ − n ·(D∇c)− = n − vn(c+ − c −) (3) Ω

vn = n ·(D̃ ∇u| + − D̃ ∇u|− )

(11)

uI = d0(ξs(θ ) + ξs″(θ ))κ + β(θ )vn

(12)

where f = l (Ω/D̅ )(F − τv = τD̅ /l , d0 = and β(θ) = Ωβ̃ (θ)D̅ /l. Also, we define β(θ) = β ξk(θ), that is, the product of the kinetic coefficient with the angle dependent part. As with the step energy, ξk(θ) = 1 in the isotropic case. Finally, in the following, we will omit the tilde of the rescaled diffusion tensor. To numerically solve the above moving-boundary problem, we make use of a phase-field model shown in Figure 1 and described in the Supporting Information (SI). The phase-field equations are solved using the adaptive nonlinear multigrid algorithm and the solver BSAM.25 Equations are discretized using finite differences, and we use an implicit Crank−Nicolson scheme in time. A nonlinear multigrid method is used to solve the discrete system (see SI for details), and the mesh adapts near the moving boundary to allow to allow for merging of side branches of growing dendritic lobes (see Figure 1 above and Figure S2 of the SI for an example of grid adaptation). The numerical simulation of a complex interface evolution problem such as this one is only feasible through an adaptive grid scheme, like the one adopted here. Numerical Results. Taking the values for the capillary parameter, the kinetic coefficient, and the anisotropies from Table 1, and selecting the phase-field parameters accordingly (see SI), we study the growth of a crystallite from a small circular nucleus, for a given deposition flux f and diffusion anisotropy δ. The goal is to study the interaction of diffusion anisotropy with different values of the deposition flux, for 4-fold and 6-fold symmetry. The specific values have been chosen to obtain a representative portrait of the different morphologies. We start by describing the evolution in two extreme cases; one corresponds to isotropic diffusion (δ = 0), and the other corresponds to a very anisotropic diffusion, D|| = 3 D⊥ (δ = 2

where Ω is the area occupied by an adatom. The additional equation necessary to solve for the adatom density can be obtained by invoking linear attachment kinetics, which gives vn = k+(c+ − ceq) + k −(c − − ceq) (4) Ω where the kinetic constants k± give a measure of the ease with which adatoms can attach to a growing edge. Since attachment can, in general, involve breaking of existing bonds between the graphene edge and the underlying substrate, these constants will depend on the edge orientation. Experiments suggest the graphene edges are “permeable” to the movement of adatoms;24 we can therefore assume that c+ = c− = cI, which leads, by combining eqs 3 and 4, to the following boundary conditions at the growing edge: 0 0 cI = ceq + ceq Γ(ξs(θ ) + ξs″(θ ))κ + β (̃ θ )vn

0⎞ ⎟⎟R(ψ )−1 D⊥ ⎠

(5)

vn = n ·(D∇c)+ − n ·(D∇c)− (6) Ω where β̃(θ) = ((k+ + k−)Ω)−1 is a kinetic coefficient. In general, the graphene nuclei may not be perfectly oriented with respect to the high-symmetry directions of the underlying copper surface. Since we consider the growth and shape evolution of individual graphene nuclei in this work, we align the coordinate axes such that the x-and y-axes coincide with the B

c0eq/τ),

2

ΩΓc0eq/l,

dx.doi.org/10.1021/nl4033928 | Nano Lett. XXXX, XXX, XXX−XXX

Nano Letters

Letter

Figure 1. Schematic representation of graphene growth. We assume a constant deposition flux, F, and desorption time constant, τv. The diffusion on the Cu facet can, in general, be anisotropic reflecting its symmetry. In the upper-left corner we present a typical numerical mesh for this simulation, and in the box there are given the equations for the phase-field model; see SI for details on the parameters and the auxiliary functions.

Table 1. Parameters Used To Compute the Phase-Field Constants and To Build the Anisotropy Functions parameter

symbol

value

capillary parameter reference kinetic coefficient reference deposition flux kinetic anisotropy (n-fold) step energy anisotropy (n-fold)

d0 β0 f0 εk,n εs,n

6 × 10−4 5.54 × 10−3 0.677 0.08 0.001

Figure 2. Evolution of selected cases. On top, 4-fold case, for δ = 0 and δ = 0.5. Similarly for the bottom two rows, 6-fold case. Flux fixed to f = 4f 0. The angle is ψ = 0 for the 4-fold case and ψ = π/6 for the 6-fold. There is a small line energy anisotropy, εs,6 = 0.001, and the kinetic anisotropy has a strength εk,n = 0.08; the interface thickness is ε = 0.00125 in the 4-fold case and ε = 0.0025 in the 6-fold case.

polygons with concave sides, the anisotropic figures have taken more exotic shapes, with constrictions in the direction of fast diffusion. To elucidate further the dynamics previously described, we have performed a systematic study of the interaction of diffusion anisotropy with the deposition flux. In all of the examples we have focused on the kinetic anisotropy21 and reduced the step energy anisotropy to a small value. We have selected two different symmetries for the kinetic anisotropy, 4fold and a 6-fold symmetry. In Figure 3 we show the results for 4-fold kinetic anisotropy. As we can see, for islands of comparable size we obtain very different results, which depend of course in the history that these have undergone. For smaller values of flux, f, we still can see the main dendrites and an incipient sidebranching, which looks very irregular, as expected from the theory and previous results for this model; see the discussion below. For higher values of f we see very regular square shapes, which become closer to the kinetic Wulff shape, a convex square, the higher the value of f is. We have taken values of f that show clearly this transition. On the other axis we have four values of diffusion anisotropy, δ. The first one is δ = 0, intended to be a template for comparison with the isotropic case, and then a small value, δ = 0.08, to see how critical the effect of increasing δ is. Finally, we have two large values, δ = 0.33 and δ = 0.5, which correspond to D|| = 2D⊥ and D|| = 3D⊥, respectively. The principal axes of diffusion are not rotated, exactly as in Figure 2. The low δ and high f dynamics is dominated by the merging of very developed dendrites with a strong sidebranching, as in Figure 2, but as we reduce the value of f we find that for high values of δ there is not a clear dendrite in the direction of small kinetic coefficient,

0.5). This evolution is portrayed in Figure 2, for 4-fold symmetry without rotation of the axes and for 6-fold symmetry with the principal diffusion axes rotated π/6. In this and in the subsequent examples, the initial condition is a small circular island with radius r = 0.3 and zero supersaturation everywhere. At first, for t = 0.25 we have a small island which has been shaped mostly by diffusion. We see that it largely retains a circular shape in the absence of diffusion anisotropy, but in the anisotropic cases it takes the shape of a more or less distorted ellipse, with its major axis pointing in the direction of fast diffusion. For t = 0.35 the island has approximately doubled its area, and we see the dendrites budding in the directions in which the kinetic coefficient is the smallest. At t = 0.40 we see that we have a dramatic change in the evolution. On the one hand, we observe that the time scale of the evolution has changed, as now it is much faster. We see that some dendrites have developed from the original buds. While in the cases without diffusion anisotropy the islands have evolved into symmetric shapes, with dendrites that become fat very soon, we see that when diffusion is anisotropic the dendrites in the direction of higher diffusion are greatly suppressed. Finally, at t = 0.45 the islands have already reached a large size, which we will compare with the results found for different values of the deposition flux below. We see that most of the dendrites have already merged with each other, as sidebranching has developed. We see still small parcels of substrate within the island in certain cases, which correspond mostly to the places where the grooves of the dendrites were. Also, while the cases with diffusion isotropy have taken the shape of star C

dx.doi.org/10.1021/nl4033928 | Nano Lett. XXXX, XXX, XXX−XXX

Nano Letters

Letter

Figure 3. Summary of the shapes of islands of comparable size for different values of the diffusion anisotropy parameter δ and the deposition flux f. The angle is ψ = 0. There is a small line energy anisotropy, εs,6 = 0.001, and the kinetic anisotropy has a strength εk,4 = 0.08. The interface thickness in this case is ε = 0.00125.

Figure 4. Summary of the shapes of islands of comparable size for different values of the diffusion anisotropy parameter δ and the deposition flux f. The angle is ψ = π/6. There is a small line energy anisotropy, εs,6 = 0.001, and the kinetic anisotropy has a strength εk,6 = 0.08. The interface thickness has a value ε = 0.0025.

but in fact there are many dendrites. This is due to the fact that the much slower diffusion at the top and the bottom of the island make the island very unstable at early stages of development, and hence a Mullins−Sekerka instability26 develops, thus creating many dendrites that compete with each other. For the smallest value of f and for δ = 0.5 we see this in a very clear way. The top and the bottom of the island have grown many dendrites, whereas the sides remain essentially unperturbed. We see this example at a smaller scale for δ = 0.08, since we see that the sidebranching in the direction of fast diffusion is less developed than in the other direction. We can also observe a general trend: The dendrites, that in all cases develop in the direction of small kinetic coefficient, become more pointed as the flux is decreased, for islands of comparable size. We have performed a similar analysis for the 6-fold case; the results are displayed in Figure 4. In that case we have the same values of flux, f, and anisotropy, δ, as in Figure 3, and as in Figure 2 the diffusion axes have been rotated π/6 in the anisotropic case. We observe similar features to the ones present in the 4-fold case, but now the fast diffusion axis is in one of the small kinetic coefficient directions, that is, in one of the directions where dendrites would normally develop. As we can see, the effect of the suppression of growth in the direction of fast diffusion is stronger the higher the value of δ is. We also see that the shape of the island in the δ = 0.5 case transitions from an approximately square shape for the highest value of f to a hourglass shape and finally to a 4-fold shape as the value of f is reduced. We also see again that the different dendrites merge into an island through the growth of sidebranching, which occurs clearly in the slow diffusion direction. As before, this effect is also observed in the small anisotropy (δ = 0.08) case. While in Figure 3 we saw a competition stage between different cells for small values of f and high values of δ, this is not observed here. We believe that this is due to the different symmetries, and the fact that the earlier development of the

kinetic dendrites, which are closer and hence do not allow so many cells in between to be developed. The results for alternatively rotated diffusion axes are shown in the SI (see Figure S1). Comparison to Experiments. Graphene chemical vapor deposition was performed in a hot wall, tube furnace system at low pressures. A 25 μm thick Cu foil was used as substrate for the synthesis of graphene. The foil was electropolished (see SI for details) to remove concave grooves on the Cu surface left by the flat rolling process utilized to manufacture the foils and any other residual roughness.17,27 This is critical for enabling the growth of single layer graphene. The electropolished Cu foil was then placed inside a quartz tube (3.6 cm in diameter and 1 m length) that was loaded into the horizontal cylindrical furnace (heating zone length, 47 cm). The quartz tube was evacuated using a rotary pump with a base pressure of 10−3 mbar, and 100 sccm of Ar were introduced into the chamber for 30 min to replace air before heating the furnace to annealing/ growth temperature of 1030 °C. During the annealing stage, the flow rate of H2 gas was 5 sccm. In the growth stage, the H2 and CH4 flow rates were regulated to systematically investigate the effect of the ratio of the partial pressures. At the end of growth stage, the CH4 was shut off, and the furnace was allowed to cool down naturally with a hydrogen and argon flow. Figure 5 shows the comparison between the experimentally observed graphene nuclei shapes and the theoretical predictions based on the phase field simulations. Different nuclei shapes can be obtained on different copper facets at different pressures of methane and hydrogen. In particular, as we progressively increase the ratio between methane and hydrogen partial pressures (which we simulate by varying the adatom flux), the shape of the nuclei evolve from compact (parts a,d,g,j) to dendritic (parts c,f,i,l).13,28 The graphene nuclei on Cu (100) exhibit square/rectangular shapes (Figure 5a) which evolve into four lobe shapes with wide dendritic extensions (Figure 5b) while increasing the methane partial pressure and decreasing the hydrogen partial pressure, or by reducing the adatom flux in D

dx.doi.org/10.1021/nl4033928 | Nano Lett. XXXX, XXX, XXX−XXX

Nano Letters

Letter

Figure 5. Graphene nuclei shapes obtained on different copper facets at different partial pressures of methane and hydrogen at the temperature of 1030 °C. Parts a, b, and c refer to Cu (100) obtained after 6, 4, and 2 min of exposure to methane and hydrogen gases with partial pressure ratios of 3 × 10−4, 0.028, and 0.81 respectively; the scale bars are 10, 4, and 2 μm, respectively. Parts d, e, and f refer to Cu (310) obtained after 20, 10, and 2 min of exposure to methane and hydrogen gases, with partial pressure ratios of 3 × 10−4, 7.7 × 10−3, and 0.028, respectively; the scale bars are 10, 10, and 4 μm, respectively. Parts g, h, and i refer to Cu (221) obtained after 7, 6, and 3 min of exposure to methane and hydrogen, with partial pressure ratios of 1.6 × 10−5, 0.028, and 0.035 respectively; the scale bars are 4, 4, and 2 μm, respectively. Parts j, l, and m refer to Cu (111) obtained after 7, 4, and 3 min of exposure to methane and hydrogen, with partial pressure ratios of of 3 × 10−6, 1.6 × 10−5, and 0.026, respectively; the scale bars are 4, 2, and 2 μm, respectively. Corresponding shapes from simulations along with the parameters used are presented next to the experimental images.

(310) is reported here for the first time. Since this facet lacks a 4 or 6-fold axis of rotation, we expect the diffusion constant to be asymmetric. By choosing n = 6 (see Table 1) and δ = 0.33, we are able to simulate the sequence of shapes observed in experiments. The nuclei on the high-index facet Cu (221) exhibit a nearly rectangular shape with two grooves that give it a four lobe character (Figure 5g) which progressively evolves into a butterfly-like shape (Figure 5h,i) with increasing methane partial pressure and decreasing hydrogen partial pressure. The final shape (Figure 5i), obtained for the highest methane to hydrogen partial pressures ratios considered here, shows a more complex butterfly shape with dendritic contours. The nuclei shapes on this surface, reported here for the first time, reflects the 2-fold symmetry of the underlying Cu facet (221). The phase field simulations can capture all aspects of shape evolution by choosing n = 4 (see Table 1 and SI) and δ = 0.67, reflecting the fact that the diffusion tensor on this 2-fold surface has two distinct eigenvalues. Our simulations are also in excellent agreement with observations of graphene growth on different copper facets reported by other investigators (see Figure S3 of the SI). Conclusions. Using a phase field method that incorporates kinetic anisotropy, step energy anisotropy, and diffusion anisotropy, we have found an interesting interaction of anisotropic diffusion with kinetic anisotropy, exemplified in a transition, for high-index surfaces, from a regime in which the island grows from the sidebranching of the main dendrites to a different regime where there is a competition between different dendrites in the slow diffusion direction. This competition is more pronounced as the deposition rate is decreased, reflecting the fact that it is a diffusion-limited instability. We have shown that the approach we have developed reproduces all of the key

our simulations. Pronounced cusp-like shapes with narrow dendritic extensions are visible in Figure 5c, obtained for the highest methane to hydrogen partial pressures ratio considered in this study. Therefore, as previously observed,7,9,10 the shape of the nuclei changes for different gas pressures but preserves 4fold symmetry which reflects the underlying crystallinity of Cu. These results are well-reproduced by our phase-field simulations with 6-fold symmetry in the edge and kinetic anisotropies and δ = 0 as the diffusion tensor for this facet with a 4-fold symmetry axes is isotropic. The nuclei on Cu (111) exhibit clear hexagonal shapes (Figure 5j), with increasing dendritic contours along the six vertices (Figure 5k,l) with an increase in methane partial pressure and decrease in hydrogen partial pressure. The final shape (Figure 5l), obtained for the highest methane to hydrogen partial pressure ratio considered in this study, shows a wide dendritic contour. As previously observed,14 the nuclei on Cu (111) shows a hexagonal symmetry at different growth conditions reflecting the symmetry of the underlying Cu facet. The phase-field simulations with 6-fold symmetry in the edge and kinetic anisotropies and δ = 0 (since the diffusion tensor for this facet is symmetric) reproduce the trends observed in experiments very well. Next we consider cases where the diffusion tensor of the facet is anisotropic. Graphene nuclei on one such case, Cu (310), exhibit a polygonal shape with vertex angles close to 120° (Figure 5d) which evolves into a nearly hexagonal shape (Figure 5e) while increasing the methane partial pressure and decreasing the hydrogen partial pressure. The vertices of the hexagon then evolve into a shape with six cusps, two of which are shallow compared to the others (Figure 5f) for the highest methane to hydrogen partial pressures ratio considered in this study. This complex sequence of shape transitions on Cu facet E

dx.doi.org/10.1021/nl4033928 | Nano Lett. XXXX, XXX, XXX−XXX

Nano Letters

Letter

(13) Wu, B.; Geng, D.; Xu, Z.; Guo, Y.; Huang, L.; Xue, Y.; Chen, J.; Yu, G.; Liu, Y. NPG Asia Mater. 2013, 5, e36. (14) Zhang, Y.; Zhang, L.; Kim, P.; Ge, M.; Li, Z.; Zhou, C. Nano Lett. 2012, 12, 2810−2816. (15) Fan, L.; Zou, J.; Li, Z.; Li, X.; Wang, K.; Wei, J.; Zhong, M.; Wu, D.; Xu, Z.; Zhu, H. Nanotechnology 2012, 23, 115605. (16) Kim, H.; Saiz, E.; Chhowalla, M.; Mattevi, C. New J. Phys. 2013, 15, 053012. (17) Kim, H.; Mattevi, C.; Calvo, M. R.; Oberg, J. C.; Artiglia, L.; Agnoli, S.; Hirjibehedin, C. F.; Chhowalla, M.; Saiz, E. ACS Nano 2012, 6, 3614−3623. (18) Loginova, E.; Bartelt, N. C.; Feibelman, P. J.; McCarty, K. F. New J. Phys. 2008, 10, 093026. (19) Mattevi, C.; Kim, H.; Chhowalla, M. J. Mater. Chem. 2011, 21, 3324−3334. (20) Wang, L.; Zhang, X.; Chan, H. L. W.; Yan, F.; Ding, F. J. Am. Chem. Soc. 2013, 135, 4476−4482. (21) Artyukhov, V. I.; Liu, Y.; Yakobson, B. I. Proc. Natl. Acad. Sci. 2012, 109, 15136−15140. (22) Li, M.; Hannon, J. B.; Tromp, R. M.; Sun, J.; Li, J.; Shenoy, V. B.; Chason, E. Phys. Rev. B 2013, 88, 041402. (23) Krug, J. In Multiscale Modeling in Epitaxial Growth; Voigt, A., Ed.; ISNM International Series of Numerical Mathematics; Birkhauser: Basel, 2005; Vol. 149; pp 69−95. (24) Li, X.; Magnuson, C. W.; Venugopal, A.; Tromp, R. M.; Hannon, J. B.; Vogel, E. M.; Colombo, L.; Ruoff, R. S. J. Am. Chem. Soc. 2011, 133, 2816−2819. (25) Wise, S.; Kim, J.; Lowengrub, J. J. Comput. Phys. 2007, 226, 414−446. (26) Mullins, W. W.; Sekerka, R. F. J. Appl. Phys. 1963, 34, 323. (27) Luo, Z.; Lu, Y.; Singer, D. W.; Berck, M. E.; Somers, L. A.; Goldsmith, B. R.; Johnson, A. T. C. Chem. Mater. 2011, 23, 1441− 1447. (28) Vlassiouk, I.; Regmi, M.; Fulvio, P. F.; Dai, S.; Datskos, P.; Eres, G.; Smirnov, S. ACS Nano 2011, 5, 6069−6076. (29) Li, Q.; Chou, H.; Zhong, J.-H.; Liu, J.-Y.; Dolocan, A.; Zhang, J.; Zhou, Y.; Ruoff, R. S.; Chen, S.; Cai, W. Nano Lett. 2013, 13, 486− 490.

features of experimentally observed shape transitions observed on different facets. Our model provides a means to explore the large parameter space that governs graphene growth, with results that are complementary to first principles or atomistic simulations, but this study is not without limitations. In particular, while experiments are carried out by chemical vapor deposition of CH4 in the presence of H2, many other intermediate products of the catalytic decomposition can also be present29 during growth. Here we have simply treated the growth of graphene nuclei by means of an effective “adatom” flux. Therefore, more effort is needed to extend this simple one-species model to include the dynamics of many ephemeral intermediate states.

■

ASSOCIATED CONTENT

S Supporting Information *

Details on the phase-field model, additional simulations with different symmetries, details on the grid generation, experimental details on electropolishing of Cu foil, and finally an additional comparison of simulations with other observed graphene crystallites reported in the literature. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS E.M. acknowledges partial support from the Balsells Foundation. E.M. and J.L. acknowledge partial support from the NSF Division of Mathematical Sciences through grant DMS1217303. C.M. acknowledges the award of a Royal Society University Research Fellowship by the UK Royal Society. V.B.S. acknowledges support from ARO (W911NF-11-1-0171) and NSF (CMMI 1308396 and DMS 1306179).

■

REFERENCES

(1) Li, X.; Cai, W.; An, J.; Kim, S.; Nah, J.; Yang, D.; Piner, R.; Velamakanni, A.; Jung, I.; Tutuc, E.; Banerjee, S. K.; Colombo, L.; Ruoff, R. S. Science 2009, 324, 1312−1314. (2) Song, J.; Kam, F.-Y.; Png, R.-Q.; Seah, W.-L.; Zhuo, J.-M.; Lim, G.-K.; Ho, P. K. H.; Chua, L.-L. Nat. Nanotechnol. 2013, 8, 356−362. (3) Yazyev, O. V.; Louie, S. G. Phys. Rev. B 2010, 81, 195420. (4) Grantab, R.; Shenoy, V. B.; Ruoff, R. S. Science 2010, 330, 946− 948. (5) Bagri, A.; Kim, S.-P.; Ruoff, R. S.; Shenoy, V. B. Nano Lett. 2011, 11, 3917−3921. (6) Wood, J. D.; Schmucker, S. W.; Lyons, A. S.; Pop, E.; Lyding, J. W. Nano Lett. 2011, 11, 4547−4554. (7) Wofford, J. M.; Nie, S.; McCarty, K. F.; Bartelt, N. C.; Dubon, O. D. Nano Lett. 2010, 10, 4890−4896. (8) Zhang, B.; Lee, W. H.; Piner, R.; Kholmanov, I.; Wu, Y.; Li, H.; Ji, H.; Ruoff, R. S. ACS Nano 2012, 6, 2471−2476. (9) Rasool, H. I.; Song, E. B.; Mecklenburg, M.; Regan, B. C.; Wang, K. L.; Weiller, B. H.; Gimzewski, J. K. J. Am. Chem. Soc. 2011, 133, 12536−12543. (10) Wang, H.; Wang, G.; Bao, P.; Yang, S.; Zhu, W.; Xie, X.; Zhang, W.-J. J. Am. Chem. Soc. 2012, 134, 18476−18476. (11) Jacobberger, R. M.; Arnold, M. S. Chem. Mater. 2013, 25, 871− 877. (12) Nie, S.; Wofford, J. M.; Bartelt, N. C.; Dubon, O. D.; McCarty, K. F. Phys. Rev. B 2011, 84, 155425. F

dx.doi.org/10.1021/nl4033928 | Nano Lett. XXXX, XXX, XXX−XXX