Modeling the Decay of Nanopatterns: A ... - American Chemical Society

4606

J. Phys. Chem. C 2007, 111, 4606-4613

Modeling the Decay of Nanopatterns: A Comparative Study between a Continuum Description and a Discrete Monte Carlo Approach Marcos F. Castez* and Ezequiel V. Albano Instituto de InVestigaciones Fisicoquı´micas Teo´ ricas y Aplicadas (INIFTA), UNLP, CONICET, Sucursal 4, Casilla de Correo 16, (1900) La Plata, Argentina ReceiVed: NoVember 30, 2006; In Final Form: February 7, 2007

The surface-diffusion-driven decay of 2-dimensional periodic nanopatterns is studied by means of a kinetic Monte Carlo model. Activation energy barriers are computed by using a harmonic approach. Accounting for the proper rates of all processes involved allows us to implement a real time algorithm. Thus, the relationship between real time and Monte Carlo time is discussed. By using this discrete approach, we recover the most relevant results expected from the linear theory of surface diffusion for initial surfaces within the small-slope approximation. We also focus on the decay kinetics of periodic nanopatterns (sinusoidal and rectangular) starting from nonequilibrium states far from the small-slope approximation. In this case, we found a nonexponential decay as well as the emergence of several interesting nanostructures (overhangs, nanoislands, nanovoids, etc.), which depend on the geometrical properties (e.g., the aspect ratio) of the initial pattern. We compare results obtained by using the proposed discrete model with the expectations of the standard continuous theory of surface diffusion, and we show that many qualitative aspects are common to both approaches, in spite of the important differences between both descriptions.

1. Introduction During the last decades, there has been great interest, from both the theoretical and the experimental points of view, in the study of the evolution of micro/nanostructures.1-10 Far from decreasing, this interest has even grown and the changing shapes and interface evolution in the nanoworld have become a very active field of research during the last years.11-16 It is expected that the spread of experimental techniques that allow greater control of materials at the nanoscale and the release of more powerful computers will lead to a deeper comparison between experiments and computer simulations in the near future, further stimulating the interest in surface science at the nanoscale. Within this broad scenario, in this paper, we focused on the evolution of 2-dimensional nanostructures driven by surface diffusion currents. The understanding of the processes involved has a great technological importance, because it is useful not only for the development of nanofabrication mechanisms11 but also for tracking the problem of the stability of nanostructures. In fact, current technologies allow us to achieve a high control on superficial nanostructures, e.g., modern techniques can be used to grow specific nanostructures resolved at the atomic scale.17 Nevertheless, such capabilities in nanofabrication must confront the fact that the obtained structures are often unstable, strongly limiting the field of possible applications. Although the continuous theory of surface evolution driven by surface diffusion, and particularly the decay of periodically corrugated surfaces, is a classical subject dating back to pioneering work by W. W. Mullins18,19 and C. Herring,20,21 up to now, the interest has been mostly focused on the study of the linear theory, which is expected to be relevant only in the so-called limit of small slopes, which is fulfilled when the ratio A/λ f 0, where A and λ are the amplitude and the wavelength * To whom correspondence should be addressed.

of the corrugation, respectively. However, it is worth mentioning that the presence of nonlinearities, the discrete nature of matter, and the crystalline structure of both the substrate and the adsorbate are not considered by the linear continuous theory. Within this context, many studies on the kinetics of profile relaxation and decay have been published in the last two decades. Monte Carlo models, rate equation schemes, statistics of jumps and level set methods,22,23 etc. have been applied in this research field from the theoretical side; moreover, a great deal of experimental work has been carried out too.11,24-28 Nevertheless, in spite of all these efforts, there are many open questions regarding the kinetics of the decay of periodic profiles and, with the exception of the well-understood linear case, a comprehensive theoretical description is still lacking. In this work, we present a kinetic Monte Carlo model for the study of the kinetics of the decay of initially far-from equilibrium periodic patterns in two dimensions. We also compare our findings obtained by studying a stochastic discrete model with expectations from a continuous deterministic vectorial equation for the evolution of surfaces driven by surface diffusion. This paper is organized as follows: In section 2, we discuss the basic physics of the proposed kinetic Monte Carlo model and other issues concerning its implementation. In section 3, we discuss and present results on the continuous, curvaturedriven decay of sinusoidal profiles. Numerical results obtained from the proposed model are discussed in section 4, both for the decay of sinusoidal profiles (4.1) and for the evolution of rectangular patterns (4.2). Finally, section 5 presents a summary and concluding remarks. 2. Model Description The studied system consists of a 2-dimensional geometry of sides Lx × Ly in a triangular lattice on which particles (atoms

10.1021/jp0682311 CCC: $37.00 © 2007 American Chemical Society Published on Web 03/08/2007

Modeling the Decay of Nanopatterns

J. Phys. Chem. C, Vol. 111, No. 12, 2007 4607

or molecules) can diffuse. Each site of the lattice can be in one of two possible states: occupied (by an atom or molecule) or void. A particle in the (i, j) position (1 e i e Lx, 1 e j e Ly) of such lattice can, in a single diffusion event, jump into one of its six neighboring sites, provided that such a site is void. In this way, the dynamics of the system is conservative and the total number of particles NT remains fixed during the whole evolution. We imposed periodic boundary conditions in the x direction and free boundary conditions in the y direction, which implies a rupture of the translational invariance at the lines j ) 1 and j ) Ly. These lines will be called the bottom and the top of the system, respectively. This geometry attempts to describe the case, for example, of the diffusion of particles on a terrace between two rectilinear steps in an fcc solid.29 Also, the proposed geometry is useful for the study of confined particles between two parallel walls. It is worth mentioning that the proposed lattice symmetry and the boundary conditions are fully consistent with those that could be expected in the diffusion of small molecules on the surface of carbon nanotubes functionalized at the ends.30 Nevertheless, modeling such a system would require a careful parametrization and probably the use of some ab initio procedure to get relevant activation energies to achieve reasonable prediction capabilities. Since our aim is to perform a general description of the physical situation, instead of focusing on a particular example, we shall employ a simple harmonic model (described below) to obtain activation barriers. So, we expect that our results may be useful in order to provide a general qualitative (and semiquantitative) knowledge on the subject. Let Eb be the bond energy between two nearest neighbor inf particles and Esup p (Ep ) the bond energy for a particle at the top (bottom) of the system j ) Ly (j ) 1).31 Let nij be the occupation number on the lattice (nij ) 0 if the lattice site is empty and nij ) 1 if it is occupied), and then the Hamiltonian reads

