Temperature Effects in the Initial Stages of Heteroepitaxial Film Growth

Kinetic Monte Carlo simulations of a model of thin film heteroepitaxy are performed to investigate the effects of the deposition temperature in the initial growth stages. Broad ranges of the rates of surface processes are used to model materials with several activation energies and several temperature changes, in conditions of larger diffusivity on the substrate in comparison with other film layers. When films with the same coverage are compared, the roughness increases with the deposition temperature in the regimes of island growth, coalescence, and initial formation of the continuous films. Concomitantly, the position of the minimum of the autocorrelation function is displaced to larger sizes. These apparently universal trends are consequences of the formation of wider and taller islands, and are observed with or without Ehrlich-Schwöebel barriers for adatom diffusion at step edges. The roughness increase with temperature qualitatively matches the observations of recent works on the deposition of inorganic and organic materials. In thicker films, simulations with some parameter sets show the decrease of roughness with temperature. In these cases, a re-entrance of roughness may be observed in the initial formation of the continuous films.


Introduction
Models of thin solid film deposition can be used to determine the main mechanisms that control their morphology and, consequently, help to design the most suitable conditions to obtain materials with the desired features [1][2][3][4][5]. Some modeling approaches are based on the comparison of quantities obtained from microscopy images, such as the surface roughness, which is usually defined as the root mean square deviation of the local heights, or the spatial correlations of those heights [2,5]. In simple growth processes, these quantities follow universal scaling laws that are related to the symmetries of stochastic growth equations [6,7]. However, this usually fails to describe the initial stages of deposition, the evolution of mounded deposits, and other small scale features.
In vapor deposition processes, the substrate temperature is one of the main parameters to determine the film properties [1]. In models of homoepitaxial deposition, the temperature effects on the roughness were widely studied [5]. If Ehrlich-Schwöebel (ES) [8] barriers for step edge descent are neglected, simulations show a decrease of surface roughness with temperature [9][10][11][12]. However, if ES barriers and other mechanisms such as downward funelling are included in the models, different features may arise [5]. An example is the nonmonotonic temperature dependence of the roughness, which is experimentally observed in the deposition of metals, such as Ag/Ag(100) and Cu/Cu(100) [13][14][15][16][17]. In this so-called re-entrant growth, the temperature increase leads to the roughness increase up to a certain temperature, but subsequently leads to a decrease in the roughness. Simulations of related models also predicts a roughness increase with the temperature or other types of re-entrant behavior [18,19].
In heteroepitaxial deposition, the nucleation and growth of isolated islands is usually observed, in contrast with the initial layer-by-layer growth of homoepitaxial deposition.
The increase of the surface roughness with the temperature was recently observed in the initial growth stages of some materials. In organic film growth, two examples are ≈1-12 nm cobalt phthalocyanine films on a SiO 2 /Si(001) substrate [20] and ≈5-15 nm MAPbBr 3 films (MA= CH 3 NH 3 ) on SiO 2 and ITO substrates [21]. The same is observed in hot wall CdTe deposition on Si(001) substrates when the coverages are of some nanometers and the height profiles are inherited from the island configurations [22].
Atomistic models were also used to describe heteroepitaxial deposition of several materials [23][24][25][26][27][28][29]. Regarding temperature effects, we highlight recent simulations with a focus on the height-to-radius aspect ratio of Ag islands [26,27]. However, as far as we know, no previous model has explained the effects on the surface roughness as observed in the experimental studies mentioned above.
The aim of this work is to investigate the effects of the deposition temperature on the morphology of heteroepitaxial films in their initial growth stages. We perform kinetic Monte Carlo (kMC) simulations of a recently proposed model that represents its main microscopic processes [30], whose rates depend on the deposited material, the substrate, and the temperature. The kMC approach is useful to connect such features at atomic/molecular scale with the morphological features at much larger scales (several nanometers or, possibly, micrometers), although with simplifications of the microscopic interactions. The general trend observed here is the increase of roughness with the deposition temperature when films with the same thickness are compared, in parallel with recent experiments. Moreover, the island width increases in these conditions, which is quantitatively confirmed by the changes in the height autocorrelation function.

