Numerical study of discrete models in the class of the nonlinear molecular beam epitaxy equation

Scaling of surface roughness in film deposition with height-dependent step edge barriers

We perform kinetic Monte Carlo simulations of film growth in simple cubic lattices with solid-on-solid conditions, Ehrlich-Schwoebel (ES) barriers at step edges, and a kinetic barrier related to the hidden off-plane diffusion at multilayer steps. Broad ranges of the diffusion-to-deposition ratio R, detachment probability per lateral neighbor, ε, and monolayer step crossing probability P=exp[-E_{ES}/(k_{B}T)] are studied. Without the ES barrier, four possible scaling regimes are shown as the coverage θ increases: nearly layer-by-layer growth with damped roughness oscillations; kinetic roughening in the Villain-Lai-Das Sarma (VLDS) universality class when the roughness is W∼1 (in lattice units); unstable roughening with mound nucleation and growth, where slopes of logW×logθ plots reach values larger than 0.5; and asymptotic statistical growth with W=θ^{1/2} solely due to the kinetic barrier at multilayer steps. If the ES barrier is present, the layer-by-layer growth crosses over directly to the unstable regime, with no transient VLDS scaling. However, in simulations up to θ=10^{4} (typical of films with a few micrometers), low temperatures (small R, ε, or P) may suppress the two or three initial regimes, while high temperatures and P∼1 produce smooth surfaces at all thicknesses. These crossovers help to explain proposals of nonuniversal exponents in previous works. We define a smooth film thickness θ_{c} where W=1 and show that VLDS scaling at that point indicates negligible ES barriers, while rapidly increasing roughness indicates a small ES barrier (E_{ES}∼k_{B}T). θ_{c} scales as ∼exp(const×P^{2/3}) if the other parameters are kept fixed, which represents a high sensitivity on the ES barrier. The analysis of recent experimental data in the light of our results distinguishes cases where E_{ES}/(k_{B}T) is negligible, ∼1, or ≪1.

Height distributions in interface growth: The role of the averaging process

Height distributions (HDs) are key quantities to uncover universality and geometry-dependence in evolving interfaces. To quantitatively characterize HDs, one uses adimensional ratios of their first central moments (m_{n}) or cumulants (κ_{n}), especially the skewness S and kurtosis K, whose accurate estimate demands an averaging over all L^{d} points of the height profile at a given time, in translation-invariant interfaces, and over N independent samples. One way of doing this is by calculating m_{n}(t) [or κ_{n}(t)] for each sample and then carrying out an average of them for the N interfaces, with S and K being calculated only at the end. Another approach consists in directly calculating the ratios for each interface and, then, averaging the N values. It turns out, however, that S and K for the growth regime HDs display strong finite-size and -time effects when estimated from these "interface statistics," as already observed in some previous works and clearly shown here, through extensive simulations of several discrete growth models belonging to the EW and KPZ classes on one- and two-dimensional substrates of sizes L=const. and L∼t. Importantly, I demonstrate that with "1-point statistics," i.e., by calculating m_{n}(t) [or κ_{n}(t)] once for all NL^{d} heights together, these corrections become very weak, so that S and K attain values very close to the asymptotic ones already at short times and for small L's. However, I find that this "1-point" (1-pt) approach fails in uncovering the universality of the HDs in the steady-state regime (SSR) of systems whose average height, h[over ¯], is a fluctuating variable. In fact, as demonstrated here, in this regime the 1-pt height evolves as h(t)=h[over ¯](t)+s_{λ}A^{1/2}L^{α}ζ+⋯-where P(ζ) is the underlying SSR HD-and the fluctuations in h[over ¯] yield S_{1-pt}∼t^{-1/2} and K_{1-pt}∼t^{-1}. Nonetheless, by analyzing P(h-h[over ¯]), the cumulants of P(ζ) can be accurately determined. I also show that different, but universal, asymptotic values for S and K (related, so, to different HDs) can be found from the "interface statistics" in the SSR. This reveals the importance of employing the various complementary approaches to reliably determine the universality class of a given system through its different HDs.

Modeling of nonequilibrium surface growth by a limited-mobility model with distributed diffusion length

