Model and Simulations of the Epitaxial Growth of Graphene on Non-Planar 6H-SiC Surfaces
aa r X i v : . [ c ond - m a t . m e s - h a ll ] S e p Model and Simulations of the Epitaxial Growth of Graphene on Non-Planar 6H-SiCSurfaces
Fan Ming and Andrew Zangwill ∗ School of Physics, Georgia Institute of Technology, Atlanta, GA 30332, USA (Dated: March 5, 2018)We study step flow growth of epitaxial graphene on 6H-SiC using a one dimensional kineticMonte Carlo model. The model parameters are effective energy barriers for the nucleation andpropagation of graphene at the SiC steps. When the model is applied to graphene growth on vicinalsurfaces, a strip width distribution is used to characterize the surface morphology. Additional kineticprocesses are included to study graphene growth on SiC nano-facets. Our main result is that theoriginal nano-facet is fractured into several nano-facets during graphene growth. This phenomenonis characterized by the angle at which the fractured nano-facet is oriented with respect to thebasal plane. The distribution of this angle across the surface is found to be related to the stripwidth distribution for vicinal surfaces. As the terrace propagation barrier decreases, the fractureangle distribution changes continously from two-sided Gaussian to one-sided power-law. Using thisdistribution, it will be possible to extract energy barriers from experiments and interpret the growthmorphology quantitatively.
PACS numbers: 68.35.-p, 68.55.-a, 81.15.Aa
I. INTRODUCTION
Epitaxial graphene grown by thermal decomposition ofSiC is an attractive enabling technology for integratinggraphene into silicon microelectronics . However, scal-able production of high quality graphene film is still chal-lenging. Progress in this field requires understanding ofthe growth mechanism. Unfortunately, most experimen-tal results on epitaxial graphene growth are qualitative.The same is true for the interpretation of these results.First-principle approaches to this subject are difficultdue to the complexity of the graphene/SiC interface and an unknown step structure. Previously, the presentauthors introduced a simple one dimensional kineticMonte Carlo (KMC) model to study the epitaxial growthof graphene on vicinal surfaces of 6H-SiC . Based on astudy of graphene strip width distributions, the simula-tion results showed two distinct growth regimes domi-nated by “coalescence” processes and “climb-over” pro-cesses, respectively.Besides growth on SiC vicinal substrates, it has beenshown recently that graphene grows spontaneously onnano-facetted SiC substrates . In this context, a“nano-facet” refers to a stable structure composed ofseveral closely-spaced triple bilayer steps on the flat SiCsurface. Suprisingly, perhaps, growth on a nano-facettedsurface often leads to better quality graphene films on theadjacent terrace . Growth on nano-facets also plays animportant role in graphene ribbon growth . The purposeof this paper is to extend our previous study for vicinalsubstrates to include growth on nano-facetted substrates.Commercially available SiC substrates exhibt roughand scratched surfaces. An effective method to produceatomically flat vicinal surfaces is high temperature H etching. For on-axis 6H-SiC substrates, this method pro-duces vicinal surfaces with triple bilayer steps. This typeof microstep formation is related to etching kinetics and energy differences between different basal planes . Thiscontrasts with off-axis 6H-SiC substrates where the H etched surface often shows periodic structures of nano-facets. Depending on whether the substrate is tilted to-ward the h i or h i direction, the nano-facet makesan angle of ≈ ◦ or 13 − ◦ from the basal plane .Besides H etching, controled nano-facets can be achievedby direct plasma etching of SiC surfaces .Nano-facets can also form spontaneously when SiC isheated close to the graphitization temperature in non-vacuum environments. This type of step bunching is ob-served for both on- and off-axis SiC substrates. At highergrowth pressure, a reduced silicon sublimation rate leadsto a higher graphitization temperature. It is likely thatat these elevated temperatures, the SiC vicinal surfaceof triple bilayer steps is further reconstructed into nano-facets to minimize the surface free energy. The forma-tion of graphene begins at these nano-facets where siliconatoms have fewer bonds . Given the previous discus-sion, SiC step bunching to produce nano-facets may oc-cur before graphene formation. This paper only considersthe formation of graphene on pre-existing nano-facets re-gardless of their origin. We also neglect the differencebetween (1¯10 n ) and (11¯2 n ) nano-facets.The structure of this paper is as follows. Section II is abrief review of our previous results for graphene growthon vicinal surfaces. In section III, we introduce our newmodel for graphene growth on a SiC surface with nano-facets. Section IV focuses on the distribution of fractureangles which can be compared directly with experimentsto extract effective energy barriers. II. STRIP WIDTH DISTRIBUTION
Our original model applies to graphene growth on avicinal surface of 6H-SiC composed of triple bilayer steps.
FIG. 1. Graphene growth kinetic processes on a vicinal sur-face.
From mass conservation, one graphene unit is obtainedby decomposing one SiC step. Step flow growth from SiCtriple bilayer steps produces strips of graphene parallel tothe SiC step edges. This allows us to model the growthas one-dimensional. The detail of the atomic structureof the graphene/SiC interface does not play a role here,so we do not consider it in the model. Two model energybarriers are necessary: one is the barrier to graphene nu-cleation at a step edge (Fig. 1(a)-(b)). The other is thebarrier to graphene strip propagation. Nucleation occursat a rate r nuc = ν exp( − E nuc /kT ), where ν ≈ s − is an attempt frequency and T is the substrate tem-perature. Propagation occurs (Fig. 1(b)-(c)) at a rate r prop = ν exp( − E prop /kT ). The two barriers are “ef-fective” energy parameters which account for Si atomssublimation, C atoms re-crystallization and subsequentgraphene growth along the step edge .After a series of propagation steps, two outcomes arepossible for each graphene strip. One is that a growinggraphene strip runs into another SiC step. We allow thegrowth to continue onto the adjacent upper terrace. Thisis called a “climb-over” process (Fig. 1(d)-(e)). The otheroutcome is that a growing graphene strip meets anothergraphene strip on the upper terrace. In this case, thetwo graphene strips connect to each other to form onecontinuous strip. We call this a “coalescence” process(Fig. 1(f)-(g)). Both climb-over and coalescence producea step which is covered by a continuous graphene layer. Asecond graphene layer can nucleate at such a step underthe first layer. We allow the subsequent growth of thesecond layer graphene to continue in the same way asthe first layer.Our major finding from this model is the cross-over ofthe graphene strip width distribution shown in Fig. 2.This cross-over is the result of competition between theclimb-over and coalescence processes. We define ∆ E = E nuc − E prop and total coverage Θ = P i i Θ i , where Θ i is the graphene coverage of layer i. When the nucleation barrier is small compared to the propagation barrier, i.e.∆ E ≤
0, all SiC steps decompose almost simultaneouslyat the beginning of growth, and the strip width distri-bution is Poisson (Fig. 2(a)). In this case, coalescenceprocesses occur in a short time interval right before thesubstrate is fully covered by graphene. Increasing thenucleation barrier relative to the propagation barrier re-duces the number of graphene strips on the surface andmakes “climb-over” processes more dominant than co-alescences. This decreases the peak magnitude in thestrip width distribution and shifts the peak position tothe right in Fig. 2(b) and (c). Eventually, when ∆ E isbig enough, the step edges become indistinguishable asthe graphene strip grows. The strip width distribution isthen uniform(Fig. 2(d)). III. MODELING GROWTH ON A NANO-FACET
In our modeling, a nano-facet is a group of triple bi-layer SiC steps with a fixed spacing between adjacentsteps (Fig. 3). Hence, the graphene growth proceedsalso in step flow mode which can be modeled as one-dimensional. The nano-facet can have different angleswith respect to the basal plane depending on both theorientation of the substrate and the graphene growthconditions. Adjusting the spacing width in the model al-lows different initial nano-facet angles. We assume thatgrowth starts at the bottom of the nano-facet, convertingone step into one graphene unit. The nano-facet nucle-ation process occurs at a rate r nuc = ν exp( − E nuc /kT )(Fig. 3(a)-(b)). The nucleated graphene unit propa-gates upward along the nano-facet at a rate r prop = ν exp( − E prop /kT ) (Fig. 3(b)-(c)). We keep the value of E prop = 0 because incomplete graphene coverage on thenano-facet is rarely observed experimentally. In otherwords, growth on the nano-facet is very fast.The nano-facet nucleation and propagation processesdefined here are very similar to those for vicinal surfaces.As soon as a nano-facet propagation event occurs, a sec-ond graphene layer can nucleate immediately under thefirst layer (Fig. 3(f)) and the growth of the second layercan continue in the same way as the first layer (Fig. 3(f)-(g)). This contrasts with our previous model for vic-inal surface growth where a second layer of graphenecan only grow at a step covered by a continous graphenestrip (see Fig. 1(h)). When the propagation on the facetreaches the top junction between the nano-facet and theSiC(0001) basal plane, the graphene growth is allowedto continue on the (0001) basal plane but with a slowerpropagation rate r ′ prop = ν exp( − E ′ prop /kT ) (Fig. 3(d)-(e)). This is due to the experimental observation thatseveral layers of graphene often grow on a nano-facet be-fore the graphene growth propagates on the (0001) basalplane . We focus on the early growth stage where thegraphene layers of each macrostep grow independently,and no coalescences or climb-overs are allowed.Fig. 4 shows that a particular choice of KMC model pa- FIG. 2. Graphene strip width distribution ρ ( s ) for different∆ E and total coverages with φ = 0 . ◦ . Different color linescorrespond to Θ = 0 . . . . W = 200. rameters produces a simulated morphology very similarto a transmission electron microscope image of graphenegrowth on a non-planar SiC(0001) surface. The TEMimage was taken of a sample which was prepared at 1325 ◦ C for 90 min, with step heights of 5 −
15 nm. Both theTEM and the simulation image show a sharp ending forall graphene layers at the bottom of the nano-facet, wherethe graphene growth starts. Fig. 4 also shows that, fromthe left to the right side, graphene layers grow continu-ously over the junction between the high index nano-facetand SiC(0001). The top graphene layer is the longest,
FIG. 3. Graphene growth kinetics processes on a nano-facet.The blocks specified by letter A and B are used to calculatethe fracture angle θ .FIG. 4. (a) KMC simulation snapshot. Total surface coverageΘ = 0 . T = 1800 K, E nuc /kT = 7 . E ′ prop /kT is 3.9 and7.1 for the top graphene layer and the interface graphene lay-ers, respectively. (b) Transmission electron microscope image,adapted from Robinson et al. . while the other layers have a similar length. In what fol-lows, we call all the graphene layers other than the topone “interface graphene layers”.Fig. 5 illustrates the formation process for the sur-face morphology in Fig. 4. The original nano-facet isshown in Fig. 5(a). The first graphene layer nucleates atthe bottom junction of the nano-facet, and propagatesquickly up the slope until it reaches the top basal plane.This is shown in Fig. 5(b). Fig. 5(c) shows that the firstgraphene layer continues to grow onto the terrace and asmall triple bilayer nano-facet is fractured from the orig-inal nano-facet (This process was previously shown inFig. 3(d)-(e)). However, the terrace propagation rate isslower than the nano-facet propagation rate. Therefore,Fig. 5(c) also shows that the first graphene layer can-not grow very long on the terrace before the second layeris nucleated under the first layer and quickly covers thelower nano-facet.The terrace propagation rate for the second graphenelayer is slower than the first layer due to the increasing FIG. 5. Formation processes for the surface morphology inFig. 4. difficulty for Si atoms to escape from SiC steps coveredby graphene. This leads to a possible scenario, shown inFig. 5(d), that before the second graphene layer growson the terrace, a third graphene layer is nucleated andcatches up to the second layer growth at the top of thelower nano-facet. Later, when the second graphene layerdoes grow onto the terrace like the first layer, anothertriple bilayer step is decomposed from the lower nano-facet. This process is the same for the third layer. How-ever, assuming that the second and third graphene layershave a similar terrace propagation rate, Fig. 5(e) showsthat a nano-facet of one unit cell high is fractured fromthe original nano-facet by growing the second and thirdgraphene layers together on the terrace.Two small nano-facets are now fractured from the ori-gianl nano-facet, with heights corresponding to one andtwo triple bilayer steps, respectively. This is the sametype of surface morphology seen in Fig. 4(b). We notethat, in our model, graphene grows by decomposing SiCtriple bilayer steps, so the graphene layer always has oneside attached to a nano-facet. To emphasize this, inFig. 5, we draw the graphene layers as slightly curvedat the top of a nano-facet.
IV. FRACTURE ANGLE DISTRIBUTION
To quantify the surface morphology, we study the frac-ture angle θ made by the graphene layers on the nano-facet with respect to the basal plane. In the simulations, θ is quantified as follows. If the nano-facet height is h , wecount the number of units from the edge of the top terrace(marked as letter A in Fig. 3(g)) to the edge of the bottomterrace (letter B in Fig. 3(g)), and define the length as L (e.g. L = 5 in Fig. 3(g)). Then θ ≡ arctan (( h + 1) /L ).The initial θ in the simulations is assumed to be θ = 30 ◦ ,which is close to the nano-facet angle observed experi-mentally. As soon as the top graphene layer growth prop- agates onto the terrace, θ starts to decrease as a func-tion of the growth time. Hence, we can use θ to charac-terize the nano-facet fracturing process due to graphenegrowth. θ can be measured directly by experiments usingthe definition above. Compared to experimental results,the KMC simulation can be used to provide some quanti-tative information about the nucleation and propagationbarriers.The initial surface for our statistical study of growth-induced fracturing is a periodic sequence of basal planesand nano-facets. We fix the height of the nano-facetsto be 22.5 nm, which corresponds to h = 30 SiC triplebilayers. We use a long terrace width and a short totalgrowth time so that no graphene strip coalescence occurs.For statistical purposes, the surface consists of 6 × alternations of SiC(0001) and nano-facets, with a vicinalangle φ = 5 . ◦ . The growth temperature is fixed at 1800K. We also fix the nucleation barrier E nuc /kT = 7 . E ′ prop as the onlyvariable.We now focus on the distribution of fracture angles ata given coverage. For simplicity only, we assume the en-ergy barriers for all graphene layers are the same. Fig. 6shows the distribution ρ ( θ ) of the normalized fracture an-gle tan ( θ ) / tan ( θ ) for different choices of terrace propa-gation barriers E ′ prop with a fixed total coverage Θ = 0 . E f = E nuc − E ′ prop , it is clear from Fig. 6that the distribution of fracture angles shows a cross-overbehavior as ∆ E f increases. For ∆ E f = 0, it is difficultfor graphene to propagate on the terrace. After all thenano-facets are fully covered by graphene, the growth onthe terrace starts almost simultaneously. This leads to aGaussian distribution of fracture angles in Fig. 6. Thisgrowth regime is very similar to graphene growth on vici-nal surfaces for ∆ E = 0, where all the steps are nucleatedat the same time and a Poisson strip width distributionis found.More interestingly, as ∆ E f increases, Fig.6 shows thatthe distribution of fracture angles becomes more andmore one-sided. Eventually, when ∆ E f = 5 .
8, the dis-tribution follows an intriguing power-law form, in whichthe probability is zero for θ values smaller than the min-imum θ given in the figure. We do not fully under-stand the power-law probability distribution. Neverthe-less, the change in fracture angle distribution can stillbe related to the strip width distribution for vicinal sur-faces in Fig. 2. For growth on vicinal surfaces, whenthe nucleation barrier is high, graphene is nucleated atsome steps much earlier than at the others, leading toa shift of peak position to the right side in the stripwidth distribution (Fig. 2(b)). In other words, longerstrips start to dominate the distribution. Similarly, forgrowth on nano-facets, when ∆ E f >
0, the nano-facetsare fractured one after another in a wide time range. Inthe regime where ∆ E f is big enough, the distributionis dominated by smaller values of θ made by longer topgraphene layers. This leads to a one-sided distribution of θ . FIG. 6. Fracture angle distribution for different terrace prop-agation barriers. The red dashed line is a power-law fit, andthe black dashed line is a Gaussian fit.FIG. 7. θ m as a function of KMC time. Dashed curves arethe fit results according to Eq. 1. Red, green, blue and blackdashed curves give E ′′ prop /kT = 2 . , . , . .
0, respec-tively.
Given the above discussion, similar to the strip widthdistribution, the cross-over behavior of the fracture angledistribution can be considered as a competition betweenthe nucleation process at the nano-facet and the propa-gation process on the terrace. When compared to exper-imental observations, the fracture angle distribution canbe used to determine ∆ E f .We define θ m as the most probable θ in the fractureangle distribution. In the following, we give a simpleanalytic treatment to show how θ m evolves as a functionof the growth time. If the experimental growth time is t E , the normalized θ m can be approximated bytan ( θ m )tan ( θ ) = h + 1 h + 1 + ν t E exp( − E ′′ prop /kT ) . (1)In this expression, h is the original nano-facet height, and E ′′ prop is an effective terrace propagation barrier. To de-rive this equation, we do not consider the initial growthtime spent on the nano-facet before the terrace growthoccurs, which gives a small reduction to t E in Eq. 1. E ′′ prop accounts for both the top graphene layer terracepropagation barrier E ′ prop and the nucleation barrier forinterface graphene layers. Increasing the number of in-terface graphene layers will reduce the relative horizontaldistance between the top and the bottom terraces at thefractured nano-facet (i.e. the horizontal distance betweenblock A and B in Fig.3(g)). Both the initial growth timespent on the nano-facet and the nucleation of interfacegraphene layers have the effect of increasing θ . There-fore, it is expected that E ′′ prop > E ′ prop . We treat E ′′ prop as a fitting parameter for our analytic model to see howit is related to E ′ prop .Fig. 7 shows the simulation result for the normal-ized θ m as a function of the KMC time τ K for differ-ent choices of terrace propagation barriers. Here wedefine ∆ = E ′ prop /kT . The KMC time is defined as τ K = t E r prop . Dashed lines are fitted curves according toEq. 1. Despite the simplicity of Eq. 1, the model curvesagree quite well with the KMC simulations for smaller ∆and later growth times. This is likely caused by the ad-justment of initial growth time spent on the nano-facet.Smaller ∆ indicates a faster graphene propagation on theterrace, which leads to a relatively shorter growth timespent on the nano-facet before the terrace growth occurs.Similarly, at later growth times when the initial growthtime on the nano-facet becomes insignificant, the fittedcurve also agrees better.As we increase the terrace propagation barrier E ′ prop ,the fitting parameter E ′′ prop is also found to increase.In fact, we find an excellent linear relationship between E ′ prop and E ′′ prop : E ′′ prop = E ′ prop + 0 .
03 eV. The rela-tively small correction to E ′ prop suggests that the majorcontribution of E ′′ prop still comes from E ′ prop for the topgraphene layer. The contribution from the nucleationbarrier for interface graphene layers acts as an additionalenergy barrier ∆ E prop = 0 .
03 eV. Therefore, with Eq. 1,the time evolution of the fracture angle θ m measured inan experiment can be used to obtain a good estimate of E ′ prop . Combined with the result for ∆ E f from the anal-ysis of fracture angle distribution, both E nuc and E ′ prop can be extracted experimentally.Most experiments of graphene growth on SiC nano-facets are done at a temperature between 1600 K and1800 K. We tested ∆ E prop in this range, but no tem-perature dependence was found. ∆ E prop also does notseem to be dependent on the original nano-facet height h . As we increase the nano-facet height from 22.5 nmto 100 nm, the change of ∆ E prop is less than 10%, wellwithin the accuracy of E ′′ prop as obtained from Fig. 7.This suggests that ∆ E prop is only a function of the en-ergy barriers in the model, and very weakly dependent onthe growth temperature or the initial nano-facet heightin the experimental parameter range. V. CONCLUSION
We have presented a model for graphene growth onnon-planar nano-facetted 6H-SiC substrates based on aprevious model of growth on vicinal substrates. The sim-ulation produces a surface morphology very similar toobservations. A description of the formation process forthis type of surface morphology is also provided. Forgraphene growth on nano-facetted SiC substrates, a frac- ture angle can be used to characterize the growth-inducedfracture of a nano-facet. The distribution of fracture an-gles is found to be related to the graphene strip width dis-tribution for growth on vicinal substrates. As the terracepropagation barrier decreases, the fracture angle distri-bution deviates from a Gaussian and eventually becomesa power-law. Our analytic result for the most probablefracture angle agrees well with the KMC simulations.
VI. ACKNOWLEDGEMENT
We thank Ed Conrad, Phil First and Walt de Heer foruseful discussions. Fan Ming was supported by the MR-SEC program of the National Science Foundation underGrant No. DMR-0820382. ∗ [email protected] C. Berger, Z. M. Song, T. B. Li, X. B. Li, A. Y. Ogbazghi,R. Feng, Z. T. Dai, A. N. Marchenkov, E. H. Conrad, P. N.First, and W. A. de Heer, J. Phys. Chem. B , 19912(2004). P. N. First, W. A. de Heer, T. Seyller, C. Berger, J. A.Stroscio, and J. S. Moon, MRS. Bull. , 296 (2010). J. Hass, W. A. de Heer, and E. H. Conrad, J. Phys.: Cond.Matter , 323202 (2008). R. M. Tromp and J. B. Hannon, Phys. Rev. Lett. ,106104 (2009). F. Ming and A. Zangwill, arXiv:1011.4096 (2011). K. V. Emtsev, A. Bostwick, K. Horn, J. Jobst, G. L. Kel-logg, L. Ley, J. L. McChesney, T. Ohta, S. A. Reshanov,J. Rohrl, E. Rotenberg, A. K. Schmid, D. Waldmann, H. B.Weber, and T. Seyller, Nat. Mater. , 203 (2009). C. Virojanadara, R. Yakimova, A. A. Zakharov, and L. I.Johansson, J. Phys. D: Appl. Phys. , 374010 (2010). J. Robinson, X. J. Weng, K. Trumbull, R. Cavalero,M. Wetherington, E. Frantz, M. LaBella, Z. Hughes,M. Fanton, and D. Snyder, Acs Nano , 153 (2010). M. Sprinkle, M. Ruan, Y. Hu, J. Hankinson, M. Rubio-Roy, B. Zhang, X. Wu, C. Berger, and W. A. de Heer,Nat. Nanotech. , 727 (2010). K. Hayashi, K. Morita, S. Mizuno, H. Tochihara, andS. Tanaka, Surf. Sci. , 566 (2009). A. Nakajima, H. Yokoya, Y. Furukawa, and H. Yonezu, J.Appl. Phys. , 104919 (2005). S. Nie, C. D. Lee, R. M. Feenstra, Y. Ke, R. P. Devaty,W. J. Choyke, C. K. Inoki, T. S. Kuan, and G. Gu, Surf.Sci. , 2936 (2008). H. Nakagawa, S. Tanaka, and I. Suemune, Phys. Rev. Lett. , 226107 (2003). W. Norimatsu and M. Kusunoki, Physica E , 691 (2010). S. Tanaka, K. Morita, and H. Hibino, Phys. Rev. B ,041406 (2010). T. Ohta, N. C. Bartelt, S. Nie, K. Thurmer, and G. L.Kellogg, Phys. Rev. B81