Deposition Model
The deposit had a simple cubic structure, which is computationally favorable for the simulation of wide and thick deposits. The lattice constant was the length unit of the model and the lateral size was set as L = 1024, which is large enough to avoid finite size effects. The flat substrate was located at z = 0, the deposit was formed at z ≥ 1, and periodic boundary conditions were considered in the x and y directions. Each atom or molecule was assumed to occupy a single lattice site; for simplicity, hereafter, we only use the term adatom for the deposited species. Solid-on-solid conditions were considered, so the deposits had no pores. The set of adatoms with the same (x, y) position was termed a column of the deposit. The height variable h(x, y) was the value of z of the topmost adatom in that column.
The atomic flux was collimated in the −z direction with a rate of F atoms per substrate site per unit time. A column (x, y) was randomly chosen for each incident atom, which adsorbs as it lands at the top of that column.
Desorption was neglected in this model. In general, this was not restrictive in homoepitaxial growth modeling [5]. In heteroepitaxial growth, this condition was fulfilled if the desorption rate was much smaller than the incident flux F, which depends on the control of growth conditions and on the material properties.
Only the adatoms at the top of the L 2 columns and z ≥ 1 executed surface diffusion (substrate atoms are immobile). When an adatom attempted to hop, the target direction was randomly chosen among the four NN columns, ±x or ±y.
The rates of intralayer adatom hops (i.e., hops without change in the z position) depended only on the number n of lateral nearest neighbors (NNs). The numbers of hops per unit time are: where D S and D A are the hopping rates of adatoms without lateral NNs on the substrate and on the deposit, respectively, which are expected to have activated Arrhenius forms: where ν is a frequency, E S > 0 and E A > 0 are activation energies, k B is the Boltzmann constant, and T is the substrate temperature. The same is expected for the factor , which is interpreted as a detachment probability per lateral NN: where E B > 0 is a bond energy that represents the interaction with each lateral NN. The rates of interlayer adatom hops are reduced in comparison with those of Equation (1). This is accomplished by multiplying those rates by a factor (probability) denoted as P hop . Due to the solid-on-solid constraint, an adatom hop in the model is always performed by removing it from the top of the current column and inserting it at the top of an NN column, as shown in Figure 1a (left). If the height difference between the columns is ∆h NN > 1, the model hop represents a sequence of the two processes shown in Figure 1a: (i) the step edge crossing (from column 1 to column 2), whose activation energy (ES barrier) is denoted as E ES ; (ii) the random motion along a vertical terrace (in column 2). During that motion, it is possible that the adatom returns to the initial position in column 1 before it reaches the top of column 2 (this return is not shown in Figure 1a). Thus, P hop has a contribution of E ES and a contribution of a possible return. Accounting for these contributions, Ref. [18] obtained the overall probability of executing the hop as: , with probability 1 − P hop , the adatom remains at the current position. The possible hops of some adatoms are illustrated in Figure 1b, with the corresponding rates and, in the case of interlayer hops, the reducing factors P hop for those rates.

Model Justification and Dimensionless Parameters
The model studied here is an extension of widely used models of metal and semiconductor growth [5]. Similar mechanisms of adatom hopping were already proposed for deposition of molecular films, including cases of large molecules such as C 60 [31,32].
The model was introduced in Ref. [30] to describe the transition from island growth to continuous film formation in heteroepitaxial deposition. It was shown that the adatom mobility on the substrate and the ES barrier are the main factors that control the roughness and the correlations in that transition. The results were used to interpret atomic force microscopy data of CdTe deposits at constant substrate temperature; possible applications to other materials were also suggested [33,34]. However, the effects of increasing the deposition temperature were not addressed in that work.
The deposition conditions may be characterized by a set of four dimensionless parameters. The diffusion-to-deposition ratios on terraces are defined as: whereν = ν/F. They may be interpreted as the average numbers of hops of adatoms on terraces during the average time 1/F of deposition of one atomic layer. The other dimensionless parameters are the detachment probability and the probability P of crossing a monolayer step edge (∆h NN = 1 in Equation (4)). For a given deposited material and a given substrate, the set of activation energies is constant. If the atomic/molecular flux remains constant, the four dimensionless parameters R S , R A , , and P increase as the temperature increases.

