Height fluctuations in homoepitaxial thin film growth: A numerical study

We report on the investigation of height distributions (HDs) and spatial covariances of two-dimensional surfaces obtained from extensive numerical simulations of the celebrated Clarke-Vvedensky (CV) model for homoepitaxial thin film growth. In this model, the effect of temperature, deposition flux, and strengths of atom-atom interactions are encoded in two parameters: the diffusion to deposition ratio $R=D/F$ and $\varepsilon$, which is related to the probability of an adatom"breaking"a lateral bond. We demonstrate that the HDs present a strong dependence on both $R$ and $\varepsilon$, and even after the deposition of $10^5$ monolayers (MLs) they are still far from the asymptotics in some cases. For instance, the temporal evolution of the HDs' skewness (kurtosis) displays a pronounced minimum (maximum), for small $R$ and $\varepsilon$, and only at long times it passes to increase (decrease) toward its asymptotic value. However, it is hard to determine whether they converge to a single value or different nonuniversal ones. For large $R$ and/or $\varepsilon$, on the other hand, these quantities clearly converge to the values expected for the Villain-Lai-Das Sarma (VLDS) universality class. A similar behavior is observed in the spatial covariances, but with weaker finite-time effects, so that rescaled curves of them collapse quite well with the one for the VLDS class at long times. Simulations of a model with limited mobility of particles, which captures some essential features of the CV model in the limit of irreversible aggregation ($\varepsilon=0$), reveal a similar scenario. Overall, these results point out that the study of fluctuations in homoepitaxial thin films' surfaces can be a very difficult task and shall be performed very carefully, once typical experimental films have $\lesssim 10^4$ MLs, so that their HDs and covariances can be in the realm of transient regimes.