H ) Eb

∑

〈i,j,i′,j′〉

Lx

nijni′j′ +

2Esup p

niLy + ∑ i)1

Lx

2Einf p

ni1 ∑ i)1

(1)

where 〈...〉 denotes a nearest neighbor restricted summation. Although standard Monte Carlo algorithms such as Metropolis, Kawasaki, and Glauber dynamics34-37 are well suited for the study of equilibrium properties of complex systems, it is also well-known that they exhibit a certain degree of arbitrariness for the study of time-dependent physical observables, which is related to the transition probabilities that follow from the master equation describing the Markov process.38 In fact, it has been often observed that physically different results are obtained by switching among these different approaches.39,40 As has been pointed out by Weinberg and co-workers, the dynamics of thermally excited processes can be described by a transition amplitude of the form

ωe(cini, cfin) ∝ exp [-Eact(cini, cfin)/kBT]

(2)

where kB is the Boltzmann constant, T is the absolute temperature, and Eact(cini, cfin) is the energy barrier for a particle to hop, which changes the configuration from cini to cfin. Like Metropolis, Kawasaki, or Glauber samplings, the transition probability given by eq 2 satisfies the detailed balance condition for a Boltzmannian equilibrium state. But, in addition, ωe is a well-defined physical quantity: it is the success probability for thermally activated barrier crossing. The values of the barrier Eact(cini, cfin) used in our simulations are obtained, following to

Figure 1. Activation barriers are obtained from the intersection between the harmonic potentials given by eqs 3 and 4 as a function of the hopping coordinate ζ. More details in the text.

Weinberg et al.,40 by using a simple harmonic model for the energy landscape, as is depicted in Figure 1. In this harmonic model, a coordinate for hopping ζ is introduced, and the following expressions are assumed for the energies of the initial and final states as functions of ζ

k Eini(ζ) ) ζ2 + E0ini 2

(3)

k Efin(ζ) ) (ζ - a)2 + E0fin 2

(4)

where k is the force constant of the harmonic wells, a is the lattice constant, and E0ini and E0fin are the minimum values of the parabolic energy landscapes given by eqs 3 and 4, respectively. These values correspond to the energy that follows from the Hamiltonian given by eq 1, yielding, for a particle at site (i, j) inf E0 ) zEb + 2Esup p δjLy + 2Ep δj1

(5)

where z is the number of occupied nearest-neighbor sites and δ is the Kronecker delta function. Then, the activation energy act Einiffin is given by the difference in energy between the intersection point of Eini(ζ) and Efin(ζ), and E0ini (see Figure 1) act ) Einiffin

(

)

0 0 Efin - Eini 1 + 2 2

2

(6)

where we have called ) ka2. As we have already discussed, we do not attempt to parametrize the model for any particular system, so we take care that not only the selected parameters are within the relevant order of magnitude to account for the involved processes but also that they are suitable to describe physical situations of interest in the field of nanotechnological surface science. Thus, we have taken a ) 0.3 nm, ) 1eV, Eb sup ) -0.1eV, Einf p ) -0.1 eV, and Ep ) 0.1 eV. The negative inf value of Ep assures an attractive interaction at the bottom of the system, whereas the positive value of Esup p makes the top boundary repulsive. In this way, the bottom has the role of a substrate for growth, whereas no significant nucleation of particles is expected to take place at the top of the system. This