Basic Quantities
The dimensionless film thickness is denoted as: This is the average number of deposited atoms per column at time t (number of deposited monolayers).
The substrate covering thickness θ c is defined as the first thickness in which 99% of the layer z = 1 is filled with adatoms. It corresponds to almost complete filling of the layer in contact with the substrate, so we expect that almost all islands have already coalesced at θ c and the continuous film will soon be formed. The value of θ c is expected to be roughly proportional to the elongation transition thickness, in which the substrate is predominantly covered by elongated structures [35].
The height fluctuations are characterized by the roughness: where the overbars denote a spatial average and the angular brackets denote an average over different configurations with a given thickness. The auto-correlation function measures the spatial correlations in the height fluctuations. At a distance s, it is defined as [36]: where the configurational average is taken over different initial positions r 0 , different orientations of s (directions x and y), and different deposits. W 2 is the square roughness defined in Equation (7) without taking the square root of the spatial averageh 2 . For any time and film thickness, Γ = 1 at s = 0. The position of the first minimum of the autocorrelation function is denoted as s M . In mounded surfaces, it is an approximation for the average distance between the top of the hills and the bottom of the valleys [36]; similar interpretation is possible in deposits with islands.
The quantities θ c , W, and s M are dimensionless and are given in units of the lattice constant throughout this work.

Simulation Parameters
We assume the following relations between the activation energies: Relation (a) implies that the adatom mobility on the substrate (z = 1) is higher than the mobility on the deposited material (z ≥ 2). In general, we choose parameters consistent with R S R A , which leads to the formation of islands with two or more layers in the beginning of the deposition, as shown in Refs. [30,37].
Relation (b) is obeyed for most metals and semiconductors, usually with E ES much smaller than both E A and E B [5]. However, recent simulation works considered E ES of the same order as or larger than E B [26,31]. For this reason, most of our simulations consider P > , but P ≈ is also considered in some parameter sets.
In previous kMC simulations of homoepitaxial or heteroepitaxial deposition, the frequency ν is set on the order 10 12 s −1 [5,23,24,27]. Fluxes on the order of 1-100 monolayers per second are used in the experimental works of interest here [20][21][22]. Most of our simulations are performed withν = ν/F between 10 10 and 10 12 , which is expected to represent a broad range of possible physical conditions.
In the simulation with a given set of physical parameters (activation energies and F), the temperature remains constant during the deposition. Table 1 shows the values of activation energies,ν, and the two temperatures in which the deposition was simulated, which are denoted as T 1 and T 2 . The deposition at the lowest temperature is always represented by a single letter (A to E). The deposition at the highest temperature is represented by a letter (A to E) followed by a number (1 to 3). Below, we justify the choice of these parameters sets. Table 1. Sets of activation energies and temperatures used in the simulations. Observe that deposition is simulated at constant temperature, T 1 or T 2 , not at varying temperatures. We first defined five initial deposition conditions, labeled from A to E, which correspond to systems without ES barriers (A,B) or with them (C,D,E). Each of these conditions has constant values of the dimensionless parameters R S , R A , , and P, which are the quantities directly used by the kMC algorithm. The values of the dimensionless parameters are listed in Table 2. However, each of these sets may correspond to several different combinations of the physical parameters, namely the activations energies, the temperature, and the ratio ν/F between hopping frequency and flux (Equations (3)-(5)). Thus, for each set (A to E), we chose up to three possible sets of physical parameters (shown in Table 1). In each case, a final deposition temperature was also chosen and used to determine the final values of the dimensionless parameters (labeled by a letter and a number), as shown in Table 2. The maximal film thickness simulated here was θ = 640, which exceeds the substrate covering thickness θ c in all cases. For each parameter set, 10 deposits were generated. The fluctuations of the average quantities among the different deposits of each set were small.
The simulations were implemented with the kinetic Monte Carlo algorithm described in detail in Ref. [37].