Kinetic Monte Carlo (KMC) simulations are a well-established numerical tool to investigate the time-dependent surface morphology in molecular beam epitaxy experiments. In parallel, simplified approaches such as limited mobility (LM) models characterized by a fixed diffusion length have been studied. Here we investigate an extended LM model to gain deeper insight into the role of diffusional processes concerning the growth morphology. Our model is based on the stochastic transition rules of the Das Sarma-Tamborena model but differs from the latter via a variable diffusion length. A first guess for this length can be extracted from the saturation value of the mean-squared displacement calculated from short KMC simulations. Comparing the resulting surface morphologies in the sub- and multilayer growth regime to those obtained from KMC simulations, we find deviations which can be cured by adding fluctuations to the diffusion length. This mimics the stochastic nature of particle diffusion on a substrate, an aspect which is usually neglected in LM models. We propose to add fluctuations to the diffusion length by choosing this quantity for each adsorbed particle from a Gaussian distribution, where the variance of the distribution serves as a fitting parameter. We show that the diffusional fluctuations have a huge impact on cluster properties during submonolayer growth as well as on the surface profile in the high coverage regime. The analysis of the surface morphologies on one- and two-dimensional substrates during sub- and multilayer growth shows that the LM model can produce structures that are indistinguishable to the ones from KMC simulations at arbitrary growth conditions.

Effects of a kinetic barrier on limited-mobility interface growth models

The role played by a kinetic barrier originated by out-of-plane step edge diffusion, introduced by Leal et al. [J. Phys.: Condens. Matter 23, 292201 (2011)JCOMEL0953-898410.1088/0953-8984/23/29/292201], is investigated in the Wolf-Villain and Das Sarma-Tamborenea models with short-range diffusion. Using large-scale simulations, we observe that this barrier is sufficient to produce growth instability, forming quasiregular mounds in one and two dimensions. The characteristic surface length saturates quickly indicating a uncorrelated growth of the three-dimensional structures, which is also confirmed by a growth exponent β=1/2. The out-of-plane particle current shows a large reduction of the downward flux in the presence of the kinetic barrier enhancing, consequently, the net upward diffusion and the formation of three-dimensional self-assembled structures.

Local roughness exponent in the nonlinear molecular-beam-epitaxy universality class in one dimension

We report local roughness exponents, α_{loc}, for three interface growth models in one dimension which are believed to belong to the nonlinear molecular-beam-epitaxy (nMBE) universality class represented by the Villain-Lais-Das Sarma (VLDS) stochastic equation. We applied an optimum detrended fluctuation analysis (ODFA) [Luis et al., Phys. Rev. E 95, 042801 (2017)2470-004510.1103/PhysRevE.95.042801] and compared the outcomes with standard detrending methods. We observe in all investigated models that ODFA outperforms the standard methods providing exponents in the narrow interval α_{loc}^{}∈[0.96,0.98] quantitatively consistent with two-loop renormalization group predictions for the VLDS equation. In particular, these exponent values are calculated for the Clarke-Vvdensky and Das Sarma-Tamborenea models characterized by very strong corrections to the scaling, for which large deviations of these values had been reported. Our results strongly support the absence of anomalous scaling in the nMBE universality class and the existence of corrections in the form α_{loc}^{}=1-ε of the one-loop renormalization group analysis of the VLDS equation.

Crossover from compact to branched films in electrodeposition with surface diffusion

We study a model for thin film electrodeposition in which instability development by preferential adsorption and reduction of cations at surface peaks competes with surface relaxation by diffusion of the adsorbates. The model considers cations moving in a supported electrolyte, adsorption and reduction when they reach the film surface, and consequent production of mobile particles that execute activated surface diffusion, which is represented by a sequence of random hops to neighboring lattice sites with a maximum of G hop attempts (G≫1), a detachment probability ε<1 per neighboring particle, and a no-desorption condition. Computer simulations show the formation of a compact wetting layer followed by the growth of branched deposits. The maximal thickness z_{c} of that layer increases with G but is weakly affected by ε. A scaling approach describes the crossover from smooth film growth to unstable growth and predicts z_{c}∼G^{γ}, with γ=1/[2(1-ν)]≈0.43, where ν≈0.30 is the inverse of the dynamical exponent of the Villain-Lai-Das Sarma equation that describes the initial roughening. Using previous results for related deposition models, the thickness z_{c} can be predicted as a function of an activation energy for terrace surface diffusion and the temperature, and the small effects of the parameter ε are justified. These predictions are confirmed by the numerical results with good accuracy. We discuss possible applications, with a particular focus on the growth of multifuncional structures with stacking layers of different porosity.