4608 J. Phys. Chem. C, Vol. 111, No. 12, 2007

Castez and Albano

situation is similar to the geometry often employed for the study of confined magnetic materials in the presence of competing surface magnetic fields at the confinement walls. The implementation of the algorithmic of the kinetic Monte Carlo model, for a single main loop, can be summarized as follows: • All possible transitions in the system, involving diffusion events from an occupied site to any unoccupied next neighbor site, are identified and their corresponding transition rates are computed. For a transition from state ini to state fin that rate reads act Winiffin ) ν0 exp[-Einiffin /kBT]

(7)

where ν0 is the effective vibration frequency (taken to be 5.21 × 1011 Hz for all the cases treated in this paper). • A diffusion event is determined by randomly choosing among all of the possible jumps weighted by their relative rate of occurrence given by eq 7. This procedure accounts for most of the CPU time and must be optimized as much as possible. In this sense, we implemented a binary tree process that is known to be the most efficient search method developed to date.41 • The transition is performed and the system is updated. When all of the precedent steps have been performed, we say that an iteration of the algorithm has been accomplished. A number of iterations equal to the number of particles in the system (NT) correspond to a Monte Carlo time step (MCS) in the model. Such definition is, to a certain extent, arbitrary, and the important question is: How is the “real time” related to the Monte Carlo time? Although most of the simulational work avoids such question by simply reporting the data by using MCS as unit time, relevant kinetic data require a deeper knowledge about the relationship between both time scales, a subject that has been carefully discussed in refs 39, 40, and 42. In the kinetic Monte Carlo approach, that connection is established in the following way: Let µ be a label for all possible transition events in the system. The transition rate Pµ ) ν0 exp(-(Eµ)/(kBT)) has the dimension of a frequency, and its reciprocal is the residence time for the particle involved in the move in such a specific event. Since those transition probabilities for all possible events are independent, the overall probability per unit time for the system to suffer a transition is obtained by adding all possible transition rates, namely P ) ∑µ Pµ. Then, the quantity 1/P is the residence time for the system in a specific configuration, and consequently, it is the time associated with one iteration (notice that NT of such iterations will correspond to a MCS). Of course, 1/P must be understood as the average time among events. A more realistic description that recovers the stochastic nature of the time interval among events could be obtained by associating a random time variable equal to -1/P ln R with an iteration (assuming exponentially distributed events), where R is a random number uniformly distributed in (0, 1]. Nevertheless, as in this work we do not attempt a statistical characterization of time intervals between events, we speed up the computation simply by taking the average value 1/P as the time per iteration. Figure 2 shows an example of the relationship between real time and MCS as obtained when a cluster evaporates on the 2-dimensional lattice. The “real time” scale is obtained in the computer simulations by adding 1/P, after each iteration, to a variable initialized as equal to zero at the beginning of the simulation, as discussed above. The nonlinear character of this relationship is evident and shows that a precise connection between real time and MCS, as provided by the kinetic Monte Carlo approach, is necessary in kinetic studies.

Figure 2. Plot showing the dependence of the real time on the number of Monte Carlo steps, as obtained during the 2-dimensional evaporation of an initially triangular cluster of Nt ) 15 particles (data correspond to a single run). Typical snapshots that correspond to different stages in the evolution (indicated by the arrows) are also shown.

Even by starting from a single-cluster configuration, one has a rich scenario of possible kinetic evolution. In fact, particles can evaporate from the cluster, voids can appear, and new clusters, unconnected with the original one, can also grow during the dynamic evolution process, etc. In this way, certain programming effort is necessary in order to distinguish between the main cluster (connected to the bottom of the system by a path of nearest neighbors) and other clusters in the system. Let hi be the position of the highest occupied site at the ith column on the main cluster. We shall characterize fluctuations in the interface by a standard estimator such as the interface roughness, (W(Lx, t))

W(Lx, t) )

x∑ 1

Lx

Lx i)1

(hi(t) - hh(t))2

(8)

Lx where hh(t) ) 1/Lx ∑i)1 hi(t). Clearly, the interface roughness is an averaged quantity. In order to characterize the interface, we also used the global width (wg(t)), which is defined as

wg(t) ) Sup{hi(t)} - Inf{hi(t)}

(9)

This nonaveraged quantity is more sensitive to the fluctuations. One has to recognize that by using W(t, Lx) and wg(t) in order to characterize the interface fluctuations, we are losing the information about the formation of voids and overhangs during the evolution. Due to the lack of a proper fluctuating estimator suitable to deal with voids and/or overhangs, in this paper, we present several snapshots of single runs obtained during the simulations, which not only complement the statistical information but also show that overhang formation is the rule rather than the exception upon the decay of patterns in those systems that are out of the small-slope approximation. 3. Brief Overview of the Continuum Description of Surface Diffusion Mediated Relaxation From a mesoscopic point of view, surface evolution mediated by surface diffusion is a well-established topic, dating back to work by Mullins18,19 and Herring.20,21 One of the key contributions from their work is the understanding of the relationship between surface diffusion currents and a geometrical property of the surface: the local curvature. In fact, within this mesoscopic approach, the surface diffusion flux is proportional to