Growth with Negligible ES Barrier
Here, we consider the case P = 1 in Eqution (4). This is reasonable, for instance, to describe the hot wall deposition of CdTe on a flexible substrate [30]. Figure 2a,b shows top views and cross-sections of deposits grown with parameter sets A and A1, respectively. Here and in other deposition conditions, we observed three regimes. First, islands grow in width and height, but leaving large gaps between them; this is the island growth regime, which is clear at θ = 5 in Figure 2b. As the gaps between the islands are filled, they coalesce, forming large structures; this is the island coalescence regime, which is observed at θ = 5 in Figure 2a and θ = 20 in Figure 2b. Finally, a continuous film formed when the substrate was fully covered, which occured shortly after the substrate covering thickness θ c .
The temperature increase from A to A1 was represented by increasing R A by a factor of 15, whereas R S and increase by smaller factors, and P remains constant. The temperature change led to the formation of wider and taller islands, which delayed their coalescence and formation of the continuous film. Quantitatively, we obtained θ c ∼10 and θ c ∼40 in sets A and A1, respectively. After full filling of the substrate, the continuous films had smooth surfaces because, in the absence of step edge barriers, the mobile adatoms could easily reach high coordination points and aggregate there [11]. Similar results were obtained when deposits grown with parameters B and B1 were compared.  These results confirm that θ c increases with temperature, as suggested by the images in Figure 2. As the deposition began, W increased with thickness because the island heights were increasing and there were large areas of almost uncovered substrate. The maximal roughness was attained in a thickness very close to θ c , i.e., when most of the substrate was already covered. The existence of this maximal roughness is typical of deposition with negligible ES barriers, as shown in Ref. [30].  When deposits with the same thickness are compared, Figure 3a shows that the temperature increase leads to the increase of W in the regimes of island growth and coalescence. This is related to the formation of taller islands at higher temperatures. The same trend remains in the initial stages of continuous film formation. This occurs because the fluctuations in the average island heights also increase with the temperature, so longer times are necessary for the adatom diffusion to suppress those fluctuations after coalescence. Figure 3b shows the autocorrelation function Γ of deposits grown with parameter sets B and B1 for three thicknesses. When deposits of the same thickness were compared, the position s M of the minimum of Γ increased with temperature. This quantitatively confirms that the islands became wider.
Long after the formation of the continuous film, Figure 3a shows that the roughness in sets A and B exceeded those of A1 and B1, respectively. In these conditions, the memory on the initial island configuration was lost, so the relation between W and the temperature was the same as that shown in studies of homoepitaxial deposition [11,38]. In this regime, the larger values of s M at higher temperatures were related to longer lateral correlation lengths [11,38].

Growth with ES Barriers
Figure 4a-c shows top views and cross-sections of deposits grown with parameter sets D, D2, and D3, respectively. From D to D2, the temperature increase was represented by increasing R S and by a factor 2.8, R A by a factor 8.1, and P by a factor 2. From D to D3, R S and R A increased by larger factors. Both changes led to the formation of wider and taller islands, while island coalescence was visibly delayed. The same was observed in a change from D to D1. Figure 5a-c shows the evolution of the roughness W in sets C, D, and E, respectively. In each plot, the data for three possible conditions of higher temperatures are shown. The corresponding substrate covering thicknesses θ c are indicated. When deposits with the same thickness were compared, W increased with deposition temperature in the regimes of island growth, island coalescence, and the initial stages of continuous film formation.  However, this scenario may change in larger thicknesses. This is the case of the change from C to C2 in Figure 5a, in which the roughness at the highest temperature (C2) was smaller than that of the lowest temperature (C) for θ 400. Figure 6a,b shows crosssections of 640-layer films grown with parameters C and C2, respectively. In the former, W increased approximately as θ 1/2 , which is characteristic of random deposition [6] and explains the large fluctuations at short scales. In the latter, the small ES barrier facilitated the smoothening when the continuous film was formed, so only the long length scale fluctuation remained. This feature was typical of smallν, which indicates relatively large fluxes; see Table 1.   At the lowest temperatures of each panel (C, D, and E), s M was small (always < 10) and slowly increased as the deposit grew. At the highest temperatures (C3, D2, or E1), s M is displaced to much larger values. This is related to the formation of wider islands; for instance, comparing the images of deposits grown with sets D and D2 (Figure 4a,b). In all other sets studied here, when deposits of the same thickness were compared, s M also increased with the temperature in the regimes of island growth and coalescence.

