Model for computing kinetics of the graphene edge epitaxial growth on copper

1.2 MB·

by  Mikhail Khenner

Model for computing kinetics of the graphene edge epitaxial growth on copper
Mikhail Khenner1,2
1Department of Mathematics, Western Kentucky University, Bowling Green, KY 42101
2Applied Physics Institute, Western Kentucky University, Bowling Green, KY 42101
(Dated: July 1, 2016)
A basic kinetic model that incorporates a coupled dynamics of the carbon atoms and dimers on
a copper surface is used to compute growth of a single-layer graphene island. The speed of the
island’sedgeadvancementonCu[111]andCu[100]surfacesiscomputedasafunctionofthegrowth
temperature and pressure. Spatially resolved concentration profiles of the atoms and dimers are
determined, and the contributions provided by these species to the growth speed are discussed.
6 Island growth in the conditions of a thermal cycling is studied.
1
0
2
Journal information: Physical Review E 93, 062806 (2016); the ab initio calculations [4], of the classical BCF-type

n DOI: 10.1103/PhysRevE.93.062806 modeling[9]totwointeractinganddiffusingspeciesthat
u feed growth of the graphene island edge. This results in
J a coupled PDE problem for the concentration fields on

9 the substrate. Both PDEs are also coupled, through the
I. INTRODUCTION
2 boundary conditions, to an ODE for the position of the