the gradient of the local curvature. Of course, in the derivation of such dependence, there is not only an implicit assumption of smoothness for the surface, as is required by a proper definition of curvature but also an underlying local equilibrium hypothesis to have well-defined thermodynamical coefficients such as surface tension. Already in the pioneering work by Mullins, most attention has been drawn to the case of interfaces in the small-slope approximation (this was mostly due to the resultant simplicity of the equations rather than to physical motivations). Subsequently, also a great deal of attention in the related literature was focused on the following Langevin equation, directly obtained by applying Mullins’s ideas in the small-slope case43-46

∂h(x, t) ∂4h ) -K 4 + η(x, t) ∂t ∂x

Figure 3. Sketch of the evolution of the interface following the normal outward vector according to eq 12, where surface diffusion currents are proportional to the local curvature gradient.

(10)

where h is a single-valued function that describes the interface, K is a constant proportional to the diffusion coefficient, and η is a Gaussian white noise, introduced to capture to some extent the stochastic nature of the diffusion process. By a direct substitution into the linear eq 10 in the absence of noise (i.e., for η ) 0), we can see that at time t a Fourier mode of the initial interface A sin (kx) evolves into

h(x, t) ) A exp(-Kk4t) sin(kx)

(11)

Equation 11 is an interesting result since it shows that every Fourier mode performs an exponential decay whose lifetime depends on temperature through the constant K and on the wavelength λ. In fact, eq 11 shows that the decay constant, in this small-slope case, is proportional to λ-4. The underlying hypothesis of smoothness in thermodynamical quantities such as the surface tension leads to the broadly accepted belief that the applicability of this approach would be restricted to temperatures above the roughening temperatureTR,47 since it is known that the surface tension is discontinuous for temperatures belowTR. The general case in which interfaces are not under the smallslope approximation has been considered with some detail in a recent paper,48 in which a vectorial, stochastic difference equation for the evolution of 2-dimensional interfaces has been studied in a more general context. Let us consider only the deterministic evolution of the model described in ref 48, in the absence of incident flow (i.e., surface diffusion as the only relevant process). Under this restriction, such difference equation reads

δb(s, r t) ) - K

∂2C dtn b ∂s2

(12)

where δr b(s, t) is the local update of the interface when time is advanced in dt, C is the local curvature, s is the arc length parameter, and b n is the local unitary outward normal vector. Evolution proceeds in a Huygens’ construction fashion, familiar from classical optics,49 sketched in Figure 3. Although our present kinetic Monte Carlo model is discrete and stochastic in nature, we shall use, nevertheless, eq 12 for comparison purposes, and, as will be discussed below, although both descriptions are quite different in nature, the results obtained with both models share some common characteristics, at least qualitatively. Figure 4 shows the different relaxation behavior for the decay of initially sinusoidal interfaces of the same wavelength (λ ) 1000 nm): although by starting from small amplitudes, the

Figure 4. log linear plots of the global width of the interface as a function of time. Results obtained by a numerical integration of eq 12 with an initially sinusoidal interface (with λ ) 1000 nm and K ) 1.86 × 106 nm4 s-1) for several initial amplitudes (as listed in the figure in units of nm). The departure from the pure exponential decay mode becomes evident, within the short-time regime, when initial amplitudes are increased, and the initial conditions of the interface lie far from the small-slope case.

global width wg decays exponentially with time, a nonexponential transient takes place when the initial interface is far from the small-slope approximation. As was already discussed in ref 48, this nonexponential decay occurs simultaneously with the spontaneous formation of overhangs. This fact is shown in Figure 5, where an initially sinusoidal profile assumes an overhung shape. Although we were unable to obtain a closed solution to the general nonlinear eq 12, we found that the actual evolution can very closely be approximated by the following vectorial parametric equation

(

b(p, r t) ) p - B(t) sin

(4πpλ ), A(t) sin(2πpλ ))

(13)

where p is the parameter that spans the curve and A(t) and B(t) are coefficients to be fitted from the numerical data. The comparison between the actual evolution obtained by numerically solving eq 12 and that obtained from the Ansatz (13) is given in Figure 5, where the excellent correspondence among curves is evident. 4. Numerical Results Here we present numerical results from the model described in section 2, for two different kinds of initial conditions, namely sine-like patterns (4.1) and rectangular patterns (4.2). It is worth noticing that some results are presented in length units of


Figure 5. Solid lines: Interface profiles obtained by a numerical integration of eq 12 starting from a sinusoidal interface far from the small-slope approximation. The initial sinusoidal profile has λ ) 1000 nm and A ) 1000 nm. Dashed lines: Curves obtained from the application of the Ansatz solution given in eq 13 with fitted values of A and B. Arrows are visual guides to follow the interface evolution.

