Symmetry of surface nanopatterns induced by ion-beam sputtering: the role of anisotropic surface diffusion

4.8 MB·

by  Javier Renedo

Symmetry of surface nanopatterns induced by ion-beam sputtering: the role of
anisotropic surface diffusion
Javier Renedo,1 Javier Mun˜oz-Garc´ıa,2 Mario Castro,3 and Rodolfo Cuerno2
1Instituto de Investigaci´on Tecnolo´gica (IIT), Universidad Pontificia Comillas, 28015 Madrid, Spain
2Departamento de Matema´ticas and Grupo Interdisciplinar de Sistemas Complejos (GISC),
Universidad Carlos III de Madrid, 28911 Legan´es, Spain
3GISC and Grupo de Dina´mica No Lineal (DNL), Escuela T´ecnica Superior de Ingenier´ıa (ICAI),
Universidad Pontificia Comillas, 28015 Madrid, Spain
(Dated: January 12, 2016)
6 IonBeamSputtering(IBS)isacost-effectivetechniqueabletoproduceorderednanopatternson
1 the surfaces of different materials. To date, most theoretical studies of this process have focused
0 on systems which become amorphous under irradiation, e.g. semiconductors at room temperature.
2 Thus, in spite of the large amount of experimental work on metals, or more recently on semicon-