Optimal detrended fluctuation analysis as a tool for the determination of the roughness exponent of the mounded surfaces

We present an optimal detrended fluctuation analysis (DFA) and apply it to evaluate the local roughness exponent in nonequilibrium surface growth models with mounded morphology. Our method consists in analyzing the height fluctuations computing the shortest distance of each point of the profile to a detrending curve that fits the surface within the investigated interval. We compare the optimal DFA (ODFA) with both the standard DFA and nondetrended analysis. We validate the ODFA method considering a one-dimensional model in the Kardar-Parisi-Zhang universality class starting from a mounded initial condition. We applied the methods to the Clarke-Vvedensky (CV) model in 2+1 dimensions with thermally activated surface diffusion and absence of step barriers. It is expected that this model belongs to the nonlinear molecular beam epitaxy (nMBE) universality class. However, an explicit observation of the roughness exponent in agreement with the nMBE class was still missing. The effective roughness exponent obtained with ODFA agrees with the value expected for the nMBE class, whereas using the other methods it does not agree. We also characterize the transient anomalous scaling of the CV model and obtained that the corresponding exponent is in agreement with the value reported for other nMBE models with weaker corrections to the scaling.

Effects of film growth kinetics on grain coarsening and grain shape

We study models of grain nucleation and coarsening during the deposition of a thin film using numerical simulations and scaling approaches. The incorporation of new particles in the film is determined by lattice growth models in three different universality classes, with no effect of the grain structure. The first model of grain coarsening is similar to that proposed by Saito and Omura [Phys. Rev. E 84, 021601 (2011)PLEEE81539-375510.1103/PhysRevE.84.021601], in which nucleation occurs only at the substrate, and the grain boundary evolution at the film surface is determined by a probabilistic competition of neighboring grains. The surface grain density has a power-law decay, with an exponent related to the dynamical exponent of the underlying growth kinetics, and the average radius of gyration scales with the film thickness with the same exponent. This model is extended by allowing nucleation of new grains during the deposition, with constant but small rates. The surface grain density crosses over from the initial power law decay to a saturation; at the crossover, the time, grain mass, and surface grain density are estimated as a function of the nucleation rate. The distributions of grain mass, height, and radius of gyration show remarkable power law decays, similar to other systems with coarsening and particle injection, with exponents also related to the dynamical exponent. The scaling of the radius of gyration with the height h relative to the base of the grain show clearly different exponents in growth dominated by surface tension and growth dominated by surface diffusion; thus it may be interesting for investigating the effects of kinetic roughening on grain morphology. In growth dominated by surface diffusion, the increase of grain size with temperature is observed.

Universality and dependence on initial conditions in the class of the nonlinear molecular beam epitaxy equation

We report extensive numerical simulations of growth models belonging to the nonlinear molecular beam epitaxy (nMBE) class, on flat (fixed-size) and expanding substrates (ES). In both d=1+1 and 2+1, we find that growth regime height distributions (HDs), and spatial and temporal covariances are universal, but are dependent on the initial conditions, while the critical exponents are the same for flat and ES systems. Thus, the nMBE class does split into subclasses, as does the Kardar-Parisi-Zhang (KPZ) class. Applying the "KPZ ansatz" to nMBE models, we estimate the cumulants of the 1+1 HDs. Spatial covariance for the flat subclass is hallmarked by a minimum, which is not present in the ES one. Temporal correlations are shown to decay following well-known conjectures.

Smoothening in thin-film deposition on rough substrates