I. INTRODUCTION
The list of applications of thin film deposition is very vast, underlying a number of recent important technological advances [1][2][3]. In most of the techniques and relevant conditions for thin film production, molecules coming from a vapor or a liquid phase adsorb onto the substrate or film surface (with a flux F per adsorption site) and then passes to diffuse at the surface. This thermally activated diffusion process may depend on the local environment of the adatoms, on the diffusion mechanism (e.g, by hopping or by exchange), on how adatoms interact at the surface, and so on [2,3]. However, in most systems, it can be realistically modeled by simply considering the substrate and film surfaces as lattices of adsorption sites and diffusion occurring by the hopping of adatoms to nearest neighbor (NN) sites [4].
In the case of homoepitaxial growth -which is the one considered in this work -, no distinction shall exist between the diffusion on the substrate or film surfaces. The hopping rate H is expected to follow an Arrhenius form H = νe −E/(k B T ) , where ν can be considered as a constant frequency and E is the activation barrier for diffusion. Usually, in the modeling of these processes, E is separated into a site-independent terrace diffusion barrier (E d ), and a contribution (E n ) which shall depends on the local environment of the adatom, i.e, the height configuration around the actual site and the tar-get ones [3,4], as well as on the materials involved. However, quite often E n is considered as depending only on the initial neighborhood of the adatom, through a bondcounting approximation. Namely, if n is the number of intra-layer NN atoms of the adatom before hopping, so E n = nE N N , where E N N > 0 is the NN bond energy. This initial value approximation, usually referred to in the literature as the Clarke-Vvedensky (CV) [5] model, obeys detailed balance and has become very popular as a generic and somewhat realistic model for homoepitaxial thin film growth, especially for molecular beam epitaxy (MBE) [3,6]. In fact, the CV model has been widely employed in the modeling of island nucleation and growth during submonolayer deposition [7][8][9][10], as well as in the study of multilayer growth [11,12]. Furthermore, different modifications of the CV model have been considered as, e.g, the addition of an Ehrlich-Schwoebel (ES) [13,14] step-edge barrier, which is key to explain mound formation [4,[15][16][17].
The kinetic roughening process of surfaces simulated with the CV model (and those from homoepitaxial MBE growth in general) is expected to belong to the Villain-Lai-Das Sarma (VLDS) [18,19] universality class, aka nonlinear MBE class [6]. Namely, the films' height field h( x, t) is expected to evolve in the long time limit according to the VLDS equation [18,19] ∂h( x, t) ∂t which describes nonlinear growth processes dominated by surface diffusion, as is the case in CV and typical MBE films. In such equation, ν 4 and λ 4 are phenomenological parameters, while η is a Gaussian white noise. Some renormalization approaches have indeed demonstrated that the dynamic scaling of surfaces from CV model follows VLDS class [20]. Namely, their squared width w 2 scales in time, w 2 ∼ t 2β , and with the size l of the observation window, w 2 ∼ l 2α (for a given time), with the same exponents α and β from the VLDS equation [21]. From a numerical side, however, the situation is more controversial. While some Monte Carlo simulation studies of CV-related models on two-dimensional (2D) substrates -which is the case relevant for thin film deposition and the single one discussed here -indicated nonuniversal temperature-dependent exponents [22], other works have confirmed VLDS scaling for a broad range of model parameters [17,[23][24][25].
A number of recent works have demonstrated, notwithstanding, that universality in surface growth goes far beyond the scaling exponents. Importantly, the (1point) height distributions (HDs) and (2-point) spatial covariances, which are central quantities determining the statistics of fluctuating surfaces, are believed now to display universal behaviors, beyond an interesting dependence on geometry. In fact, this has been widely established for systems belonging to the Kardar-Parisi-Zhang (KPZ) [26] class in both 1D [27][28][29][30][31][32][33][34][35][36][37] and 2D [38][39][40][41] substrates, whereas a single study exists indicating the same for the VLDS class [42], as well as another one for linear universality classes [43]. In face of this, HDs and spatial covariances have been used in several recent works as a tool for determining the universality class of real thin film surfaces [44][45][46][47][48][49]. Despite this appealing application, as far as we know, HDs and covariances have never been investigated for realistic models for thin film growth. For instance, all results for the VLDS class reported in [42] were obtained from simulations of "toymodels" (with limited mobility of adatoms). Therefore, the aim of the present work is two-fold: (i ) try to confirm the universality of HDs and covariances for the VLDS class in 2D through the analysis of these quantities for the CV model; and, conversely, (ii ) investigate the effects of temperature (T ), deposition flux (F ) and bond energy strengths on such quantities. As will be demonstrated in the following, the CV HDs display strong and longliving corrections for low T (and/or large F ), turning almost impossible to firmly establish their universality class in feasible deposition times. At high T , however, they agree with that for other VLDS models. The covariances present mild finite-time corrections and always agree with the VLDS one at long time. We simulate also a simplified (limited mobility) model [24] which mimics some features of the CV model in the limit of irreversible aggregation, and the same behavior was observed there.
The remainder of the paper is organized as follows. In Sec. II the investigated models are defined, as well as the quantities used in the analyses of their surfaces. Results for the HDs and spatial covariances are presented in Secs. III and IV, respectively. Section V summarizes our final discussions and conclusions.

II. MODELS AND QUANTITIES OF INTEREST
The Clarke-Vvedensky (CV) [5] model is a version of the classical solid-on-solid (SOS) model by Gilmer and Bennema [50], for the regime of complete condensation, when particle desorption from the surface is negligible. The SOS condition implies that overhangs are not allowed at the surface, so that the films formed are compact. In the CV model, particles are deposited at random positions of the substrate with a homogeneous flux F per adsorption site. Here, we consider the substrate as a square lattice of lateral size L, with periodic boundary conditions, which is an approximation for films growing in (001) orientation. The growth is performed on flat substrates. Once deposited, an adsorbed particlethe adatom -passes to diffuse at the surface through the hopping to nearest-neighbor (NN) sites. The total hopping rate is given by H = νe −(E d +nE N N )/k B T , where n is the number of intra-terrace NN atoms of a given adatom before hopping. The frequency ν was considered as ν = k B T /(2π ) in the original CV formulation, as expected from transition theory, with being the Planck's constant, k B the Boltzmann's constant and T the temperature. However, it is quite common to consider ν simply as a constant (ν ∼ 10 12 − 10 13 Hz) [3,4] and we will adopt this simpler definition here. In order to model the growth of a specific system with the CV model, we should also specify reliable values for the activation energy barriers E d and E N N for such a system. In the present work, notwithstanding, we are interested in investigating general properties of the height fluctuations of CV films for broad ranges of F , T and energies. So, instead of dealing with all parameters of the model explicitly, it is convenient to rewrite the hopping rate as H = Dε n , where D = νe −E d /k B T is the hopping rate of free adatoms and ε = e −E N N /k B T is related to the probability of an adatom "breaking a bond" with a lateral NN one. Since the growth is controlled by the ratio H/F = (D/F )ε n , it is interesting to define also R = D/F , so that R and ε will be the model parameters.
We investigate also a limited mobility SOS model introduced by Aarão Reis [24], where again particles are randomly deposited onto a flat square lattice substrate, but only the freshly deposited adatom is allowed to diffuse at the surface. Such adatom can diffuse (by randomly hopping to NN sites, similarly to CV model) while it is free and once it arrives (by deposition or diffusion) at a site with at least one in-plane NN atom, it irreversibly aggregates there. If no one of such sites is found after G diffusion steps, the adatom stops diffusing and permanently aggregates at its final position. Therefore, G is the only parameter in this lateral aggregation of diffusing particles (LADP) model. As demonstrated in [24], by tuning G this model can yield surfaces with the same roughness properties of the CV model with irreversible aggregation, i.e., with ε = 0. Moreover, strong numerical evidence that this model belongs to the VLDS class has been provided in [24,51].
In both CV and LADP models, the time unity will be defined as the deposition of one monolayer (ML) of particles. Here, times up to t = 100000 will be analyzed, corresponding to the deposition of 10 5 MLs. We remark that this number of MLs is at least one order of magnitude larger than those considered in typical experiments, as well as in recent numerical works on these and related models [24,25,51]. Results for lateral substrate sizes L = 1024 (and L = 512 in some few cases) will be presented below. For each set of model parameters, at least one hundred of different films were grown, totaling 26×10 6 points in the statistics. Some less accurate simulations for L = 256 and L = 2048 were also performed to confirm that finite-size effects are negligible in our data.
According to the so-called "KPZ ansatz", the height at a given point of the films' surfaces, during the transient growth regime, is expected to evolves as [27,52] h where the growth velocity is always v ∞ = 1 for the models considered here, Γ is a model-dependent parameter setting the roughness amplitude, and β is the universal growth exponent. Moreover, χ is a random variable fluctuating according to a probability density function P (χ), i.e., the underlying height distribution (HD), which is expected to be universal, but dependent on the surface geometry. The films analyzed here are flat and, hence, we will compare our results with previous ones for this geometry. The nth centered moment h n m of the HDs is defined as h n m = h −h n , where the bars indicate averages over the heights of each film surface, at a given time, while · · · denotes the average over different films. In order to characterize the HDs (and other distributions as well) one usually analyzes adimensional ratios of their moments, being the skewness S = h 3 m / h 2 3/2 m = χ 3 m / χ 2 3/2 m and the kurto- The two-point spatial statistics of the films' surfaces can be investigated through the correlator whereh ≡ h − h . From dynamic scaling theory, one expects that C S (r, t) w 2 F (r/ξ), where ξ is the correlation length parallel to surface, with F (x) being an, in principle, universal and geometry-dependent spatial covariance. Once it is usually hard to obtain accurate estimates of ξ and the covariances obtained here for the CV model display a pronounced minimum at r = r min (due to a modulated oscillatory decreasing behavior, analogous to that found in other models with dynamics dominated by surface diffusion [42,43]), we will analyze the rescaled covariances by plotting C S (r, t)/w 2 against r/r min . In this way, all rescaled curves coincide at r/r min = 0 (where C S /w 2 = 1) and have minima at r/r min = 1, so that they shall collapse whenever they have the same shape.

III. RESULTS FOR THE HEIGHT DISTRIBUTIONS
In this section, we focus on the behavior of the height distributions (HDs) during the transient growth regime. Initially, the case of irreversible aggregation is discussed, which is followed by an analysis of the full CV model.

A. CV model with irreversible aggregation
Before analyzing the full CV model, it is interesting to start considering the limiting case ε = 0. This corresponds to strong intra-terrace adatom-adatom interac- where effectively only free adatoms can move and there is only a single parameter, R, to be considered. Although this seems an oversimplification, we notice that this "irreversible aggregation CV" (IACV) model has already been investigated for multilayer growth [24,53]. Moreover, it is a very common practice to assume that aggregation is irreversible in studies of island nucleation and growth during the submonolayer regime, where the model considered here is a kind of multilayer version of the situation with "critical nucleus" i * = 1 [3,4,[54][55][56][57].
Let us start investigating the effect of R on the amplitude Γ in the "KPZ ansatz" (Eq. 2). From such equation, one expects that h n m = (Γt) nβ χ n m + . . ., for n 2, and thence the quantity g n ≡ h n m /t nβ should converge to a constant g n → Γ nβ χ n m at long times. As demonstrated by Aarão Reis [24], the surface roughness W = h 2 m for the IACV model follows the scaling W = (L α /R 1/2 )f (ξ/L), with the lateral correlation length ξ = (Rt) 1/z , where z = α/β is the dynamic exponent. Since the scaling function f (x) behaves as f (x) ∼ x α in the growth regime, we shall have h 2 m At 2β /R 1−2β , with β 0.1975 being the exponent of the 2D VLDS class [21]. In fact, rescaled curves of h 2 m R 1−2β /t 2β versus t for not so large R tend to a plateau, where they collapse reasonably well, as shows Fig. 1(a). This confirms the scaling behavior and allow us to estimate the amplitude A 12(2) (see the extrapolation in the insertion of Fig. 1). Moreover, this shows that Γ 2β χ 2 m = A/R 1−2β and, by assuming that χ 2 m is universal (i.e., independent of R) one has Γ ∼ R −γ , with the exponent γ = 1/(2β) − 1 ≈ 1.5316. If this is the case, rescaled curves of h n m [R γ /t] nβ for the higher moments should also converge to constant values as t → ∞. An example of this is shown in Fig.  1(b), for the fourth moment (n = 4), where one sees that curves for different R's do not collapse so well as those for n = 2 [in Fig. 1(a)]. It is noteworthy, however, that such rescaled curves do not present clear plateaus even for the smaller R's, showing that the higher moments have stronger finite-time corrections. This possibly explains the absence of a good data collapse and suggests that the plateau will occur for B ≡ h 4 m R 2−4β /t 4β 500, when t → ∞. We notice that this is consistent with the 2D VLDS class, since B/A 2 = h 4 m / h 2 2 m = K + 3 and a kurtosis K ≈ 0 was found in [42] for the VLDS HDs, so that B = 3A 2 is expected to be found approximately in the range [300, 600]. For large R, we do not observe any plateau or collapse in Fig. 1, certainly because the moments are still far from the VLDS scaling regime [ h n m ∼ t 0.1975n ]. The temporal evolution of the skewness S and kurtosis K of the HDs are depicted in Figs. 2(a) and 2(b), respectively. As one can see, these quantities present very strong finite-time effects for small R, with very pronounced minima (maxima) in S (K), and they seem to be far from their asymptotics even after the deposition of 10 5 MLs. For R 10 5 , however, such non-monotonic convergence gives place to an initial oscillatory behavior in S and K, due to a layer-by-layer growth at (relatively) short times, and their values converge to the ones numerically estimated for other models belonging to 2D VLDS class in flat geometry: |S V LDS | ≈ 0.13 and K V LDS ≈ 0.00 [42]. This suggests that for any R the HDs shall agree with the VLDS one, but huge deposition times are needed to observe this when R is not so large. In fact, the ratios for R = 10 4 also agree reasonably well with S V LDS and K V LDS at long times (see Fig. 2).
As demonstrated in [24], by setting G ≈ 0.28R 0.58 in the LADP model, it is capable of reproducing the roughness behavior of the IACV model with parameter R. Namely, the curves of W × t for both models collapse. This leads us to inquire whether the same thing happens with the higher moments and their ratios and, quite interestingly, the answer is positive. In fact, rescaled curves of h 2 m and h 4 m for the LADP model with G = 1, 4, 15, 59 and 222 -which would mimic the IACV surfaces for R = 10, 10 2 , 10 3 , 10 4 and 10 5 , respectively -are also shown in Figs. 1(a) and 1(b) and they remarkably kurtosis K for the IACV (circles) and the LADP (plus symbols) models, for the same parameters and color scheme from Fig. 1. The dashed lines represent the universal values of S and K for the 2D VLDS class estimated in Ref. [42]. The arrows indicate the times at which the covariances displayed in Fig. 6(a) were measured. collapse with the IACV ones. A similar agreement is observed also in S and K for both models, for the same parameters, as seen in Fig. 2. This demonstrates that not only the variance of the HDs coincide, but the full HDs are approximately the same for both models, at a given time, for that choice of parameters. Therefore, since simulations of the LADP model are much faster than those for the IACV one, we can use it to investigate the HDs' behavior for larger R. For instance, figures 1 and 2 also show the moments and ratios, respectively, for the LADP model with G = 846, corresponding to the IACV model with R = 10 6 . We remark that this would be extremely hard to simulate with IACV model for the times, sizes and number of samples we are considering. These results are in consonance with the ones for smaller R's and, particularly, indicate that for R > 10 5 the skewness passes to converge to S V LDS from above. The finite-time effect in K is also stronger than that observed in R = 10 5 . This is indeed expected because by increasing R the duration of the transient layer-by-layer regime also increases. Therefore, although the LADP model does not display such regime, it seems that its surfaces take similar times as the IACV ones to attain the HDs' asymptotic regime.

B. CV model with reversible aggregation
Now, we consider the original CV model, with a nonnull probability ε for an adatom detaches from a lateral neighbor. Once again, it is interesting to analyze how the nonuniversal parameter Γ in Eq. 2 depends on R and ε. In Ref. [25], numerical evidence was provided that during the growth regime the surface roughness behaves as where γ = 1/(2β) − 1 as in the previous subsection. This means that plots of h n m [R γ (0.025 + ε)/t] nβ versus t should converge to a single plateau as t → ∞, what is somewhat confirmed here, as shows Fig. 3. For the second moment a reasonable collapse is found for long In comparison with the amplitude A above, we have A /A = 0.025 2β ≈ 0.23, which is consistent with our numerical estimate A /A = 0.21 (7). From the discussion in the previous subsection, for the 2D VLDS class, a plateau at B = 3A 2 approximately in the range [10,30] is expected in the fourth moment, which is indeed observed in Fig. 3(b). This is in agreement with the results above for the IACV model, indicating that the finitetime effects are stronger in the higher-order moments. It is important to remark that the simple dependence of Γ on R and ε presented here, can also present corrections. For instance, rescaled curves for R = 10 [not shown] do not collapse so well with those depicted in Fig.  3. Moreover, in general, the collapses can be considerably improved by including (somewhat arbitrary) logarithmic corrections in Γ [not shown]. Anyhow, a firm conclusion about the existence and form of such corrections would require simulations for larger R's and then larger substrate sizes and much longer deposition times than those considered here, which are not feasible with our current computer resources.

Figures 4(a) and 4(b) present the variation in time of
the HDs' skewness S and kurtosis K for ε = 0.01 and several values of R, where a behavior analogous to that found for the IACV model (Fig. 2) is observed. Namely, for small R severe finite-time effects appear, hallmarked by minima (maxima) in S (K), and it is hard to determine whether these ratios converge to the same or different asymptotic values. For the larger R's, on the other hand, there is a faster convergence and the asymptotic values agree quite well with those for the VLDS class. Since the ratio R = D/F determines how much the free adatoms diffuse at the surface before the deposition of a new one, these results and those from the previous subsection demonstrate that when such diffusivity is low the HDs suffer with strong corrections, which decreases as R increases. Note that the low adatom diffusivity for small R yields the nucleation of a large number of 2D islands during the submonolayer regime, over which the subsequent 3D islands grow. Once an adatom is deposited over one of such 3D islands, it has a small probability of escaping from it, for small R, so that such islands tend to grow fast vertically. This yields deep valleys between them in the films' surfaces, which certainly explains the rising of HDs with large negative skewness at (relatively) short times. Since the 3D structures tend to form plateaus at their centers, this explains also the large positive kurtosis. As such structures coalesce, however, the deep valleys tend to diminish and then |S| and K start decreasing. By increasing R, fewer islands are nucleated in the submonolayer regime, forming fewer valleys, or even a layer-bylayer growth appears at short times for large enough R. In any case, this prevents the advent of strong finite-time corrections in the HDs. This scenario is qualitatively confirmed in Figs. 5(a) and 5(b), where the temporal evolution of characteristic 1D cross-sections of films' surfaces are compared for R = 10 and R = 1000, respectively. In the former case, we find that surfaces are indeed featured by deep valleys at short times, which tends to decrease for very large t, while much more smooth morphologies are found for R = 1000.
According to the reasoning above, we might expect also to find a decrease in the finite-time effects by increasing ε, for a given R, because a large detachment rate of adatoms from steps also contributes to prevent the formation of deep valleys in the films' surfaces. This is indeed observed in Fig. 5(c), where typical 1D profiles for CV films for R = 10 are shown for the same deposition time (1000 MLs) and different values of ε. A quantitative confirmation of this is given in Figs. 4(c) and 4(d), where S and K as function of time are depicted for CV films with R = 10 3 and several values of ε. There, one sees that the minima (maxima) in S (K) curves decrease as ε increases and for ε 0.05 they converge to the VLDS ratios. This strongly suggests that for the smaller ε's they shall also converge to there, but taking much longer times.
Results for R = 10, however, indicates that this is not the case for large values of ε. In fact, curves of S for ε 0.10 converge to nonuniversal ε-dependent values at long times [see 4(e)]. We remark that this is somewhat expected because the case ε = 1, where adatoms do not interact with their NN ones (since E N N = 0), corresponds to a kind random deposition. As a matter of fact, although adatoms are still diffusing, they only stop moving when they are earthed by another adatom, due to a random deposition or a random diffusion towards it. These random processes alone do not generate correlations in the system and, thus, the surface roughness increases as W ∼ t 1/2 and the HDs are Gaussian. This is indeed confirmed in Fig. 4 (c), where the skewness for ε = 1 converges to S = 0. Moreover, the kurtosis for this parameter also converges to K = 0 [see Fig. 4(f)]. The same is true for larger values of R. We can see in Fig.  4(c) that the curve of S for R = 1000 and = 0.2 seems do not converge to S V LDS or it will take a very long time to arrive at that. In all panels, the dashed line is the universal covariance curve for the 2D VLDS class, estimated in [42].

IV. RESULTS FOR THE SPATIAL COVARIANCES
Now, we discuss the two-point spatial covariance function (Eq. 3) of CV films' surfaces. Once again, let us start with the case of irreversible aggregation. Figure  6(a) shows rescaled covariance curves for the IACV model with R = 1000 and the LADP one with G = 15, for the times indicated by the arrows in Fig. 2(a). We observe that for short times, when the HDs are quite different from the VLDS one -i.e., the values of S and K are far from S V LDS and K V LDS [see Figs. 2(a) and (b)]the covariances are also different from the VLDS curve, which was numerically estimated by us in [42]. However, as time increases, the covariances converge towards the VLDS one, presenting a reasonable collapse with it for the longest time considered here. This makes clear that the same sort of finite-time effects observed in the HDs are also present in the covariances. Figure 6(b) presents data for the IACV model with R = 10 4 and its LADP counterpart -for which the HDs agree well with the VLDS one already for t 10 4 -and, in fact, a quite good collapse of the rescaled covariances with the VLDS one is observed for t = 10 4 . Therefore, similarly to what happens in the HDs, as the adatom diffusivity increases the finite-time corrections decrease, uncovering the VLDS universality in the system. It is important to emphasize also the good agreement between the IACV and LADP covariances for the parameters G ≈ 0.28R 0. 58 . This indicates that the entire statistics of IACV films' surfaces are captured by the simplified LADP model.
The results for the CV model are shown in Figs. 6(c) and 6(d) for ε = 0.01 and several R's, and R = 10 3 and several ε's, respectively. In this case, the time dependence is similar to that just discussed for the irreversible aggregation [ Fig. 6(a)], so, only results for the longest simulation time are presented in Figs. 6(c) and 6(d). There, one can see that in general the curves for the CV model with very different parameters collapse quite well among them, as well as with the VLDS curve. Especially for R = 10, although a slight deviation is observed in Fig.  6(c) for small r/r m , the overall behavior is quite close to the VLDS one. This demonstrates that the finite-time corrections in the covariances are milder than those found in the HDs. Furthermore, the agreement with the VLDS covariances provides a strong confirmation that the CV model belongs to VLDS class.

V. CONCLUSIONS
We have presented an extensive analysis of the 1-point and 2-point statistics of surfaces of the CV model for homoepitaxial thin film growth, for broad ranges of model parameters (R ∈ [10, 10 5 ] and ε ∈ [0, 1]), covering a wide range of temperature (and/or deposition flux) and energy strengths. For instance, if one assumes that ν = 10 13 s −1 and E d = 1 eV, which are close to the values expected in the growth of several semiconductor, metal and organic films [3,4], for a flux F = 1 ML/s we have simulated temperatures approximately in the interval T = 420 − 630 K. Conversely, for a typical temperature value, e.g., T = 420 K, one has the flux in the interval F = 0.0001 − 1.0 ML/s. In any case, these parameters are consistent with those commonly used in actual thin film deposition experiments [3,4].
We have also simulated the simplified LADP model [24], where the mobility of adatoms is limited to the freshly deposited one and found that, by appropriately tuning the parameter in such model, it is able to produce surfaces with the same HDs and covariances of the CV model with irreversible aggregation (IACV). This is a rather interesting finding, once the already known fact that both models have the same (actually, very similar) roughness evolution, as demonstrated in [24], does not necessarily imply that their surfaces would have the same 1-point and 2-point statistics. Our results, however, demonstrate that this is indeed the case. As pointed out above and stressed recently in Refs. [53,58], the existence of this kind of simplified model reproducing the surface features of a more realistic and complex one (the IACV model here) is very important since this allows us to investigate regimes of the latter model, which would not be computationally accessible in a feasible amount of time, through simulations of the former one.
Most of our results indicate that the HDs and covariances asymptotically agree with those obtained in Ref. [42] in simulations of simplified models belonging to the 2D VLDS class. Substantially, this provides a strong confirmation that the CV model indeed belongs to the VLDS class and, conversely, confirms the universality of the HD and covariance for this class.
It turns out, however, that severe finite-time effects are observed in these quantities, especially in the HDs, for small R (corresponding to low T and/or high F ). Namely, when the adatom mobility is low when compared with deposition; what leads to the formation of deep valleys in the films' surfaces, which yield the deviations. For a given value of R, the atom-atom energy E N N plays also an important role in this matter, enhancing the finite-time effects when it is large (i.e., when ε is small). Interestingly, this behavior is contrary to the one found in the surface roughness, W , where VLDS scaling is observed already at short times for small R, but the time it takes to appear increases with R, due to the initial layer-by-layer regime which appears for large R. For this reason, a clear scaling W ∼ t 0.1975 is not observed in simulations of the CV model for large R, even for the deposition of 10 5 MLs, as shown in Figs. 1(a) and 3(a). Similar conclusions have been reported in [24,25,51]. Hence, HDs and covariances are complementary measures to W in such systems, with the former (later) working better for large (small) R's. For this reason, it is very important to investigate all these quantities together in order to determine the universality class of a given growing film. Moreover, a careful finite-time analysis is imperative in face of the strong effects observed here, which may be present also in other realistic models, as well as in experiments on homoepitaxial thin film growth.