nanometers (nm), whereas in other cases, the use of lattice units is more convenient. In any case, it is a trivial task to switch among units: for lengths along the x axis one has to multiply the length in lattice units times the lattice constant a ) 0.3 nm, whereas for distances along the y axis an additional factor x3/2 appears due to the geometry of the triangular lattice. Since in this paper we focus on the study of surface diffusion as the main relaxation process, we are restricted to the lowtemperature regime because at higher temperatures evaporationcondensation processes make a non-negligible contribution. Moreover, we have assumed a conserved number of particles in our system, which only is physically realizable at low temperatures, since desorption/evaporation cannot be ignored at high temperatures. It should also be noted that extremely low temperatures are not suitable for our study, since in that case the system spends most of the time cycling between generating/ annihilating very few excitations, and, in this way, in order to obtain relevant information about global changes occurring in the system, one needs to spend prohibitively long CPU times. So, most of the reported results are within the temperature range 100-400 K. 4.1. Decay of Sinusoidal Patterns. In this section, we present numerical results from our kinetic Monte Carlo model for the decay of initially sinusoidal profiles. For profiles within the small-slope approximation (for instance, sinusoidal profiles with an amplitude/wavelength ratio smaller than 0.1), we found an exponential decay of the pattern amplitudes. Moreover, the decaying constant depends on the wavelength as κ ∝ λ-4, i.e., in complete agreement with the continuous linear theory expectations (see eq 11). All of these results are summarized in Figure 6, for two different temperatures. These results are essentially independent of the temperature within the considered range 100 e T < 400, although for the higher temperatures (i.e., near 400 K) the exponent in the dependence between κ and λ is somewhat lower than 4 (in absolute value). This latter result can easily be interpreted because, at higher temperatures, the evaporation-condensation process starts to become relevant and, as it is well-known, when the decay of the profile is dominated by evaporation-condensation processes, the exponent in the power-law relation between κ and λ is no longer 4 but 2.4,50 For initial sinusoidal profiles far from the small-slope approximation, we obtain a spontaneous overhang formation

Castez and Albano

Figure 6. Exponential decay of W2 as a function of time, obtained for two different temperatures (200 K above, 300 K below) and for different wavelengths as listed in the figures. Axes were properly rescaled to obtain data collapse. Data corresponding to the small-slope case (the aspect ratio is A/λ ) 0.1). The insets show log-log plots of the decaying constant as a function of the wavelength λ and the full lines have slope -4.

Figure 7. Snapshots of the interface at successive times (from 0 to 500 MCS) showing the decay of a sinusoidal profile far from the smallslope approximation. Notice the qualitative agreement with the decay obtained in the continuous approach, shown in Figure 5. Results obtained at T ) 300 K by using geometrical parameters in the initial wave given by (in lattice units) A ) 20 and λ ) 20.

in the intermediate state, before the decay of the interface into the equilibrium nearly flat shape. This can be observed in the typical snapshots (taken every 100 MCS and obtained at T ) 300 K) shown in Figure 7, which correspond to the decay of a sinusoidal interface with A ) 20 and λ ) 20 (in lattice units) so that A/λ ) 1 . 0.1. The shape of the interface at different stages is very similar to that predicted by the continuous eq 12 (see Figure 5), although of course, some differences become also evident due to the stochastic nature of the Monte Carlo model and some faceting that can be observed in the snapshots of Figure 7. This faceting is more pronounced at lower temperatures, as can be observed in the snapshots shown in Figure 8, where the decay of the profile is followed at 100 K. Although the interface is rather faceted in this case, as expected for the very-low-temperature regime, the qualitative behavior is the same as that observed at higher temperatures (e.g., T ) 300 K in Figure 7), in the sense that overhung profiles also appear in the intermediate stages of the temporal evolution. In Figure 9, we compare the evolution of the global width wg as a function of time as obtained by using the kinetic (discrete) Monte Carlo model with the predictions due to the


Figure 8. Snapshots of the interface at successive times (from 0 to 3000 MCS) showing the decay of a sinusoidal profile far from the small-slope approximation and for very low temperature (100 K). In spite of the fact that the interface is rather faceted, one can observe the spontaneous formation of overhangs in the transient regime (at the top right), in qualitative agreement with the expectations of the continuous approach. The geometrical parameters in the initial wave are (in lattice units) A ) 10 and λ ) 10.

Figure 9. Interface global width as a function of time for initially sinusoidal profiles (λ ) 50 lattice units) with different amplitudes (all of them indicated in lattice units at the legends). On the left: Evolution obtained from the discrete kinetic Monte Carlo model. On the right: Evolution obtained by a numerical integration of the continuous eq 12.