The evolution of the surface roughness W of a thin film deposited on a rough substrate is studied with a model of temperature-activated adatom diffusion, irreversible lateral aggregation, and no step energy barrier, in which the main parameter is the ratio R of diffusion and deposition rates. At sufficiently low temperatures (R≲10), the average number of adatom steps after adsorption is very small, thus W monotonically increases with time t due to an approximately uncorrelated deposition at short times. If the temperature is not very low (R∼10(3) or larger), smoothening occurs at short times and the Villain-Lai-Das Sarma (VLDS) growth equation governs the long time roughening, which is attained after a crossover time t(c) that increases with the correlation length ξ(i) of the substrate. Scaling arguments predict the dependence of t(c) on temperature and on the substrate production time and the scaling relation for the difference between the roughness of films deposited on rough and flat substrates, in good agreement with numerical results. The effect of temperature is not a direct extension of previous results on flat substrates because the short wavelength fluctuations delay the formation of terraces. For this reason, the effective energy obtained from the dependence of t(c) on R is 40% of the energy of activated adatom diffusion. A scaling law for the initial smoothening is proposed as W/W(i)=Ψ(t/t(c1)), with a crossover time t(c1)≡R(-θ)ξ(i)(z), where W(i) is the substrate roughness, θ≈0.4, and z is the VLDS dynamical exponent. It provides good data collapse if W is not very small and is suggested to be tested experimentally.

Transition from compact to porous films in deposition with temperature-activated diffusion

We study a thin-film growth model with temperature activated diffusion of adsorbed particles, allowing for the formation of overhangs and pores, but without detachment of adatoms or clusters from the deposit. Simulations in one-dimensional substrates are performed for several values of the diffusion-to-deposition ratio R of adatoms with a single bond and of the detachment probability ε per additional nearest neighbor, respectively, with activation energies are E(s) and E(b). If R and ε independently vary, regimes of low and high porosity are separated at 0.075≤ε(c)≤0.09, with vanishingly small porosity below that point and finite porosity for larger ε. Alternatively, for fixed values of E(s) and E(b) and varying temperature, the porosity has a minimum at T(c), and a nontrivial regime in which it increases with temperature is observed above that point. This is related to the large mobility of adatoms, resembling features of equilibrium surface roughening. In this high-temperature region, the deposit has the structure of a critical percolation cluster due to the nondesorption. The pores are regions enclosed by blobs of the corresponding percolating backbone, thus the distribution of pore size s is expected to scale as s(-τ̃) with τ̃≈1.45, in reasonable agreement with numerical estimates. Roughening of the outer interface of the deposits suggests Villain-Lai-Das Sarma scaling below the transition. Above the transition, the roughness exponent α≈0.35 is consistent with the percolation backbone structure via the relation α=2-d(B), where d(B) is the backbone fractal dimension.

Skewness in (1+1)-dimensional Kardar-Parisi-Zhang-type growth

We use the (1+1)-dimensional Kardar-Parisi-Zhang equation driven by a Gaussian white noise and employ the dynamic renormalization-group of Yakhot and Orszag without rescaling [J. Sci. Comput. 1, 3 (1986)]. Hence we calculate the second- and third-order moments of height distribution using the diagrammatic method in the large-scale and long-time limits. The moments so calculated lead to the value S=0.3237 for the skewness. This value is comparable with numerical and experimental estimates.

Nonuniversal effects in mixing correlated-growth processes with randomness: interplay between bulk morphology and surface roughening

To construct continuum stochastic growth equations for competitive nonequilibrium surface-growth processes of the type RD+X that mixes random deposition (RD) with a correlated-growth process X, we use a simplex decomposition of the height field. A distinction between growth processes X that do and do not create voids in the bulk leads to the definition of the effective probability p(eff) of the process X that is a measurable property of the bulk morphology and depends on the activation probability p of X in the competitive process RD+X. The bulk morphology is reflected in the surface roughening via nonuniversal prefactors in the universal scaling of the surface width that scales in p(eff). The equation and the resulting scaling are derived for X in either a Kardar-Parisi-Zhang or Edwards-Wilkinson universality class in (1+1) dimensions and are illustrated by an example of X being a ballistic deposition. We obtain full data collapse on its corresponding universal scaling function for all p∈(0;1]. We outline the generalizations to (1+n) dimensions and to many-component competitive growth processes.

Normal dynamic scaling in the class of the nonlinear molecular-beam-epitaxy equation