n ductors at high temperatures, such experimental contexts have received relatively little theoretical
a attention. These systems are characterized by transport mechanisms, e.g. surface diffusion, which
J are anisotropic as a reflection of the crystalline structure not being overruled by the irradiation.
1 Here,wegeneralizeapreviouscontinuumtheoryofIBSatnormalincidence,inordertoaccountfor
1 anisotropic surface diffusion. We explore systematically our generalized model in order to under-
standtheroleofanisotropyinthespaceorderingpropertiesoftheresultingpatterns. Inparticular,

] we derive a height equationwhichpredicts morphologicaltransitions amonghexagonal andrectan-
l
l gular patterns as a function of system parameters and employ an angular correlation function to
a assessthesepatternsymmetries. Bysuitablychoosingexperimentalconditions,itisfoundthatone
h
might be able to experimentally control the type of order displayed by the patterns produced.
-
s
e PACSnumbers: 79.20.Rf,68.35.Ct,81.16.Rf,05.45.-a
m
.
t I. INTRODUCTION TheBHequationanditsgeneralizationsweresimilarly
a
m derivedasinRef.5,byaddingtogetherphysically-diverse
contributionsintoasingleequationforthetargetheight.
- Ion-beam sputtering (IBS) is a technique employed to
d efficiently nanostructure surfaces:1 a solid target is bom- Alternatively, as shown in Refs. 10 and 11, one can de-
n scribe the dynamics of two different fields, the surface
barded with energetic ions, which erode material induc-
o height and the density of material (e.g. adatoms, adva-
ingself-organizedpatternformationatthetargetsurface.
c cancies) subject to transport at the surface. This ap-
[ There is a wide technological interest in this technique,
proach describes surface dynamics successfully in many
sinceitallowstoobtainorderednanostructureswithcon-
1 different contexts, from granular matter12 to epitaxial
trolled roughness, wavelength, and orientation.2,3 More-
v growth.13 In the IBS context, it enables improvements,9
3 over,itisscalable,cost-efficient,andcanbeusedinmany most notably by coupling different physical mechanisms
3 materials, including semiconductors, metals, and insula-
in a natural way. For instance irradiation is expected
5 tors. Onechallengethatstilllimitsthewidespreaduseof
to influence surface diffusion and be reflected in the
2 IBSisthelackofanunifiedtheoreticalframeworkwhich
corresponding terms in the height equation, typically
0
guides experimental designs.
. as a linear high-order derivative term. However, di-
1 In this regard, continuum models have been relatively rect expansion of Sigmund’s contribution in the erosion
0
successful in describing the dynamical behavior of these velocity to such linear14,15 or nonlinear orders16,17 are
6
1 nanostructures, typically in terms of macroscopic vari- affected by consistency issues with respect to pattern
: ables like the target surface height. For materials which formation.18–20 Such type of issues do not occur in two-
v
areorbecomeamorphousunderlowenergy(E (cid:39)1keV) field formulations.9 Thus, the KS equation was consis-
i
X IBS, like semiconductors,4 Bradley and Harper (BH) pi- tently generalized into the so-called extended KS (eKS)
r oneered this approach through a linear continuum the- model for IBS.21,22 For normal incidence conditions,
a ory which explains the formation of ripples and their this model has been studied for one-dimensional (1D)
orientation,5 based on Sigmund’s theory of sputtering6 systems,23 and for 2D systems and rotating targets.24
and Mullins’ thermal surface diffusion.7 The success of Oblique incidence is studied in Ref. 25. While be-
this model to account for the origin of the patterns ing a phenomenological approximation of fuller hydro-
triggered an intense activity and further generalizations. dynamic descriptions,26 two-field modeling provides a
In particular, relevant nonlinear corrections were identi- generic framework which allows to modify the interface
fied in Ref. 8, leading to an equation of the Kuramoto- equationwhenimprovedmodelsoferosionand/ortrans-
Sivashinsky (KS) type.9 Importantly, nonlinearities were portareconsidered. Todate, thetwo-fieldmodeland/or
seen to moderate the pattern-forming linear instability the eKS equation have been (semi)quantitatively vali-
and eventually stabilize the surface morphology.
2
dated in several IBS experiments.27–31 In this paper we put forward a two-field model of IBS
nanopatterning under conditions in which anisotropies
The scenario just described focuses almost exclusively
to surface transport are relevant. As a basis for fur-
on targets for which the crystalline structure is over-
ther studies, our goal is to demonstrate non-trivial ef-
ruled by the IBS process. However, there are impor-
fects arising already within the simplest anisotropic sce-
tantinstancesinwhichthisisnotthecase,mostnotably
metals2,32 and semiconductors at high temperature.33–35 narios,whichwillmotivateourchoicesinthemodelingof
both, transport and irradiation-related mechanisms. As
In both cases, the strong dependence of the diffusivities
a result, we obtain a generalization of the eKS equation,
of ad-atoms and vacancies with the crystallographic di-
whichisintegratednumericallyfornormalionincidence.
rection can play a crucial role in the pattern formation
Our results show that anisotropic surface diffusion has
process. For metals, the surface is not amorphized af-
non-trivial effects, and allows to reproduce nanopatterns
ter ion impact. For semiconductors, increasing temper-
withdifferentlocalorderingstructuresinmonoelemental
atures above the recrystallization value analogously re-
systems, from hexagonal to square, akin to those experi-
stores dynamical dominance of crystalline anisotropies.
mentally reported for IBS of metals.30,42
AsdiscussedinRefs.32and36formetallicsystems, two
This paper is organized as follows. Our generalized
regimes can be distinguished: (i) diffusive regime, when
two-field model with anisotropic diffusion is put forward
pattern formation is governed by thermal surface diffu-
in Section II. In principle, the model holds for arbitrar-
sion, typically for intermediate temperatures and rela-
ily oblique ion incidence. However, in order to isolate
tively low ion fluxes, and (ii) erosive regime, when pat-
the effect of anisotropy in diffusion, rather than in irra-
tern formation is controlled by the direction of the ion
diation, we then restrict ourselves to normal incidence.
beam, usually for very high or very low temperatures
For this case we derive an equivalent interface equation
and for large enough ion fluxes. For instance, the diffu-
which generalizes the eKS model. This novel nonlinear
sive regime allows for anisotropic ripple formation under
equation is studied numerically in Section III, where the
isotropic,normalincidenceconditions,andingeneralim-
effect of each one of the parameters which control the
pliesthatboththeripplewavelengthandorientationare
controlled by temperature.32 This behavior can not be system behavior is discussed in detail. Finally, Section
IVcontainsourconclusionsandanoutlookonfuturede-
explained using previous models of IBS for amorphous
velopments. Some details on our modeling are provided
targets, in which the diffusive terms are isotropic. A
in the Appendix.
generalization of the linear BH model to anisotropic ma-
terials was proposed in Refs. 37 and 38. Some proper-
ties observed in IBS of metals could thus be described,
II. GENERALIZED TWO-FIELD MODEL
butinthisformulationsurfacetransportdoesnotcouple
with erosion in a natural way. Likewise, with a focus on
strongly kinetic effects, previous two-field10,11 and one- A. Derivation
field models33–35 have described crystalline anisotropies,
but only at nonlinear order. However, in principle under A two-field model is a system of two coupled partial
theseconditionssurfacediffusioncurrentsneedtoinclude differential equations describing the temporal evolution
anisotropic linear terms,39,40 which account for e.g. the of two important macroscopic variables.9 The first vari-
direction-dependence of barriers to adatom/advacancy ablecorrespondstotheheightofthebombardedsurface,
diffusion on terraces, along step edges, etc.41 h(x,t), at substrate position x=(x,y) and time t. The
second one describes the thickness (which, for a fixed
In view of the above, there is a need for studies in
atomicvolume,isproportionaltothedensity)ofthemo-
which crystalline anisotropies to material transport are
bile surface adatoms layer, R(x,t). For semiconductors
systematicallyaddressed,forsurfacesundergoinglowen-
at room temperature, irradiation creates an amorphous
ergy IBS. Already the simplest scenario of anisotropic
layer with a thickness of the order of the ion range,4
linear surface diffusion can lead to non-trivial modifica-
within which transport can be described in terms of vis-
tions of pattern properties, even if possibly not modi-
fying other, such as stability phase diagrams.39,40 For cous flow.26,43,44 However, for metals or for semiconduc-
instance, recent experiments with gold targets30,42 have tors at high temperatures, such amorphization does not
take place,32,36 so that the surface layer on which trans-
obtained highly ordered nanodot patterns by sequential
port occurs can be assumed to have roughly an atomic
ion-beam sputtering (SIBS). The procedure consisted in
thickness. In such cases, the dynamics of h and R are
sputteringundernormalincidenceapre-patternedripple
coupled by mass conservation as
structure previously obtained by oblique bombardment.
When the initial surface is flat and not pre-patterned,
∂R
a more disordered dot pattern is obtained, which still ∂t =(1−φ)Γex−Γad−∇·J, (1)
shows square in-plane order.42 Although the (isotropic)
∂h
eKSmodelreproducesmanyoftheexperimentalproper- =−Γ +Γ . (2)
∂t ex ad
tiesoftheensuingnanobeadpattern,30,42 itisnotableto
predictthissquaresymmetry,beinglimitedtodescribing Here, Γ is a function that describes the rate at which
ex
more isotropic, hexagonal order. targetatomsareexcavated(locallydecreasingh)andcan
3
become mobile (locally increasing R), while Γ models diffusion and IBS, at this stage ES-related mechanisms
ad
the rate of atom addition back to the solid (increasing h are left for further work. They are expected to play a
and decreasing R). The parameter φ ∈ [0,1] measures significant role in patterns whose in-plane order extends
the fraction of eroded atoms that are actually sputtered to a longer range, and in which wavelength coarsening
away from the surface, while φ¯ = 1 − φ measures the is more sizeable, than in e.g. experiments on nanobead
fractionoferodedatomsthatremainsubjecttotransport formation.30,42
at the surface. Inordertoclosethesystemofequations(1)-(2),theex-
Equation(1)includesanadditionalconservedcurrent, cavationandadditionrateshavetoberelatedtotheden-
J,whichaccountsforsurfacetransportmechanisms. For sityofadatoms(R)andtothegeometryofthesubstrate
instance, this current could readily incorporate Carter- (handitsspacederivatives)themselves. Foranarbitrary
Vishnyakov (CV) contributions45 due to mass redistri- incidenceangleθandassumingthattheprojectionofthe
bution, believed to be relevant in the case of semicon- ion beam is along the xˆ direction, we consider21,22
ductors at room temperature.46 These have been em-
(cid:104) ∂h ∂2h ∂2h
ployed in a number of similar two-field models for IBS
Γ =α 1+α +α +α (5)
of compound systems, or for IBS of monoelemental tar- ex 0 1x∂x 2x∂x2 2y∂y2
getsunderconcurrentimpuritycodeposition,seee.g.Ref. (cid:16)∂h(cid:17)2 (cid:16)∂h(cid:17)2(cid:105)
+α +α ,
47 for a partial overview. For monoelemental targets, 3x ∂x 3y ∂y
CV-type effects can also be reflected in Γ .21 In any
ex
case, for metals or semiconductors at high temperatures where α0 > 0 is the excavation rate for a flat surface.
thistypeofmassredistributionis not expected toplaya In Eq. (5) the terms with coefficients α1x, α2j corre-
role, nor is, on a more mesoscopic level, surface-confined spond to the lowest linear-order approximation to the
viscous flow.26,43,44,48 The main transport mechanism is dependence of the sputtering yield on the local height
expected to be, rather, thermal surface diffusion. Mi- derivatives, as in BH’s theory,5 while those with coef-
croscopically, this is an activated process for which ener- ficients α3j characterize the corresponding lowest-order
getic barriers exist, whose depths depend on the crystal- nonlinear corrections.8 Due to the assumed geometry for
lographicdirections.49 Onamorecoarse-grainedlevel,as ion bombardment, for normal incidence (θ = 0) one has
e.g. in Mullins’ classic theory,7 surface diffusion is medi- α1x = 0, while α2x = α2y, α3x = α3y. In general, in the
ated by surface tension, which for metals is paradigmat- absence of CV-type effects, one expects α2j < 0 leading
ically anisotropic.50 Mathematically, we thus consider topatternformation(BHinstability),whilenon-zeroα3j
Fickian diffusion at the surface as described by guarantee non-exponential increase of the surface rough-
ness for long times, as mentioned above.
J =−D∇R, (3) Finally, for the local addition rate we consider22
where D ∈R2×2 is a (positive definite) diffusion tensor, (cid:104) (cid:16) ∂2h ∂2h(cid:17) (cid:105)
Γ =γ R 1+γ +γ −R , (6)
ratherthanaconstant,thatimplementsthepresenttype ad 0 2x∂x2 2y∂y2 eq
of anisotropy. Its most general form reads
where γ >0 is the nucleation rate, i.e., 1/γ represents
0 0
(cid:20)D 0 (cid:21) (cid:20)D D (cid:21) the average time in which ad-atoms incorporate to a flat
D =M(ψ) 0(cid:107) D M−1(ψ)= Dxx Dxy , (4) surface. Asdiscussedearlier,24,25 inabsenceof ion-beam
⊥ xy yy
driving Eq. (6) describes Mullins’ thermal surface diffu-
where M(ψ) is a counterclockwise rotation matrix of sion, in such a way that Req is related with the surface
angle ψ which gives the orientation of the fast diffu- concentration of mobile species, while γ2j ≥ 0 are sur-
sion direction with respect to the xˆ direction, so that face tension coefficients which, in general, can also be
D ≥D >0 without loss of generality. anisotropic.
(cid:107) ⊥
Similar models of anisotropic surface diffusion have The two-field model (1)-(6) supports a flat solution in
beenemployede.g.inmodelsofthedynamicsofvicinal51 which the surface height erodes with an uniform speed
and singular52 surfaces in epitaxy, which have been ex- and both h and R are space-independent functions.55
perimentally validated.53 Note that the present (surface Performing a standard linear stability analysis of per-
diffusion) anisotropy is independent of that induced by turbations of this solution which are periodic with wave-
ion bombardment under an oblique angle of incidence; vectork,wecanobtainthepatternwavelength,(cid:96)i,along
consequently,itisstillrelevantunderotherwiseisotropic each direction i = x,y within linear approximation.
normal incidence conditions. In principle, J could incor- Specifically, we define (cid:96)i = 2π/ki(cid:96) as the length-scales
porateadditionalterms,mostnotably(linearandnonlin- at which the dispersion relation is maximized in each
ear)contributionsdependingonthesurfaceheight,which direction. In our case,55 k(cid:96) = ((cid:15)φγ α /2R D γ )1/2,
i 0 2i eq i 2i
arerelatedwithEhrlich-Schwoebel(ES)anisotropicbar- where the parameter (cid:15)≡α /γ turns out to be small as
0 0
riers to surface diffusion.33–35,41,54 However, the mor- a consequence of the difference between the typical time
phological instability associated with these terms differs scales associated with diffusion (typically of the order of
physically from the BH instability. In order to assess ps) and the ion-beam driving (of the order of seconds
more clearly the interplay between anisotropic surface for the ion fluxes usually employed).25 This separation
4
in time scales allows one to simplify the analysis of the where K and Λ are matrices defined as
2
mathematical model [Eqs. (1)-(2)] since it allows to per-
form a multiple-scale perturbative analysis to obtain a (cid:20)K 0 (cid:21) (cid:34)λ(2) 0 (cid:35)
K= x and Λ = x . (9)
closed equation for the height. Analysis that closely fol- 0 Ky 2 0 λ(y2)
lows Ref. 25 leads to an effective nonlinear equation for
the time evolution of h, which reads
The number of independent parameters in Eqs. (8) and
(9) has reduced dramatically, the remaining ones be-
∂h ∂h (cid:88) ∂2 (cid:16)∂h(cid:17)
= γ + Ω ing ν = (cid:15)φγ α , λ(1) = −(cid:15)φγ α , K = D R γ +
∂t x∂x ij∂i∂j ∂x 0 2 0 3 i i eq 2
i=x,y (cid:15)(φR γ γ −φ¯D )α , and λ(2) =(cid:15)(φR γ γ −φ¯D )α ,
eq 0 2 i 2 i eq 0 2 i 3
+ (cid:88) (cid:20)−ν ∂2h +λ(1)(cid:16)∂h(cid:17)2(cid:21) (7) wherei=x,y. Itisimportanttonotethat,incontrastto
i∂i2 i ∂i the equation obtained in Ref. 22,24, in which only terms
i=x,y oftheform∇2hand(∇h)2 appearundernormalioninci-
− (cid:88) (cid:20)K ∂2 (cid:16)∂2h(cid:17)+λ(2) ∂2 (cid:16)∂h(cid:17)2(cid:21), dence, in the case of Eq. (8) the second-order derivatives
ijk∂i∂j ∂k2 ij ∂i∂j ∂k ∂2/∂x2 and ∂2/∂y2 are weighted by parameters that de-
i,j,k=x,y pend on the different diffusion coefficients.
Under experimental conditions leading to pattern for-
where the coefficients γ , Ω , ν , λ(j), and K depend
x ij i i ijk mation,ν >0inEq.(8). Withrespecttothecoefficients
on ion energy, flux, incidence angle, etc. through their
of the linear fourth-order derivative term in this equa-
dependencies on α and all other parameters entering
ij tion, note that they contain contributions that couple
Γ and Γ , as specified in Appendix A.
ex ad different physical mechanisms in a natural way. Thus,
Equation(7)ispartiallysimilartotheevolutionequa-
the contribution proportional to surface diffusivity and
tionobtainedinRef.22forisotropicsurfacediffusionand
surface tension is completely analogous to the form of
oblique ion incidence. However, in that case the only ge-
Mullins’ surface diffusion, although note that the coeffi-
ometrical condition responsible for breaking the x ↔ y
cients may also include ion-induced contributions which
symmetry was the non-zero value of the incidence an- are temperature-independent.25 The remaining term in
gle, in such a way that the system was symmetric under
K coupleserosion(beingproportionaltoα )withtrans-
i 2
space reflection in the y direction, but not in the x di-
port(surfacediffusivity)andsurfacetension, furtherim-
rection. In the case of Eq. (7), this same cause for space
plementingion-induceddiffusivity. Wewillconsidercon-
anisotropy is enhanced by anisotropic surface diffusion
ditions in which this fourth-order derivative term has a
and by anisotropic surface tension. As a consequence,
net smoothing effect, so that K > 0. Finally, we con-
i
not only are the x ↔ y and x ↔ −x symmetries bro-
sider the products λ(1)λ(2) to also be positive. Mathe-
ken,butthey ↔−y symmetryisbrokenaswell,nowby j
matically, thisconditionisrequiredforEq.(8)tobefree
the two latter conditions. The differences between Eq.
of so-called “cancellation modes”, known to occur un-
(7) and the one obtained in Ref. 22 will be further dis-
der appropriate conditions in related continuum models,
cussed in the following sections, for the case of normal
such as the anisotropic KS56 and eKS18,22 equations.
ion incidence.
WenextrescaleEq.(8)inordertoworkindimension-
less units. This allows to perform generic statements on
the system behavior, while at the same time it also sim-
B. Effective Equation for Normal Incidence
plifies the discussion by minimizing the number of free
parameters. Hence, we define
Having as a reference experimental behaviors those
reported in Refs. 42 and 30, in which an initial Au- (cid:16)K(cid:17)1 (cid:16)K(cid:17)1 K ν
prepatterned surface was further irradiated at normal x= 2x(cid:48), y = 2y(cid:48), t= t(cid:48), h= h(cid:48),
ν ν ν2 λ(1)
incidence, we will focus here in such condition θ = 0.
This implies24,25 α1x = 0, α2x = α2y = α2, and where K = (Kx+Ky)/2. Dropping the primes, Eq. (8)
α3x = α3y = α3, and will allow us to isolate the effects now reads
purely due to anisotropies in surface diffusion. For this
reason, we will moreover assume isotropic surface ten- ∂h = −∇2h+(∇h)2−2∇·(cid:2)A∇(cid:0)∇2h(cid:1)(cid:3) (10)
sion, namely, γ =γ =γ . As in Refs. 42 and 30, we ∂t
2x 2y 2
will also take x and y to be aligned with the substrate −2r ∇·(cid:8)B∇(cid:2)(∇h)2(cid:3)(cid:9) ,
0
directionsalongwhichsurfacediffusivitiesareoptimized.
Under these conditions, both the excavation and the ad- where the matrices A and B are
dition rates become isotropic, Eq. (7) taking the simpler
(cid:20) (cid:21)
form A= K = αx 0 ,
Kx+Ky 0 1−αx
∂h
= −ν∇2h+λ(1)(∇h)2 (8) Λ (cid:20)β 0 (cid:21)
∂t B = 2 = x ,
−∇·(cid:2)K∇(cid:0)∇2h(cid:1)(cid:3)−∇·(cid:8)Λ2∇(cid:2)(∇h)2(cid:3)(cid:9), λ(x2)+λ(y2) 0 1−βx
5
with code was implemented in MATLAB, being based on a
standard finite-difference scheme in space (for the linear
α = Kx , β = λ(x2) , r = ν (cid:16)λ(x2)+λ(y2)(cid:17). terms) and a fourth-order Runge-Kutta method for the
x Kx+Ky x λ(x2)+λ(y2) 0 λ(1) Kx+Ky time evolution, using a spatial grid with 256×256 nodes,
atimestep∆t=0.01,andaspacestep∆x=1. Thedis-
Given that three independent rescalings have been per- cretization of the nonlinear terms was based on the one
formed on Eq. (8), which depends on six independent proposedbyLamandShininRef.58. Wehaveemployed
parameters, thefinalEq.(10)dependsonthreeindepen- periodic boundary conditions and initial height values
dent constants only, αx, βx, and r0. Note that, in all the which are uniformly distributed between 0 and 0.1. Be-
physically relevant cases, r0 >0. sides inspection of the resulting surface morphologies, in
It is interesting to stress some of the features of Eq. allcaseswehavecalculatedtheglobalsurfaceroughness,
(10): (i)Asexpected,theanisotropiesareonlycausedby W, as well as the wavelengths, (cid:96) and (cid:96) , after averag-
x y
the different diffusivities Dx and Dy. This is reflected in ing over 10 realizations of the initial condition for each
thefactthattheparametersαxandβxwillgenerallytake parameter set. Additionally, we have computed the nor-
valuesdifferentfrom1/2. Therefore,theweightsforboth malizedautocorrelationfunction, R , thatallowsoneto
N
directions in the last two terms of Eq. (10) will be differ- determine the local arrangement of the patterns and is
ent. When αx = βx = 1/2 the isotropic diffusion equa- defined as59
tion for normal incidence proposed in Ref. 21 is recov- 1 1 (cid:90)
ered. (ii) The dimensionless parameter r is the squared R (x,t)= [h(x+r,t)h(r,t)−h¯2(t)]dr, (12)
0 N W L2
ratio of two length scales. One of these length scales
is computed as the ratio between the parameters of the where W is the surface roughness and h¯ is the mean
conservedKardar-Parisi-Zhang(KPZ)22 nonlinearterms height over the whole spatial grid of size L×L.
appearing in the equation, λ(2), and that of the non-
i
conserved KPZ nonlinearity that ensues, λ(1), namely,
[(λ(2)+λ(2)]/λ(1))1/2. The second length scale is set by A. Isotropic case: αx =βx =1/2
x y
the parameters of the linear terms, as [(K +K )/ν]1/2.
x y
Theparametercombinationr thusprovidesanestimate To begin with our analysis, and for the sake of later
0
oftherelativerelevanceofthevariouscontributionsthat comparison with anisotropic parameter conditions, we
compete in the dynamics of the system. first recall the results obtained in Ref. 24 for isotropic
Asjuststressed,ingeneralEq.(10)isanisotropic,thus systems under normal ion incidence, as a special case of
thepatternswillpresentdifferentwavelengthsalongeach our model in which both surface diffusivities are equal,
principal direction. These can be estimated as functions D =Dx =Dy. In this case Eq. (8) simply reduces to
of the parameters of the equation, by performing a stan-
∂h
dard linear stability analysis.25 This leads to =−ν∇2h−K∇4h+λ(1)(∇h)2−λ(2)∇2(∇h)2, (13)
∂t
(cid:114) (cid:115)
(cid:96) = 2π =2π 2Ki ≈2π 2ReqDiγ2i, (11) where K = Kx = Ky and λ(2) = λ(x2) = λ(y2). After
i k(cid:96) ν (cid:15)φγ α a rescaling which is similar to the one employed in the
i 0 2i
previous section, this equation reduces to the particular
where we have substituted the values of Ki and ν pro- case of Eq. (10) in which αx = βx = 1/2. Note that, in
vided after Eq. (9). As expected, these wavelengths co- principle, the simulations reported in Ref. 24 correspond
incide with the values obtained in Section II for the lin- to the unrescaled Eq. (13).
ear stability analysis of the full two-field model. Impor- The following features for the surface roughness and
tantly, the expressions obtained for (cid:96)i can be often used pattern wavelength were obtained in this case:24 (i) The
to perform (semi)quantitative comparisons between the surface morphology shows a short-time transient behav-
present type of continuum models and experiments at ior. During this interval W grows exponentially with
shorttimes,priortotheonsetofnon-lineareffects.28,30,57 time and a dot pattern appears whose characteristic
wavelengthisaccuratelydescribedbythelinearanalysis.
Indeed, this stage is controlled by the linear terms ν∇2h
III. RESULTS andK∇4h. (ii)Afterthislinearregime,acrossovertakes
place towards a behavior which is controlled by the con-
Thus far, we have been able to derive an effective di- served nonlinear term λ(2)∇2(∇h)2, in which the growth
mensionless equation for normal ion incidence that con- of W and (cid:96) in time can be approximated by power-laws,
tains all the physical mechanisms of the problem and with effective exponents whose values depend on equa-
depends on three free parameters only, Eq. (10). In this tion parameters. (iii) For long times, the non-conserved
sectionwestudysystematicallythisequationbyindepen- nonlinear term λ(1)(∇h)2 induces eventual saturation of
dently changing the values of each of these three param- W and (cid:96), and height disorder at large scales.
eters. Given the strong nonlinearities in the equation, Actually, the relative duration of the various dynami-
we resort to a numerical integration. Specifically, our cal regimes turns out to be controlled by the parameter
6
r .22,24 Thus, large r values correspond to the predom- obtained by measuring the distance from the origin to
0 0
inance of the conserved KPZ nonlinearity at intermedi- the first maximum of the autocorrelation function along
ate times, allowing for a stronger coarsening process and the corresponding axis. For all the values of α consid-
x
an improved order of the height values throughout the ered, we obtain (cid:96) (cid:39) (cid:96) , the only significant difference
x y
surface, namely, a smaller roughness. On the contrary, being that, for α = 0.25 and α = 0.3, the first peak
x x
small r values correspond to a non-linear regime dom- of the autocorrelation function along the y-direction is
0
inated by the KPZ nonlinearity, with a relatively short higher than the first peak along the x-direction, both
intermediate coarsening regime and a rougher surface at peakshavingthesameheightsforα =0.5. Thisimplies
x
long times. that,forα ∈(0,0.5),dotsaremorecorrelatedalongthe
x
Atthispoint,itisimportanttoremarkthatthe1Dand y-direction, although the differences are not substantial.
2DbehaviorsofEq.(13)differquitestronglywithrespect Figure2bshowsthetimeevolutionoftheglobalrough-
to the ordering properties. To begin with, note that, nessW andwavelengths,(cid:96) (solidsymbols)and(cid:96) (open
x y
given the fact that the band of unstable modes extends symbols), for the same values of α as in Fig. 2a. As for
x
downtok=0,combinedwiththeoccurrenceofthenon- the isotropic eKS model, three time regimes can be dis-
conserved KPZ nonlinearity, order in the pattern can be tinguished: Initially the roughness grows exponentially,
short-range at most.60 This does not prevent the equa- up to intermediate times when its growth rate slows
tionfromprovidingaquantitativelyaccuratedescription down; eventually it reaches a similar time-independent
of experimental patterns.28,57 Then, while relative ho- value for all α . On the other hand, the behavior of the
x
mogeneity in height values correlates positively with an pattern wavelengths with α is different. Initially both
x
enhanced“in-plane”orderforEq.(13)in1D,22thisisnot (cid:96) and (cid:96) start growing slowly, with (cid:96) being larger than
x y y
thecasein2D.Namely,forthe2Dcase,smallerr values (cid:96) for α < 0.5. This is due to the fact that the fourth-
0 x x
seemtofeatureimprovedin-planeshort-rangehexagonal order linear terms controlled by the parameter α are
x
ordering of the dot structure. Conversely, larger values expectedtoplayanimportantrolepreciselyatthesmall
ofr ,thatallowforstrongercoarsening(widercells)and spatialand time scalesat which the linear instabilityde-
0
smaller overall roughness, correspond to surfaces with velops. As a matter of fact, looking at the expressions of
poorer in-plane ordering. For a thus-far unreported ex- thewavelengthspredictedbythelinearinstabilityanaly-
plicitcomparison,seeFig.1,inwhichther =10and50 sis,Eq.(11),wecaneasilynotethat(cid:96) shouldbesmaller
0 x
casesareexplicitlyillustrated. TheresultsshowninFig- than(cid:96) ifK <K ,whichisindeedthecaseforα <0.5.
y x y x
ure 1 also show that in both cases the wavelengths along This linear transient behavior is followed by a coarsen-
each direction, (cid:96) and (cid:96) , are equal, since the surface ing process controlled by the nonlinear terms, which fi-
x y
diffusion is isotropic. nallydrivebothwavelengthstosimilarsaturationvalues.
The rich crossover behavior of Eq. (13) will be use- Thus,sincetheparametersofthenonlineartermsr and
0
ful to better understand the morphologies described by β arefixed,thefinalsurfacetopographiesareverysimi-
x
Eq. (10) as a function of parameters values. laratlongtimesforthedifferentvaluesofα . Addition-
x
ally,becausethenonlineartermsareisotropic(β =0.5)
x
the wavelengths reach similar values in both directions
B. General values of α at long times. In the next sections we study the impact
x
of the coefficients of the nonlinear terms on the system
dynamics and pattern formation and evolution.
To study the effect of non-isotropic values of the pa-
rameter α mediating the linear surface-diffusion terms
x
in Eq. (10), we have integrated numerically Eq. (10) for
r = 10, β = 0.5, and values of α ∈ [0,1]. Due C. General values of βx
0 x x 2
to the symmetry of Eq. (10) with respect to reflections
of α around the isotropic 1/2 value, the behavior for We next consider the influence on the topography
x
α ∈ [1,1] can be easily obtained from our simulations of the anisotropic, conserved nonlinearity which is con-
x 2
by simply swapping the x and y axes. trolled in Eq. (10) by the parameter β . To this end,
x
Figure2ashowsthesurfacemorphologies,thenormal- numerical integrations of Eq. (10) have been performed
ized autocorrelation function, and cross-cuts of the nor- for fixed values of αx and r0, and βx ∈ [0,21]. Similarly
malized autocorrelation function along the x and y axes, to the case of α , results for β ∈ [1,1] can be deduced
x x 2
RN and RN, respectively, at t = 1000 and for different fromthesimulationsshownnextbyswappingthexandy
x y
values of α . Note that, as discussed above, α = 0.5 axes. Figure 3a displays the morphology, the normalized
x x
corresponds to the isotropic case already studied in Ref. autocorrelation function, and the autocorrelation func-
24. Asnoticedbyinspectingtheleftcolumnofthefigure, tion along the x and y axes at t = 1000 for different
thesurfacemorphologydoesnotchangequalitativelyfor valuesofβ . Theadditionalcaseβ =0.5forthechosen
x x
different values of α . This robustness with respect to α correspondstotheisotropicsystemalreadyshownon
x x
α is further evidenced in the middle and right panels the third row of Fig. 2a. The temporal evolution of the
x
of Fig. 2a, in which RN is shown side-by-side with RN surface roughness and wavelengths are shown in figure
x
and RN. The x and y wavelengths, (cid:96) and (cid:96) , can be 3b.
y x y
7
α =0.5,β =0.5,
x x
0.4 r =10
0
r =50
0.3 0
W0.2
0.1
0
100 101 102 103
t
40 r0=10: (cid:2)x
r =10: (cid:2)
0 y
30 r =50: (cid:2)
0 x
(cid:2),xy20 r0=50: (cid:2)y
(cid:2)
10
0
100 101 102 103
t
(a) (b)
FIG. 1: (Color online) (a) Top-view surface morphologies (left column), normalized autocorrelation functions, R
N
(center column), and normalized autocorrelation functions along the x and y directions, RN and RN respectively
x y
(right column), predicted by Eq. (10) at t=1000 for α =0.5, β =0.5, and different values of r (see legends). (b)
x x 0
Temporal evolution of the roughness, W, and the wavelengths along the x and y directions, (cid:96) (solid symbols) and
x
(cid:96) (open symbols) respectively, for the same parameter values as in (a).
y
β =0.5,r =10,
x 0
0.4
α =0.25
x
α =0.30
0.3 αx=0.50
x
W0.2
0.1
100 0 101 102 103
t
20
α =0.25: (cid:2)
(cid:2)(cid:2),xy1105 αααxxx===000...233500::: (cid:2)(cid:2)(cid:2)xyx
x y
α =0.50: (cid:2)
x x
α =0.50: (cid:2)
5 x y
10 0 101 102 103
t
(a) (b)
FIG. 2: (Color online) (a) Top-view surface morphologies (left column), normalized autocorrelation functions, R
N
(center column), and normalized autocorrelation functions along the x and y directions, RN and RN respectively
x y
(right column), predicted by Eq. (10) at t=1000 for β =0.5, r =10, and different values of α (see legends). (b)
x 0 x
Temporal evolution of the roughness, W, and the wavelengths along the x and y directions, (cid:96) (solid symbols) and
x
(cid:96) (open symbols) respectively, for the same parameter values as in (a).
y
8
Ifβ issmall,theconservednonlinearityactspredomi- Thesimulationresultsfordifferentvaluesofr areshown
x 0
nantlyalongthey-axis,inducingstrongercoarseningbe- in Fig. 4a. Analogously to the isotropic case for normal
havior,hence(cid:96) becomeslargerthan(cid:96) fortimesbeyond incidencestudiedinSectionIIIAandillustratedinFigs.
y x
the linear regime, see Fig. 3b. Actually, this behavior is 1a and 1b, in the presence of anisotropic surface diffu-
associated with a change in the pattern symmetry, see sion the patterns present more coarsening and a smaller
e.g.theβ =0.2caseinFig.3a. Indeed, thelargervalue roughness when r is larger. However, the quality of in-
x 0
of(cid:96) impliesthatthedotsorcellsbecomemoreelongated plane ordering of the dots is poorer. For the parameter
y
in the y-direction, leading to the emergence of a ripple valuesconsideredinFig.4a,slightlyelongateddotsgroup
patternwithridgesparalleltoit. Thiscanbenotedboth together following square arrangements, as can be noted
in the morphology and in the autocorrelation functions, lookingatthesurfacemorphologies. However,theshort-
and is in spite of the fact that we are considering normal range square order is hindered for larger r values. This
0
incidence conditions for the ions. Note, this is a purely is also reflected in the autocorrelation function, where a
non-linear effect, as Eq. (10) is completely isotropic at moreperfectsquarepatternisrevealedforsmallervalues
linear order for this parameter condition. On the other of r .
0
hand,ifβx increases,theelongationofthedotsalongthe The temporal evolution of the roughness and wave-
y-direction is attenuated, they form arrangements with lengths for different values of r are represented in Fig.
0
a more square (rather than rectangular) symmetry, and 4b. Again three main regimes can be distinguished. The
the effect is mitigated, see Fig. 3a for βx = 0.3. At any roughness grows exponentially in the first, linear regime,
rate, for βx (cid:54)= 1/2 the isotropy of the pattern is clearly followed by power-law growth, and by saturation at very
broken. long times. As in Refs. 22 and 24, the final roughness is
With respect to the time evolution of the surface indeed smaller for large r values, the long-time config-
0
roughness,Fig.3bindicatesanunambiguousdependence urations showing more uniform height values. For such
withthevalueofβ ,whichcontrastswiththeresultsob- large r , the two wavelengths (cid:96) and (cid:96) are also larger,
x 0 x y
tained for α in the preceding section. For small values due to the longer coarsening process undergone. Note,
x
of β , saturation occurs later and the saturation value is because α = 0.5, the linear terms have the same effect
x x
larger. This seems reminiscent of results for the 1D eKS in both directions. Since β = 0.3 in the simulations
x
equation when increasing the strength of the conserved shown, and as we saw in the previous section, the wave-
KPZnonlinearitywithrespecttotheremainingtermsin length grows more in the y direction and patterns with
the equation.22 On the other hand, for short times the (cid:96) >(cid:96) are always obtained.
y x
roughness values are practically the same for all β , sug-
x
gesting that such an increase of the roughness is indeed
a nonlinear effect. Regarding the pattern wavelengths
E. Unusual patterns and order under normal
in the two directions, both grow very slowly and take
incidence: Ripples and square or hexagonal dot
similar values during the short times associated with the arrangements
linear instability. At intermediate times, both grow at
increased rates; ultimately, they reach very different sat-
AssuggestedbyFig.4a,forintermediatevaluesofβ it
uration values depending on the specific value of β . In- x
x ispossibletogeneratesurfacesforwhichthedotpatterns
deed,asalreadynotedabove,forrelativelysmallvaluesof
display short-range order with square symmetry. Such a
thisparameterthepatternwavelengthinthey-direction,
morphology is locally characterized by each single dot
(cid:96) , becomes larger than (cid:96) , as can be clearly appreciated
y x having on average four nearest neighbors located along
in Fig. 3b already for β = 0.3. For even smaller values
x the two Cartesian directions. Here we employ the nor-
of β , such as β = 0.2, (cid:96) interrupts its growth process
x x x malizedheightautocorrelationfunction,R ,toquantify
early while (cid:96) keeps growing for a long time (note that N
y the spatial order on the surface. If the morphology does
itscoarseningprocesshasnotyetstoppedatt=1000for
correspond to a pattern with such a square-symmetric
β = 0.2), resulting into very different (cid:96) > (cid:96) . This is
x y x order, the central maximum of R lies within a perfect
due to the fact that the conserved nonlinear term, which N
squareformedbyeightnearestsatellitepeaksonasquare
inducesthecoarseningprocess,isstrongeralongtheydi-
arrangement. Although the surfaces described by Eq.
rection. Insummary,theroleofβ istwofold: itmodifies
x (10) present this type of structure to a certain degree,
the local arrangement (symmetry and order) of the pat-
it is not possible to obtain a strictly square symmetry
terns and it amplifies/reduces the coarsening dynamics
due to the anisotropy introduced by the conserved non-
selectively along one of the system directions.
linearities controlled by the parameter β . Indeed, the
x
heterogeneoussurfacediffusivitiesinthetwospacedirec-
tions lead to different wavelengths, even under isotropic
D. General values of r0 (normalincidence)irradiation,producingarelativelyor-
dered array of dots, but with different typical sizes in
We continue in this section with the morphological each direction. For β ∈(0,0.5), dots are more correlated
effects of the third independent parameter in Eq. (10), withtheirneighborsalongthexdirectionthanalongthe
namely, the ratio of nonlinear to linear length scales, r . y direction, as can be noticed in the height autocorrela-
0
9
α =0.5,r =10,
x 0
0.8
β =0.2
x
0.6 βx=0.3
β =0.5
x
W0.4
0.2
100 0 101 102 103
t
100 βx=0.2: (cid:2)x
β =0.2: (cid:2)
x y
80 β =0.3: (cid:2)
x x
(cid:2)y 60 ββx==00..35:: (cid:2)(cid:2)y
(cid:2),x 40 βxx=0.5: (cid:2)xy
20
100 0 101 102 103
t
(a) (b)
FIG. 3: (Color online) (a) Top-view surface morphologies (left column), normalized autocorrelation functions, R
N
(center column), and normalized autocorrelation functions along the x and y directions, RN and RN respectively
x y
(right column), predicted by Eq. (10) at t=1000 for α =0.5, r =10, and different values of β (see legends). (b)
x 0 x
Temporal evolution of the roughness, W, and the wavelengths along the x and y directions, (cid:96) (solid symbols) and
x
(cid:96) (open symbols) respectively, for the same parameter values as in (a).
y
α =0.5,β =0.3,
x x
0.5
r =10
0
0.4 r =20
0
r =50
0.3 0
W
0.2
0.1
100 0 101 102 103
t
r =10: (cid:2)
0 x
50 r =10: (cid:2)
0 y
40 r0=20: (cid:2)x
r =20: (cid:2)
0 y
(cid:2)y30 r0=50: (cid:2)x
(cid:2),x20 r0=50: (cid:2)y
10
100 0 101 102 103
t
(a) (b)
FIG. 4: (Color online) (a) Top-view surface morphologies (left column), normalized autocorrelation functions, R
N
(center column), and normalized autocorrelation functions along the x and y directions, RN and RN respectively
x y
(right column), predicted by Eq. (10) at t=1000 for α =0.5, β =0.3, and different values of r (see legends). (b)
x x 0
Temporal evolution of the roughness, W, and the wavelengths along the x and y directions, (cid:96) (solid symbols) and
x
(cid:96) (open symbols) respectively, for the same parameter values as in (a).
y
10
tion functions obtained in the previous sections. See for with(cid:96) >(cid:96) occur. Moreover,fortheseparametervalues
y x
example Fig. 4a, where the correlation values are clearly the surface heights becomes more heterogeneous, differ-
larger along the x-axis. Recall that decreasing β in this entlocalregionspresentingdifferentaverageheights. For
x
range of values actually increases the elongation of dots evenlowervaluesofβ ,aripplestructureappears,witha
x
along the y axis. periodicityalongthey-direction. Againthismorphology
Enhancement of local square order can be achieved displays quite heterogeneous average heights in different
bringing together the previous property with the fact regions, while it still features a short-scale structure of
thatlocalorderisimprovedforrelativelysmallr values. rather elongated dots which are quite ordered along the
0
Thus, Fig. 5 displays the surface morphologies obtained x-direction. Hence, decreasing the value of β induces
x
for r = 5 and different values of β and their corre- a transition from short-range hexagonal, to square and
0 x
sponding autocorrelation maps. The symmetry of the then to rectangular ordering of the dots.
short-range order of the pattern can be identified easily Ananalogoustransitionbetweenhexagonalandsquare
intheautocorrelationmap,whichhasbeencalculatedfor patterns has been studied in Ref. 61 for the case of mag-
the (100×100) black boxes indicated. Indeed, since the netic fluids under applied magnetic fields. In this work
morphology is disordered at long distances, some of the the authors employ an angular correlation function that
localorderinformationislostwhentheheightautocorre- makes use of the discrete Fourier transform of theheight
lationfunctioniscomputedinthewholedomain. Atany field in order to characterize the (hexagonal or square)
rate,Fig.5showshowdotswithsquare-symmetricshort- symmetry of the pattern, and thus assess morpholog-
range order can actually occur for intermediate values of ical transitions under changes in external parameters.
β when the elongation along the y direction exists but Here, we define a similar angular correlation function,
x
is not excessively pronounced. but relative to the values of R , rather than those of
N
h(x,t). Specifically, the angular autocorrelation func-
βx =0.1 βx =0.25 βx =0.5 tion, P(ψ,t), which we propose to quantify the pattern
order is
1 (cid:90)
P(ψ,t)= [R (x,t)R (M(ψ)x,t)−R¯2 (t)]dx,
L2 N N N
w
(14)
where R¯ (t) is the space average of the autocorrelation
N
function, R (x,t), over a square window of lateral size
N
L , andasaboveM(ψ)is thecounterclockwise rotation
w
matrix of angle ψ. The function P(ψ,t) thus measures
theheightcorrelationateverypositionintheconsidered
domain and compares it with the result obtained at a
position which is rotated by an angle ψ. We further de-
fine the normalized angular autocorrelation function as
P (ψ,t) = P(ψ,t)/P(0◦,t). The reason for considering
N
an area of lateral size L < L is because some rotated
w
FIG. 5: Top-view surface morphologies (top row) points M(ψ)x could remain out of the considered do-
predicted by Eq. (10) at t=1000 for αx =0.5, r0 =5, main for a square grid. Besides, due to the global disor-
and different values of βx (see legends). Corresponding derofthepatternsinducedbytheKPZnon-linearity,the
normalized autocorrelation functions (bottom row) short-range order of the pattern needs to be quantified
computed over the indicated squares of size 100×100. in smaller areas.
The normalized angular autocorrelation functions cor-
Closer inspection of Fig. 5 actually suggests that up responding to the three basic morphologies shown in
to three main types of patterns can be expected for Eq. Fig. 5 are displayed in Fig. 6. The local extrema in
(10), depending on the value of β : Ripples (with a dot- P(ψ,t) signal how well correlated are the points in the
x
ted substructure) and dots with square or with hexag- morphology with those rotated by an angle ψ. In all
onal short-range order. Moreover, as we have already cases, due to the system symmetry under a 2D space in-
seen, the degree of local order of the pattern can be en- version (x,y) → (−x,−y), this function is periodic with
hancedbytuningthevalueofr . Indeed,thethreemain 180◦ period. For β =0.1 (dotted red line) ripples form
0 x
different patterns just mentioned can be clearly distin- and the only rotation that leaves the system unchanged
guished in Fig. 5, where r = 5 has been fixed. In is precisely one with ψ = 180◦, hence the maxima in
0
the case of isotropic surface diffusion (β = 0.5), dots P(ψ,t) as a function of ψ are separated by this value.
x
with (cid:96) = (cid:96) group into hexagonal short-range order, A different behavior is found for β = 0.25 (solid blue
x y x
where each dot tends to be in the center of a hexagon line), when a dot pattern with square-symmetric order
formed by the nearest neighbor dots, and local regions appears. In this case the distance between consecutive
tend to have the same average height. For intermedi- maxima of P(ψ,t) is 90◦, since the height correlation is
atevaluesβ ∈[0.25, 0.3],square-orderedelongateddots itselfmaximizedafterarotationof90◦ forapatternwith
x