continuous description (eq 12) for several initial amplitudes. Although there are some differences (for example, the curves corresponding to the Monte Carlo approach are noisy since the model is stochastic, and the asymptotic value of wg is not zero due to thermal fluctuations even in the equilibrium state), the curves are qualitatively similar and give evidence of the departure from the exponential decay that takes place when the initial amplitudes are not consistent with the small-slope approximation. 4.2. Decay of Rectangular Patterns. In this section, we present some typical results on the time evolution of initially rectangular patterns. Let Bs and Bi be the upper and lower basis of the pattern and H the height of the rectangles. The spatial period of the pattern then is Bs + Bi. More specifically, we are interested in the relaxation behavior of patterns in the largeslope case (i.e., for H g Bs, Bi), and varying the ratio between the bases Bs and Bi. In Figure 10, we show typical snapshots of the decay of the profiles of an initially rectangular pattern with Bs ) Bi ) 10, and H ) 30. It is observed that after a short transient, the profile evolves with a shape that resembles the evolution obtained when


Figure 10. Decay of a rectangular pattern in the case Bs ) Bi (from 0 to 500 MCS). Transient states are completely similar to that obtained starting from sinusoidal patterns. Temperature in the simulation is 300 K, and geometrical parameters of the initial pattern are (in lattice units) H ) 30 and Bs ) Bi ) 10.

Figure 11. Interface global width as a function of time for a rectangular pattern as initial condition (Bs ) Bi ) 10, in lattice units) with different amplitudes (all of them indicated in lattice units at the legends). Temperature in the simulation is 300 K and data are averaged over 10 independent runs. In the inset: Comparison between the decay of the global width starting from a rectangular pattern and from a sinusoidal one with the same amplitude (data in the sinusoidal case are shifted along the x axis to stress the correspondence with the rectangular case after a transient).

the initial profile was a (high slope) sinusoidal profile. The onset of coarsening between adjacent rectangles, which becomes evident in the latest state of the relaxation (see Figure 10) has also been observed in the case of sinusoidal patterns (Figure 7). Therefore, in this case again one expects that the interface may qualitatively be well described by the parametrization given by eq 13 and plotted in Figure 5. Thus, the overall behavior of the evolution is well described by applying the superposition principle, with fast filtering of the high-frequency modes present in the interface. This is a very interesting result, because the pattern of Figure 10 corresponds to a strongly nonlinear case, and then superposition ideas are not expected to hold, in principle. Furthermore, as is shown in the Figure 11, the global width increases at short times. Such initial increase in the interface width for an initially rectangular interface has been found already both in theoretical48 and in experimental24 work. Figure 11 also shows that, for small amplitudes (small H values when compared to Bs + Bi) wg increases during a short transient that is followed by an exponential decay, as is expected for a small-


Castez and Albano

Figure 12. Decay of a rectangular pattern in the case Bs , Bi (from 0 to 5000 MCS). A transient state consisting of a matrix-like array of nanoislands is clearly observed. Temperature in the simulation is 300 K, and geometrical parameters of the initial pattern are (in lattice units) H ) 70, Bs ) 4 and Bi ) 16.

Figure 13. Decay of a rectangular pattern in the case Bi , Bs (from 0 to 500 MCS). The formation of a transient state consisting of a matrixlike array of nanovoids is evident. Temperature in the simulation is 300 K, and geometrical parameters of the initial pattern are (in lattice units) H ) 70, Bs ) 16 and Bi ) 4.

amplitude sinusoidal profile. When the amplitudes get larger, there is an intermediate stage where the decay is no longer exponential; in fact for the largest H values, Figure 11 shows a rather slow decay similar to that already found for largeamplitude sine profiles (see Figure 9). The inset of Figure 11 shows a comparison between the decay of the rectangular pattern and a sinusoidal one, both with the same geometrical parameters. Notice that data from the sinusoidal case are shifted along the x axis to stress the fact that both systems follow nearly the same evolution after a transient time (of course, the transient time is different in each case). On the other hand, when Bs , Bi, the main cluster is unstable and after a short transient every rectangle in the pattern breaks down into several islands, the number of such islands being a stochastic quantity, with an average proportional to the height of the rectangles H. The average size of the islands is also proportional to the width of the rectangles Bs. In this way, the transient state is conformed by a number of small clusters arranged in a matrix-like structure, as can be seen in Figure 12. In the opposite case, when Bi , Bs, fluctuations on the sides of adjacent rectangles cause them to merge, forming voids. The situation is closely related to the case of emergence of nanoislands, and it can be thought as some kind of island-void symmetry, since both are complementary to a certain degree: The average size of the nanovoids is proportional to Bi, whereas the average number of voids formed between two rectangles is proportional to H. The whole picture is that of a matrix-like arrangement of nanovoids appearing after a short transient, as is shown in Figure 13. Of course, the described states are metastable because, waiting for a long enough time, both nanoislands and nanovoids would diffuse and coalesce forming a main cluster, with a rather compact structure, in order to decrease the energy of the system. Nevertheless, it is important to stress that such metastable structures can survive for long periods, and even can become virtually frozen by lowering the temperature once the metastable state is established.