island edge.
] Epitaxial growth of a high quality, large area single-
i
c and multi-layer graphene sheets on a transition-metal
s
- substrates is presently in the focus of the research ef- II. MODEL FORMULATION
l forts worldwide, as has been discussed in several re-
r
t view articles [1–3]. Strategies were developed to grow
m As we pointed out in the Introduction, the model is
the graphene sheets with the area of up to 1 cm2 using
aimedatcomputingthevelocity,x˙ (t),ofagrowingedge
. the chemical vapor deposition (CVD) of a hydrocarbons 0
t of a single-layer graphene island. Here x (t) is the posi-
a (such as methane, CH ) on copper, an abundant and in- 0
m 4 tion of the edge, see Fig. 1.
expensive substrate [18–20]. Alongside the experimental
- efforts, modeling of the CVD graphene growth on Cu
d
was also attempted. These studies can be divided into
n
o three groups: ab initio and Kinetic Monte Carlo (KMC)
c methods [4, 6, 7], rate equations [4, 14], and phase-field
[ methods[8]. ByassumingtheanisotropicdiffusiononCu

of the carbon atoms and their anisotropic attachment to
3
v the islands, the authors of the latter reference succeeded
6 in computing the growth of the multi-lobe graphene is-
4 lands on the substrates of different crystallographic ori-
4 entations. However, with the focus of the study on the FIG. 1: Sketch of the graphene island growing on Cu, with
6 morphologies,thekineticmapsofthegrowthratedepen- the atomic events shown.
0
denceonthecontrollableprocessparameters,suchasthe
.
1 pressure and temperature, were not computed. Since the edge grows predominantly by attachment
0 In this paper we describe a simpler, one-dimensional of the carbon atoms and dimers [4–6], let C(x,t) and
6
1 PDE model, whose purpose is to compute the growth C2(x,t) be the concentrations of the carbon atoms ad-
: rate of a single graphene island as a function of three sorbed on Cu and the carbon dimers, respectively. The
v control parameters: the crystallographic orientation of latter result from the assembly of two previously ad-
Xi the substrate, the growth temperature, and the pressure sorbed carbon atoms.
of a gas of the carbon atoms that impinge on the sub- The model is comprised of the following PDEs and
r
a strate. In the manner of Ref. [8], we abstract from the boundary conditions [10].
details of the dissociation of a hydrocarbons, assuming
• Evolution equation for the concentration C on the
that it results in the carbon gas from which the carbon
section of the copper substrate that is not yet cov-
atoms are adsorbed on Cu surface. However, we rec-
ered by the growing graphene island:
ognize, as is pointed out in the ab initio studies, that
besides the carbon atoms there is other diffusing species ∂C ∂2C
that may contribute to the island growth [3] - of which ∂t =Dc∂x2 −χcC2+F, −(cid:96)≤x≤x0(t). (1)
the carbon dimers are thought to be the most impor-
tant [4–6]. Our hybrid model can be seen as an exten- Here D is the carbon atoms diffusivity, F the ad-
c
sion, directly informed by the activation energies from sorption flux, and the sink term −χ C2 describes
c
2
the loss of the carbon atoms due to their assembly ofthecorrespondingconcentrationfromitsequilib-
into the dimers; the kinetics of this loss is recip- rium value [13]. Ω=πa2 is the atomic area, where
rocal in t, e.g. C ∼ 1/χ t, as follows from the a=7×10−9 cm is the radius of the carbon atom.
c
ODE dC/dt = −χ C2. We found that it is not
c
Equations (1)-(7) are made dimensionless by choosing
necessary to include atoms desorption in Eq. (1),
(cid:96), (cid:96)2/D and 1/Ω as the length, time, and concentration
particularly since the desorption rate has not been c
scale, respectively. Keeping the same notations for the
publishedandbecauseevenwithoutthedesorption
dimensionless variables, the dimensionless system reads:
theconcentrationisquitesmall(Fig. 2). Eq. (1)is
awell-posednonlinearPDEwithauniquesolution
for all t>0 [11, 12]. ∂C ∂2C
= −αC2+β, (8)
The boundary conditions for C are: ∂t ∂x2
∂C ∂2C
x = −(cid:96): ∂C =0, (2) ∂t2 = D ∂x22 −δC2+αC2, (9)
∂x
x˙ (t) = −R [C(x (t),t)−ΩC ]− (10)
0 c 0 eq
∂C
x = x (t): −D =K (C−C ). (3) 2R D[C (x (t),t)−ΩC ], x (0)=0.
0 c∂x c eq c2 2 0 eq2 0
∂C ∂C
The first boundary condition states that far from x = −1: =0, 2 =0; (11)
∂x ∂x
thegrapheneisland(atthecenterofthesubstrate)
∂C
thecarbonconcentrationprofileissymmetric. The x = x (t): =R (ΩC −C), (12)
0 ∂x c eq
secondonestatesthatatthegrowingedgetheflux
∂C
of the carbon atoms is proportional to the differ- 2 =R (ΩC −C ).
∂x c2 eq2 2
encebetweentheconcentrationthereandtheequi-
librium concentration [13]; the proportionality pa-
rameter Kc is the kinetic coefficient, which gives a Here the eight parameters are: α=χc(cid:96)2/Dc (the assem-
measure of the ease with which the carbon atoms bly rate of the atoms into the dimers), β = FΩ(cid:96)2/Dc
can attach to the edge. (the adsorption flux of the atoms), δ = κc2(cid:96)2/Dc (the
desorptionrateofthedimers), D =D /D (theratioof
c2 c
• Evolution equation for the concentration C2, also the diffusivities), Rc = Kc(cid:96)/Dc (the attachment rate of
on the section of the substrate not yet covered by the atoms), R =K (cid:96)/D (the attachment rate of the
c2 c2 c2
the growing graphene island: dimers),ΩC ,andΩC (thedimensionlessequilibrium
eq eq2
concentrations).
∂C ∂2C
∂t2 = Dc2 ∂x22 −κc2C2+χcC2, (4) The initial condition for C is taken in the form of a
smoothed step function with a narrow transition, in the
−(cid:96)≤x≤x (t).
0
middle of the interval, from a smaller positive value at
HereD isthediffusivityofthecarbondimers,κ x = −(cid:96) to a larger value at x = x (0) = 0. Zero initial
c2 c2 0
their desorption rate, and the source term χ C2 is conditionforC istaken, e.g. att=0thereisnodimers
c 2
due to assembly of the carbon atoms into dimers. on the substrate.
Notice that through this term the linear Eq. (4) is Apart from the multiple parameters, the system (8)-
one-way coupled to Eq. (1). (12) looks deceptively simple. However, this is the
TheboundaryconditionsforC2 mirrorthoseforC: moving-boundaryproblem,sincethepositionx0(t)ofthe
graphene edge is a priori unknown and must be found
∂C
x = −(cid:96): 2 =0, (5) alongsidewiththeconcentrations. Duetoamovingedge,
∂x anychangeintheconcentrationsgradientsneartheedge
x = x (t): −D ∂C2 =K (C −C ). (6) affects the edge growth speed, and the change in speed
0 c2 ∂x c2 2 eq2 in turn affects the concentrations near the edge and be-
yond. After focusing on the physical parameters in the
HereK istheattachmentcoefficientofthedimers
c2
next section, in Section IV the procedure for the numer-
and C the equilibrium concentration for the
eq2
ical solution of this system of equations is described.
growth by the dimers attachment.
• Equation of the edge growth:
x˙ (t) = −ΩK [C(x (t),t)−C ]− (7)
0 c 0 eq III. PHYSICAL PARAMETERS
2ΩK [C (x (t),t)−C ], x (0)=0.
c2 2 0 eq2 0
This equation states that the edge velocity is the All physical parameters are taken in the Arrhenius
sumofthecontributionsresultingfromtheattach- form, with the most recent and complete, to our knowl-
ment of the atoms and dimers, where each contri- edge, values of the activation energies [4] (see Table I).
bution is proportional to the deviation at the edge The pre-exponential factors are taken proportional to
3
k T/h, where h is Planck’s constant [14]. takes the form:
B
∂U (cid:18) 1 (cid:19)2 ∂2U
D = kBTa2e−EDc/kBT, D = kBTa2e−EDc2/kBT, ∂t = 1+x0(t) ∂ξ2 + (15)
c h c2 h (cid:18) (cid:19)
1+ξ ∂U
k T P x˙ (t) −αU2+β,
χc = Bh e−Eχ/kBT, F = √2πm0k Te−Ead/kBT, 0 1+x0(t) ∂ξ
B ∂V (cid:18) 1 (cid:19)2 ∂2V
κc2 = kBhTe−Eκ/kBT, ∂t = D 1+x0(t) ∂ξ2 + (16)
(cid:18) (cid:19)
k Ta k Ta 1+ξ ∂V
Kc = Bh e−EKc/kBT, Kc2 = Bh e−EKc2/kBT, x˙0(t) 1+x0(t) ∂ξ −δV +αU2,
Ceq = Ω−1e−ECeq/kBT, Ceq2 =Ω−1e−ECeq2/kBT. (13) x˙0(t) = −Rc[U(0,t)−ΩCeq]− (17)
2R D[V(0,t)−ΩC ], x (0)=0.
c2 eq2 0
∂U ∂V
ξ = −1: =0, =0; (18)
∂ξ ∂ξ
Cu surface E E E E E E E E E
Dc Dc2 χ κ Kc Kc2 ad Ceq Ceq2 ∂U
ξ = 0: =(1+x (t))R (ΩC −U), (19)
[111] 0.5 0.49 0.9 1.7 0.71 0.74 0.1 0.87 0.87 ∂ξ 0 c eq
[100] 1.11 0.62 0.59 1.7 1.42 1.07 0.1 0.87 0.87 ∂V
=(1+x (t))R (ΩC −V).
∂ξ 0 c2 eq2
Table 1. Activation energies (in eV).
Values for E , E , E and E were not pub-
κ ad Ceq Ceq2
lishedforgraphenegrowthoncopper. ThusintheTable
1 we adopt the generic values for E , E and E
ad Ceq Ceq2
[15, 16], and for E we adopt a value that partially cur-
κ
tails the otherwise unlimited growth of the dimers con-
centration (caused by the perpetual assembly of the car-
bonatomsintodimers),allowingthecomputationtopro-
ceed until the edge grows over the entire available sub-
strate. This value is large, thus the desorption flux is
small.
CarbongaspressureP isvariedintherange100−600
0
mTorr,m=2×10−23gisthemolecularweightofcarbon,
thetemperatureT isintheinterval873−1273K,andthe
FIG. 2: (Color online). Example dimensionless concentra-
half-width of the substrate (cid:96)=1 mm.
tions of the atoms (solid lines) and dimers (dashed lines),
shown in units of the respective dimensionless equilibrium
concentration, vs. the transformed dimensionless coordinate
along the Cu[111] substrate (see Sec. IV). ξ = 0 (or η = 0)
correspondstothegrowingedge. Timeincreasesinthedirec-
tion shown by the arrow (from the black to the red). Notice
IV. SOLUTION METHODS
thattheconcentrationsarenotfixedatξ=0(bythebound-
ary conditions (19)). Also notice that the dimers concentra-
tion is larger than the one of the atoms, which corroborates
Since solving PDEs on a time-dependent domains is
the findings in the ab initio computations [4].
difficult, we first map the interval −1 ≤ x ≤ x (t) onto
0
a fixed interval −1≤ξ ≤0 for the new space variable ξ,
using the transformation In the transformed system the edge of the graphene
island is at ξ = 0 at all times. However, the true edge
position x (t) is found from Eq. (17) and therefore the
0
ξ = x−x0(t), C(x,t)=U(ξ(x,t),t), (14) kinetics of growth is preserved. The system (15)-(19) is
1+x (t) the initial-boundary value problem for two one-way cou-
0
C (x,t)=V(ξ(x,t),t), pledPDEs(withthevariablecoefficients),whicharealso
2
coupledtothefirst-orderODEforx (t). Forthesolution
0
of this systems we adopted the classical Method of Lines
where U and V are the concentrations of the atoms and (MOL),whichconvertsthePDEsintoODEsbydiscretiz-
dimers on the fixed interval. Then the system (8)-(12) ing the space variable using finite differences. However,
4
with the realistic physical parameters from Sec. III the from the dimers exceeds by a factor of five the one from
computed concentration profiles feature a steep bound- the atoms; on Cu[100], the atoms provide a negligible
ary layer at the growing edge (see Fig. 2, also Figures contribution. Thissupportstherecentconclusionsinthe
8(a,b)). For the prediction of the edge growth rate, it is ab initio [4, 6] and experimental papers [5, 18, 19] that
crucial to resolve these layers with a high accuracy. We thegrapheneedgegrowsprimarilybythedimersattach-
foundthatthisisachievedbyamethodthatwedescribe ment.
below.
First, the space variable ξ is transformed as [17]:
γ
ξ(η)=η+ sinπη, −1≤η ≤0, |γ|<1. (20)
π
Notice that this map is invertible when the absolute
value of the parameter γ is less than one, also η(0) =
0, η(−1) = −1. The purpose of the transformation is
tomapthewould-benon-uniformcomputationalgridon
−1≤ξ ≤0(where atγ <0thegrid pointsareclustered
near ξ = 0) onto a uniform grid on −1 ≤ η ≤ 0. In all
computations we used γ =−0.95.
Next, the final transformed system is discretized in η
using the second-order finite differences with the fixed FIG. 3: (Color online.) Cu[111] surface. (a) Edge speed,
gridspacingshandh/2, andtwoODEsystemsresulting |x0(t)|, vs. the time at P0 =500 mTorr and various temper-
atures. The last point on each curve corresponds to the time
from such discretizations are solved independently and
t atwhichtheedgereachesthemidpointofthesubstrate,
in parallel using the same initial condition. Richardson final
x=−(cid:96)=−0.1 cm. Inset shows the mean edge speed vs. the
interpolation is performed after a fixed number of time
temperature (that was first mapped onto the unit interval),
steps. Thiswayaspatiallyfourth-orderaccuratesolution
wherethemeaniscalculatedoverthetimeintervalfromzero
is obtained on a coarse grid. This solution is then inter-
to t for each curve; the line is the logarithmic fit to the
final
polatedontoafinegridbeforethenextstepistaken. The data shown by squares. (b) Magnification of the T =1073 K
temporal accuracy is achieved automatically by an ODE curve from the panel (a). Insets show separately the compo-
solver. The grid refinement study was performed, which nents of this speed due to the attachment to the edge of the
indicated that using h = 0.0008 results in the needed atoms and dimers.
overall computational accuracy for all parameters values
of interest.
V. RESULTS
WebeginthissectionwiththecomparisonsofFigures3
and4,computedforgraphenegrowthonCu[111]surface,
with the corresponding Figures 5 and 6 for the growth
on Cu[100] surface.
InFigures3(a)and5(a)itcanbeseenthatthegrowth
slows down as the temperature increases, which perhaps
explains better graphene quality and larger islands at
higher growth temperatures [1, 5, 18–20]. From the in-
sets to these Figures, it appears that the slow-down is
FIG. 4: (Color online). Cu[111] surface. Edge speed vs. the
logarithmic. This is a new and important model pre- time at T = 1173 K and various pressures. Inset: the mean
diction, as the quantitative experimental results on the speed vs. the pressure (squares) and the linear fit, also at
growthspeedscalingwiththetemperaturehavenotbeen T =1173 K.
published. At each temperature the speed is nearly a
constant value for the entire duration of the simulation Figures 4 and 6 show the dependencies of the edge
(changing less than 1%, see Figures 3(b) and 5(b)). Also speed on the time and pressure at a fixed temperature.
we noticed that the speed is smaller on Cu[100] and it ThespeedincreaseslinearlywithP . Thisisanotherkey
0
slowlyandmonotonicallydecreaseswithtimeonthissur- model prediction that remains to be supported by the
face, while on Cu[111] surface the curve is S-shaped; the experiment; thequantitativeexperimentaldatawerenot
latter dynamics is somewhat similar to the one shown in published. We remark here that the computed growth
Fig. 2 of Ref. [18]. Attachment of the dimers provides speeds shown in Figures 3 - 6 exceed by a few orders of
the major contribution to the growth speed (see the in- magnitude the speeds that are reported in the experi-
sets). In the case of growth on Cu[111], the contribution mental papers. Values from the experiments seem to be
5
oftheorder10−6−10−5 cm/sforthetemperaturerange back to the parameters at 1273 K and computed until
thatweuseinthecomputations. Weconjecturethatthe the substrate overgrowth by a graphene sheet was com-
discrepancies are primarily due to the larger P values pleted. In Fig. 7 we show the growth speed, and in Fig.
0
used in our computations than the carbon partial pres- 8, the concentrations profiles. First, we notice that the
sures in the experiments, as well as because the adopted growth speed is fully reversible, e.g. after the tempera-
E value is approximate. Since the pressures of the hy- tureisquenchedfrom973Kto1273Kthespeedreturns
ad
drocarbonsortheevaporatedcarbonarenotconsistently toitsvaluepriortothecool-down. Whatisremarkableis
reported in the experimental literature, we took for P the large factor (≈40) by which the speed increases (de-
0
thesetof“growthpressure”(or“chamberpressure”)val- creases) when the temperature is decreased (increased).
ues from Ref. [1]. From the inset of Fig. 4 one can see This value can be directly compared to Fig. 3, which is
thatthespeedoftheorder10−5 cm/swouldresultwhen computedatthesameP andataconstant temperature
0
the fit is extrapolated to P ∼ 13 mTorr (for Cu[111] throughout the entire growth phase. There, the factor
0
surface, see Fig. 6, this value is 21 mTorr). These ex- by which the speed changes is 7.5 when the temperature
trapolatedvaluesareorder-of-magnitudeconsistentwith isdroppedfrom1273Kto973K.Clearly, quenchingthe
those reported in Refs. [5, 18] in conjunction with the temperature down and then up during growth results in
growth rates of the orders that we stated above. a large net increase of the growth speed (notice that the
growth is completed in 1.6 s in Fig. 7 and in 5 s in
Fig. 3). Closer examination shows that this increase is
attributed nearly entirely to the dimers; their concentra-
tionexperiencesawaymoreabruptchange(comparedto
theconcentrationoftheatoms)whenthetemperatureis
quenched up/down. This is shown in Fig. 8, where the
concentrationsareplottedbeforethecool-down,afterthe
cool-down, and after the warm-up. Such response of the
concentrations to the temperature quenches is another
indicator that the dimers are primarily responsible for
the experimentally observed growth kinetics.
FIG. 5: (Color online.) Same as Fig. 3, but for Cu[100]
surface.
FIG.7: (Coloronline). Cu[111]surface,P =500mTorr. The
0
temperature is quenched from 1273 K to 973 K and back.
FIG. 6: (Color online). Same as Fig. 4, but for Cu[100]
surface.
It is common in the experiments to employ thermal
cycles during growth or sharply decrease the tempera-
ture at the very end of the growth phase. This typically
results in better quality of the graphene layer, also its
area is enlarged [19, 22, 23]. Why this happens is not FIG.8: (Coloronline.) Cu[111]surface,P =500mTorr. (a)
0
well understood [21]. Our model is well-suited for giving Atoms concentrations at T = 1273 K and 973 K before the
some insights into this situation. We started the compu- cool-down (dashed magenta line), after the cool-down (solid
tation using the parameters at 1273 K and computed for green line) and after the warm-up (dash-dotted purple line).
some time, then instantaneously switched to the param- (b) Dimers concentrations; the lines coloring is the same as
in (a).
eters at 973 K and computed more, and finally switched
6
It was determined [18–20] that the growth slows down may be also necessary to include the atoms desorption
with time, the more so the closer the graphene islands terminEq. (1)andtheatomsanddimersde-attachment
approach each other [20]. In the cited papers the Cu rates (from the island) into Eq. (7), along with the cor-
crystallographicsurfaceisnotidentifiedthough,itisonly responding source terms in Eqs. (1) and (4). It must
statedthatthegrowthisrealizedonaCufoil. Also,since be noted though, that a time-resolved graphene growth
the observations are made when there is several growing experiments that generate a high-precision data on the
islands, as is always the case, the growth slowdown may growth rates, as well as the matching detailed descrip-
not occur were it was possible to grow a single island. tions of the plethora of the growth conditions and pa-
In our modeling, the minor decrease of the growth speed rameters, are still rare, which presents quite a challenge
is seen for Cu[100] surface, but not for Cu[111] surface. to further tuning the model.
However, it will be fairly straightforward to incorporate
another growing island into the model, which may allow Acknowledgments. The author acknowledges construc-
tomorepreciselydifferentiatebetweenthegrowthmodes tive discussions with V. Dobrokhotov (WKU Applied
on these Cu surfaces. For better predictive capability it Physics Institute).
[1] C. Mattevi, H. Kim, and M. Chhowalla, “A review of tionsofthenonlineardiffusionequationwithanonlinear
chemical vapour deposition of graphene on copper”, J. source”, J. Austral. Math. Soc. Ser: B 39, (1998) 513-
Mater. Chem. 21, (2011) 3324-3334. 527.
[2] A. Reina and J. Kong, “Graphene Growth by CVD [13] M.Uwaha,“Fluctuationandmorphologicalinstabilityof
Methods”, Ch. 7 in: R. Murali (ed.), Graphene Na- stepsinasurfacediffusionfield”,Advancesintheunder-
noelectronics: From Materials to Circuits, Springer standingofcrystalgrowthmechanisms,Ed.T.Nishinaga
Science+Business Media, 2012. et al., Elsevier (1999), pp. 31-45.
[3] H.Tetlow,J.PosthumadeBoer,I.J.Ford,D.D.Vveden- [14] A. Zangwill and D.D. Vvedensky, “Novel growth mech-
sky,J.Coraux,andL.Kantorovich,“Growthofepitaxial anism of epitaxial graphene on metals”, Nano Lett. 11,
graphene: Theory and experiment”, Phys. Reports 542, (2011) 2092-2095.
(2014) 195-295. [15] W. Hong, Z. Suo, and Z. Zhang, “Dynamics of terraces
[4] P. Wu, Y. Zhang, P. Cui, Z. Li, J. Yang, and Z. on a silicon surface due to the combined action of strain
Zhang, “Carbon dimers as the dominant feeding species and electric current”, J. Mech. Phys. Solids 6, (2008)
inepitaxialgrowthandmorphologicalphasetransitionof 267-278.
graphene on different Cu substrates”, Phys. Rev. Lett. [16] G.SchulzeIcking-Konert,M.Giesen,andH.Ibach,“De-
114, (2015) 216102. cay of Cu adatom islands on Cu (111)”, Surf. Sci. 398
[5] K. Celebi, M.T. Cole, J.W. Choi, F. Wyczisk, P. Legag- (1998) 37-48.
neux, N. Rupesinghe, J. Robertson, K.B.K. Teo, and [17] W.F. Spotz and C.F Carey, “Formulation and experi-
H.G.Park,“Evolutionarykineticsofgrapheneformation ments with high-order compact schemes for nonuniform
on copper”, Nano Lett. 13, (2013) 967-974. grids”,Intl.J.Num.Meth.Heat&FluidFlow 8,(1998)
[6] P.Wu,W.Zhang,Z.Li,andJ.Yang,“Mechanismsof 288-303.
graphenegrowthonmetalsurfaces: Theoreticalperspec- [18] X. Li, C.W. Magnuson, A. Venugopal, R.M. Tromp,
tives”, Small 10, (2014) 2136-2150. J.B. Hannon, E.M. Vogel, L. Colombo, and R.S. Ruoff,
[7] P. Gaillard, T. Chanier, L. Henrard, P. Moskovkin, and “Large-area graphene single crystals grown by low-
S. Lucas, “Multiscale simulations of the early stages of pressure chemical vapor deposition of methane on cop-
the growth of graphene on copper”, Surf. Sci. 637638, per”, J. Am. Chem. Soc. 133, (2011) 2816-2819.
(2015) 11-18. [19] Z.Yan,J.Lin,Z.Peng,Z.Sun,Y.Zhu,L.Li,C.Xiang,
[8] E. Meca, J. Lowengrub, H. Kim, C. Mattevi, and V.B. E. LoicSamuel, C. Kittrell, and J.M. Tour, “Toward the
Shenoy,“Epitaxialgraphenegrowthandshapedynamics synthesisofwafer-scalesingle-crystalgrapheneoncopper
oncopper: phase-fieldmodelingandexperiments”,Nano foils”, ACS Nano 6, (2012) 9110-9117.
Lett. 13, (2013) 56925697. [20] Z.-J. Wang, G. Weinberg, Q. Zhang, T. Lunkenbein, A.
[9] W.K.Burton,N.Cabrera,andF.C.Frank,“Thegrowth Klein-Hoffmann,M.Kurnatowska,M.Plodinec,Q.Li,L.
of crystals and the equilibrium structure of their sur- Chi, R. Schloegl, and M.-G. Willinger, “Direct observa-
faces”, Phil. Trans. R. Soc. A 243, (1951) 299-358. tionofgraphenegrowthandassociatedcoppersubstrate
[10] See A.L.-S. Chua, E. Pelucchi, A. Rudra, B. Dwir, E. dynamicsbyinsituscanningelectronmicroscopy”,ACS
Kapon, A. Zangwill, and D.D. Vvedensky, “Theory and Nano 9, (2015) 1506-1519.
experiment of step bunching on misoriented GaAs(001) [21] H.Kim,E.Saiz,M.Chhowalla,andC.Mattevi,“Model-
during metalorganic vapor-phase epitaxy”, Appl. Phys. ingoftheself-limitedgrowthincatalyticchemicalvapor
Lett. 92, (2008) 0113117 for a somewhat similar PDE depositionofgraphene”,NewJ.Phys. 15,(2013)053012.
modelofthemulti-speciescrystalgrowth(notspecialized [22] L. Gao, J.R. Guest, and N.P. Guisinger, “Epitaxial
for particularities of the graphene growth). grapheneonCu(111)”,NanoLett. 10,(2010)3512-3516.
[11] H.A. Levine, “The role of critical exponents in blowup [23] A.A.Koos,A.T.Murdock,P.Nemes-Incze,R.J.Nicholls,
theorems”, SIAM Review 32, (1990) 262-288. A.J.Pollard,S.J.Spencer,A.G.Shard,D.Roy,L.P.Biro,
[12] K. Vijayakumar, “On the integrability and exact solu- and N. Grobert, “Effects of temperature and ammonia
7
flow rate on the chemical vapour deposition growth of 16, (2014) 19446.
nitrogen-doped graphene”, Phys. Chem. Chem. Phys.