The scaling of local height fluctuations is studied numerically in lattice growth models of the class of the nonlinear stochastic equation of Villain-Lai-Das Sarma (VLDS) in substrate dimensions d=1 and 2. In d=1, the average local slopes of the conserved restricted solid-on-solid (CRSOS) models converge to a finite value in the long-time limit, with power-law corrections in time whose exponents are close to 0.1. Other VLDS models in d=1, such as that of Das Sarma and Tamborenea, show a divergence of local slopes up to 10(6) monolayers, typical of anomalous roughening, but a comparison of roughness distributions shows that they scale as the linear fourth-order growth equation in those time scales. Normal scaling is also obtained in a modified VLDS equation with instability suppression, in contrast to recent numerical works. In d=2, a CRSOS model and a model with lateral aggregation of diffusing particles show normal scaling of the local slopes, also with small correction exponents. These results consistently show that the VLDS class has normal dynamic scaling in d=1 and 2, in agreement with the theoretical predictions of Phys. Rev. Lett. 94, 166103 (2005), and they show that the apparently anomalous features observed in previous works are effects of large scaling correction terms or crossover effects.

Roughness exponents and grain shapes

In surfaces with grainy features, the local roughness w shows a crossover at a characteristic length r(c), with roughness exponent changing from α(1)≈1 to a smaller α(2). The grain shape, the choice of w or height-height correlation function (HHCF) C, and the procedure to calculate root-mean-square averages are shown to have remarkable effects on α(1). With grains of pyramidal shape, α(1) can be as low as 0.71, which is much lower than the previous prediction 0.85 for rounded grains. The same crossover is observed in the HHCF, but with initial exponent χ(1)≈0.5 for flat grains, while for some conical grains it may increase to χ(1)≈0.7. The universality class of the growth process determines the exponents α(2)=χ(2) after the crossover, but has no effect on the initial exponents α(1) and χ(1), supporting the geometric interpretation of their values. For all grain shapes and different definitions of surface roughness or HHCF, we still observe that the crossover length r(c) is an accurate estimate of the grain size. The exponents obtained in several recent experimental works on different materials are explained by those models, with some surface images qualitatively similar to our model films.

Dynamic scaling in thin-film growth with irreversible step-edge attachment

We study dynamic scaling in a model with collective diffusion (CD) of isolated atoms in terraces and irreversible aggregation at step edges. Simulations are performed in two-dimensional substrates with several diffusion to deposition ratios R identical with D/F. Data collapse of scaled roughness distributions confirms that this model is in the class of the fourth-order nonlinear growth equation by Villain, Lai, and Das Sarma (VLDS) with negligible finite-size effects, while estimates of scaling exponents show some discrepancies. This result is consistent with the prediction of a recent renormalization group approach and improves previous numerical works on related models. The roughness follows dynamic scaling as W=Lalpha/R1/2f(xi/L), with correlation length xi=(Rt)1/z, where z is the dynamic exponent. We also propose a limited mobility (LM) model where the incident atom executes up to G steps before a new atom is adsorbed, and irreversibly aggregates at step edges. This model is also shown to belong to the VLDS class. The size of the plateaus in the film surface increases as G1/2 and the lateral correlation scales as G1/2t1/z. The time evolution of the roughness reproduces that of the CD model if an equivalent parameter G approximately R2/z is chosen. This suggests the possibility of using LM models with tunable diffusion length to simulate processes with simultaneous diffusion of many atoms. A scaling approach is used to justify exponent values and dynamic relations for both models, including the significant decrease of surface roughness as R or G increases.

Multifractal behavior of the surfaces evolved with surface relaxation

A discrete model exhibiting conserved dynamics with nonconserved noise involving particles of different nature, termed as linear and nonlinear, is proposed here. The morphology of the surface has been studied with different abundances of these particles. The saturated surface, slowly evolved from a lower contribution of nonlinear particles to a higher contribution of nonlinear particles, splits into four distinct scaling regimes with three crossover lengths. Each regime is characterized by different scaling property. It is shown that when the contribution of the nonlinear particles crosses a critical value, the surface morphology shows a linear-nonlinear "phase transition." The roughness exponent in a nonlinear regime is well compared with that of the continuum nonlinear equation in a molecular beam epitaxy (MBE) class as well as a MBE motivated discrete model.

Maximal- and minimal-height distributions of fluctuating interfaces