simulations by taking the reciprocal of the sum over all possible transition rates and then computing the residence time for the system in a specific configuration. We found a nonlinear relationship between time scales, namely the real time and the Monte Carlo steps. This fact shows that a proper choice of time units in Monte Carlo models is essential for the study of kinetic properties. Results from our discrete model were compared to predictions from the continuous theory of surface evolution by surface diffusion. Particular attention was payed to both the transition between the exponential and the nonexponential decay modes when initial sinusoidal profiles depart from the small-slope approximation, and to the morphological transition in which the spontaneous emergence of overhangs occurs as a transient state. All of these phenomena, previously reported for a continuous model,48 were also observed in the present discrete model. Moreover, starting from sine-like interfaces satisfying the small-slope approximation, we recover the well-known result from the linear theory: an exponential decay of the amplitudes with a decay constant that depends on the wavelength as κ ∝ λ-4. This agreement between results obtained for two very different modeling strategies not only proves the robustness of the assumptions involved but also has its own theoretical relevance because that the relationship between continuous descriptions and discrete models is still under debate in several physical situations. We have also studied the decay of initially rectangular nanostructures, consisting of “columns” of width Bs separated by a distance Bi. Our simulations show that for Bs , Bi (Bs . Bi) a matrix-like arrangement of nanoislands (nanovoids) emerges. The average number and average size of these islands and voids depend on the geometry of the initial pattern. The same applies to the duration of the transient state, although, of course, such nanostructures can be frozen by lowering the temperature. For Bi ∼ Bs, the pattern evolves, for intermediate times, in a fashion similar to the case of an initially sinusoidal pattern with a wavelength equal to Bi + Bs. It is worth mentioning that this behavior provides evidence that the highfrequency modes of the interface are filtered fast, as could be expected from the linear theory, although we are starting from a situation very far from the small-slope approximation and, consequently, in such cases the resulting evolution is strongly nonlinear.

5. Concluding Remarks We have studied a kinetic Monte Carlo model for the relaxation of 2-dimensional nanostructures by surface diffusion. Activation barriers for atom hopping were computed by using a simple harmonic model. Real time is introduced in the

Modeling the Decay of Nanopatterns Our results show that, for both the discrete Monte Carlo approach and the continuum description, overhang formation during the intermediate stages of the equilibration process is quite frequent, so that it may be considered the rule rather than an exception. In this sense, we stress the importance of modeling this kind of system in a way that accounts for the possible formation of overhangs during the dynamic evolution. In contrast, highly restrictive approaches like those involved in the formulation of standard solid-on-solid models,43 that by definition do not allow for overhangs, cannot describe either multivalued interfaces and void formation or other interesting physical situations. To the best of our knowledge, the emergence of overhangs and, more generally, the decay of structures by surface diffusion far from the small-slope approximation have not systematically been studied yet from the experimental side. Within this context, we expect that experiments performed by working on truly 2-dimensional nanostructures or cross-sectional studies on real 3-dimensional systems, would be helpful to increase our knowledge of several technologically relevant processes at the nanoscale. Acknowledgment. This work is supported by CONICET, ANPCyT, and UNLP. We acknowledge fruitful discussions with Dr. Ezequiel P. M. Leiva and Dr. Roberto C. Salvarezza. Authors are members of the CONICET (Argentina National Research Council). References and Notes (1) Bales, G. S.; Redfield, A. C.; Zangwill, A. Phys. ReV. Lett. 1989, 62, 776. (2) Dubson, M. A.; Jeffers, A. Phys. ReV. B 1994, 49, 8347. (3) Erlebacher, J.; Aziz, M. J.; Chason, E.; Sinclair, M. B.; Floro, J. A. Phys. ReV. Lett. 2000, 84, 5800. (4) Family, F. J. Phys. A 1986, 19, L441. (5) Fu, E. S.; Johnson, M. D.; Liu, D.-J.; Weeks, J. D.; Williams, E. D. Phys. ReV. Lett. 1996, 77, 1091. (6) Israeli, N.; Kandel, D. Phys. ReV. Lett. 1998, 80, 3300. (7) Kallabis, H.; Wolf, D. E. Phys. ReV. Lett. 1997, 79, 4854. (8) Karunasiri, R. P.; Bruinsma, R.; Rudnick, J. Phys. ReV. Lett. 1989, 62, 788. (9) Maritan, A.; Toigo, F.; Koplic, J.; Banavar, J. R. Phys. ReV. Lett. 1992, 69, 3193. (10) Son, C.-S.; Kim, T. G.; Wang, X.-L.; Ogura, M. J. Cryst. Growth 2000, 221, 201. (11) Andreasen, G.; Schilardi, P. L.; Azzaroni, O.; Salvarezza, R. C. Langmuir 2002, 18, 10430. (12) Castez, M. F.; Fonticelli, M. H.; Azzaroni, O.; Salvarezza, R. C.; Solari, H. G. Appl. Phys. Lett. 2005, 87, 123104. (13) Israeli, N.; Kandel, D. Phys. ReV. Lett. 2002, 88, 116103. (14) Liu, D.; Evans, J. W. Phys. ReV. B 2002, 66, 165407.