Discussion
Before analyzing the results obtained here, it is useful to recall previous results of homoepitaxial deposition models, in which no difference in the adatom diffusivity is expected at different heights. If the ES barrier is negligible, simulations show that the film roughness decreases with temperature [9][10][11]. This occurs because larger diffusion lengths on terraces (related to R A ) and large probabilities of detachment from terrace borders (related to ) allow the adatoms to search for low energy positions at longer distances. If the ES barrier is present, the same trend may be observed for low coverages, in which the growth is nearly layer-by-layer, but changes appear in larger thicknesses, including possible re-entrant behavior [15,[17][18][19][39][40][41].
However, in the initial stages of heteroepitaxy, our simulations showed an apparently universal trend of the roughness and of the autocorrelation function, with or without ES barriers. When deposits with the same thickness were compared, in all cases, the roughness increased with the deposition temperature until formation of continuous films. The coverage range with these features then comprises regimes of island growth and coalescence (which are not quantitatively distinguished in our simulation data). Moreover, the position s M of the first minimum of the autocorrelation function increases with the temperature. These features are related to the formation of taller and wider islands. The roughness increase agrees with experimental observations in the initial growth stages of some materials [20][21][22]. The predicted behavior of s M may be checked in future studies by analyzing microscopy data.
Simple arguments explain the simultaneous increase of W and s M with the temperature. In most of our simulations, as the temperature increases, the diffusion-to-deposition ratio of adatoms on the substrate (R S ) weakly changes, but the detachment probability has larger changes. This leads to the increase of adatom mobility and, consequently, to the increase of the island width, as shown in studies of submonolayer growth [40,42,43]. Concomitantly, the larger mobility allows a larger number of adatoms to reach the island borders and move to their tops. For this reason, the islands also get taller as the temperature increases. This positive correlation between roughness and temperature is still observed after the formation of the continuous film. This suggests that the memory on the initial island configuration remains in that stage.
Our results also show that this trend may change in large thicknesses, after the continuous film is already formed. In one case, we showed that the roughness decreased with temperature in these conditions. It anticipates the possibility of re-entrant behavior, in which the variation of the roughness with the temperature is non-monotonic in a certain thickness range. Figure 8 illustrates this re-entrance when a temperature T(C ) higher than those of sets C and C1 is considered: between thicknesses θ ≈ 227 and θ ≈ 630, raising the temperature from T(C) = 350K to T(C1) = 436K leads to a slight decrease of roughness, but a subsequent raising to T C = 460K leads to a roughness increase. Observe that the physical parameters for deposition are the same in C, C1, and C .
The re-entrance of the roughness is a consequence of the different growth regimes of films that have the same thickness but that are deposited at different temperatures. At the lowest temperatures, they are in regimes of continuous film growth, in which the initial island configuration has weak influence. However, at the highest temperature, the islands were larger, taller, and took longer times to coalesce, so the larger roughness is an effect of that morphology.
Finally, note that our model considers ES barriers both for upwards and downwards motion, but this is not the reason for the observed behavior. For instance, if only barriers for downward motion were considered, the net adatom flux from the substrate to the top of the islands would increase. This would form taller islands, so the temperature effects described above would be enhanced. For these reasons, we expect that the features described here can be observed in several deposition processes on weakly interacting substrates, as long as the main assumptions of the model (e.g., the short range of the interactions) remain valid. At this point, it is interesting to recall that recent simulations of the initial stages of compound film growth also show an increase of the nanoparticle size and of the surface roughness with the temperature, which we understand are consequences of the same microscopic growth mechanisms [44].

Conclusions
The effect of the temperature in the initial stages of heteroepitaxial deposition on weakly interacting substrates was studied with kMC simulations. The deposition model accounts for the enhanced adatom diffusion on the substrate in comparison with the diffusion on the deposit [30]. We investigated the same quantities usually measured in experimental studies and in kinetic roughening theories: the surface roughness, which is defined as the root mean square fluctuation of the local heights, and the auto-correlation function of those heights [6,36].
When deposits with the same thickness are compared, we observe that an increase in the deposition temperature leads to larger roughness in the regimes of island growth, coalescence, and initial formation of a continuous film. Concomitantly, the position of the minimum of the autocorrelation function is displaced to larger sizes. These apparently universal results are consequences of the formation of wider and taller islands in those regimes and were observed with or without barriers for adatom diffusion at step edges.
Immediately after island coalescence, the same temperature effect is observed. However, after the continuous film is fully formed and the memory on the initial island configuration is lost, the opposite effect may be observed, i.e., a decrease of the roughness with increasing temperature. This was shown to lead to re-entrant behavior in some ranges of thickness and temperature.
Note that each sample simulated here is grown at a constant temperature, which is the most usual condition in experiments. Significantly different results are observed if the temperature changes during the growth of a film, possibly including apparent anomalous roughening [45]; however, this case is not considered here.
The increase of roughness with temperature was recently observed in the deposition of inorganic and organic materials [20][21][22]. Our simulations provide an interpretation of these experimental results, which are consequences of the formation of larger and taller islands in the initial growth stages. Our model also shows that measuring the autocorrelation function from microscopy images of growing films may be interesting to reinforce such interpretation.
Our results also broaden the spectrum of applications of the proposed model, which formerly described the initial stages of deposition of some materials [30]. Further extensions of the model, possibly considering the crystalline structures of specific materials, may provide quantitative descriptions of the initial growth stages, for instance with estimates of island sizes and coalescence times.