Maximal- and minimal-height distributions (MAHD, MIHD) of two-dimensional interfaces grown with the nonlinear equations of Kardar-Parisi-Zhang (KPZ, second order) and of Villain-Lai-Das Sarma (VLDS, fourth order) are shown to be different. Two universal curves may be MAHD or MIHD of each class depending on the sign of the relevant nonlinear term, which is confirmed by results of several lattice models in the KPZ and VLDS classes. The difference between MAHD and MIDH is connected with the asymmetry of the local height distribution. A simple, exactly solvable deposition-erosion model is introduced to illustrate this feature. The average extremal heights scale with the same exponent of the average roughness. In contrast to other correlated systems, generalized Gumbel distributions do not fit those MAHD and MIHD, nor those of Edwards-Wilkinson growth.

Finite-size effects in roughness distribution scaling

We study numerically finite-size corrections in scaling relations for roughness distributions of various interface growth models. The most common relation, which considers the average roughness <w(2)> as a scaling factor, is not obeyed in the steady states of a group of ballisticlike models in 2+1 dimensions, even when very large system sizes are considered. On the other hand, good collapse of the same data is obtained with a scaling relation that involves the root mean square fluctuation of the roughness, which can be explained by finite-size effects on second moments of the scaling functions. We also obtain data collapse with an alternative scaling relation that accounts for the effect of the intrinsic width, which is a constant correction term previously proposed for the scaling of <w(2)> . This illustrates how finite-size corrections can be obtained from roughness distributions scaling. However, we discard the usual interpretation that the intrinsic width is a consequence of high surface steps by analyzing data of restricted solid-on-solid models with various maximal height differences between neighboring columns. We also observe that large finite-size corrections in the roughness distributions are usually accompanied by huge corrections in height distributions and average local slopes, as well as in estimates of scaling exponents. The molecular-beam epitaxy model of Das Sarma and Tamborenea in 1+1 dimensions is a case example in which none of the proposed scaling relations work properly, while the other measured quantities do not converge to the expected asymptotic values. Thus although roughness distributions are clearly better than other quantities to determine the universality class of a growing system, it is not the final solution for this task.

Scaling in the crossover from random to correlated growth

In systems where deposition rates are high compared to diffusion, desorption, and other mechanisms that generate correlations, a crossover from random to correlated growth of surface roughness is expected at a characteristic time t0. This crossover is analyzed in lattice models via scaling arguments, with support from simulation results presented here and in other works. We argue that the amplitudes of the saturation roughness and of the saturation time t(x) scale as t0(1/2) and t0, respectively. For models with lateral aggregation, which typically are in the Kardar-Parisi-Zhang (KPZ) class, we show that t0 approximately p(-1), where p is the probability of the correlated aggregation mechanism to take place. However, t0 approximately p(-2) is obtained in solid-on-solid models with single-particle deposition attempts. This group includes models in various universality classes, with numerical examples being provided in the Edwards-Wilkinson (EW), KPZ, and Villain-Lai-Das Sarma (nonlinear molecular-beam epitaxy) classes. Most applications are for two-component models in which random deposition, with probability 1-p, competes with a correlated aggregation process with probability p. However, our approach can be extended to other systems with the same crossover, such as the generalized restricted solid-on-solid model with maximum height difference S, for large S. Moreover, the scaling approach applies to all dimensions. In the particular case of one-dimensional KPZ processes with this crossover, we show that t0 approximately nu(-1) and nu approximately lambda(2/3), where nu and lambda are the coefficients of the linear and nonlinear terms of the associated KPZ equations. The applicability of previous results to models in the EW and KPZ classes is discussed.

Numerical study of roughness distributions in nonlinear models of interface growth

We analyze the shapes of roughness distributions of discrete models in the Kardar, Parisi, and Zhang (KPZ) and in the Villain, Lai, and Das Sarma (VLDS) classes of interface growth, in one and two dimensions. Three KPZ models in d=2 confirm the expected scaling of the distribution and show a stretched exponential tail approximately as exp(-x0.8), with a significant asymmetry near the maximum. Conserved restricted solid-on-solid models belonging to the VLDS class were simulated in d=1 and d=2. The tail in d=1 has the form exp(-x2) and, in d=2, has a simple exponential decay, but is quantitatively different from the distribution of the linear fourth-order (Mullins-Herring) theory. It is not possible to fit any of the above distributions to those of 1/f(alpha) noise interfaces, in contrast with recently studied models with depinning transitions.