J. Phys. Chem. C, Vol. 111, No. 12, 2007 4613 (15) Murty, M. V. R. Phys. ReV. B 2000, 62, 17004. (16) Rieth, M. Nano-Engineering in Science and Technology; World Scientific: Singapore, 2003. (17) Kolb, M.; Ullmann, R.; Will, T. Science 1997, 275, 1097. (18) Mullins, W. W. J. Appl. Phys. 1957, 28, 333. (19) Mullins, W. W. J. Appl. Phys. 1959, 30, 77. (20) Herring, C. In Physics of Powder Metallurgy; Kingston, W. E., Ed.; McGraw-Hill Book Company Inc.: New York, 1951. (21) Herring, C. In Structure and Properties of Solid Surfaces; Gomer, R., Smith, C. S., Eds.; The University of Chicago Press: Chicago, 1952. (22) Baumann, F. H.; Chopp, D. L.; Dı´az de la Rubia, T.; Gilmer, G. H.; Greene, J. E.; Huang, H.; Kodambaka, S.; O’Sullivan, P.; Petrov, I. MRS Bull. 2001, 26 (3), 182. (23) Osher, S.; Sethian, J. A. J. Comp. Phys. 1988, 79, 12. (24) Kan, H.-C.; Shah, S.; Tadyyon-Eslami, T.; Phaneuf, R. J. Phys. ReV. Lett. 2004, 92, 146101. (25) Bonzel, H. P.; Latta, E. E. Surf. Sci. 1978, 76, 275. (26) Bonzel, H. P.; Preuss, E. Surf. Sci. 1995, 336, 209. (27) Giesen-Seibert, M.; Jentjens, R.; Poensgen, M.; Ibach, H. Phys. ReV. Lett. 1993, 71, 3521. (28) Giesen-Seibert, M.; Ibach, H. Surf. Sci. 1994, 316, 205. (29) Although a careful modelization of such system would take into account the fact that the six neighboring sites of a given site on an fcc lattice are not equivalent, contrarily to our assumptions. (30) Kim, C.; Seo, K.; Kim, B.; Park, N.; Choi, Y. S.; Park, K. A.; Lee, Y. H. Phys. ReV. B 2003, 68, 115403. (31) Negative values of these energies correspond to attractive interacinf tions. Positive (repulsive case) values of Esup p or Ep could be expected in the case, for example, of descendant steps, with associated Schwoebel barriers.32,33 (32) Schwoebel, R. L.; Shipsey, E. J. J. Appl. Phys. 1967, 37, 3682. (33) Schwoebel, R. L. J. Appl. Phys. 1969, 40, 614. (34) Glauber, R. J. J. Math. Phys. 1963, 4, 294. (35) Heermann, D. W. Computer Simulation Methods; Springer-Verlag: Heidelberg, Germany, 1989. (36) Kawasaki, K. In Phase Transitions and Critical Phenomena; Domb, C., Green, M. S., Eds.; Academic: New York, 1972; Vol. 2. (37) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087. (38) Van Kampen, N. G. Stochastic Processes in Physics and Chemistry; Elsevier Science B.V.: Amsterdam, The Netherlands, 1992. (39) Kang, H. C.; Weinberg, W. H. Phys. ReV. B 1988, 38, 11543. (40) Kang, H. C.; Weinberg, W. H. J. Chem. Phys. 1989, 90, 2824. (41) Yang, Y. Ph.D. Thesis, University of Virginia, 2000; advisor: H. N. G. Wadley. (42) Fichthorn, K. A.; Weinberg, W. H. J. Chem. Phys. 1991, 95, 1090. (43) Barabasi, A. L.; Stanley, H. E. Fractal concepts in surface growth; Cambridge University Press: Cambridge, 1995. (44) Das Sarma, S.; Tamborenea, P. Phys. ReV. Lett. 1991, 66, 325. (45) Lai, Z.-W.; Das Sarma, S. Phys. ReV. Lett. 1991, 66, 2348. (46) Wolf, D. E.; Villain, J. Europhys. Lett. 1990, 13, 389. (47) Lapujoulade, J. Surf. Sci. Rep. 1994, 20, 191. (48) Castez, M. F.; Salvarezza, R. C.; Solari, H. G. Phys. ReV. E 2006, 73, 011607. (49) Hecht, E.; Zajac, A. Optics; Addison-Wesley Publishing Company: Massachusetts, 1974. (50) Edwards, S. F.; Wilkinson, D. R. Proc. R. Soc. London A 1982, 381, 17.

Modeling the Decay of Nanopatterns: A ... - American Chemical Society

Recommend Documents