Thin film deposition with time-varying temperature

The model studied in this work was introduced in [19]. It is a solid-on-solid type model (no overhang in the film surface) and the deposit has a simple cubic lattice structure. The substrate is flat, at the xy plane, with linear size L (L × L columns), with the lattice parameter being the length unit. The column height is the z coordinate of the topmost adatom. There is an external flux of F atoms per site per unit time, measured in number of monolayers (ML) per second. An incident atom is adsorbed upon landing above a previously deposited atom or a substrate site. All adsorbed atoms with no lateral and no upper neighbor diffuse with coefficient D (number of steps per unit time). If an adatom has a lateral or an upper neighbor, then it is permanently aggregated at that position. Figure 1 illustrates the possible steps of some mobile atoms at the film surface.

Zoom In Zoom Out Reset image size

The diffusion coefficient D depends on temperature as

$$\begin{eqnarray} &&D={\nu }_{0}\exp \left (-E/{k}_{\mathrm{B}}T\right ), \end{eqnarray} \tag{ 1 }$$

where ν₀ is a jump frequency and E is the activation energy in a flat surface, which typically amounts to tenths of eV. The parameter

$$\begin{eqnarray} &&R\equiv D/F \end{eqnarray} \tag{ 2 }$$

quantitatively represents the interplay between temperature and pressure and determines the film properties. The other relevant parameter is the film thickness t, proportional to the deposition time.

In the submonolayer growth regime, this model corresponds to irreversible island growth with critical nucleus of size i = 1, which has already been applied to several systems [4]. The irreversible lateral aggregation is a reasonable assumption for low temperature or large atomic flux, although it does not respect detailed balance conditions. This model also has the advantage of reducing simulation times compared with models where all atoms are mobile. The aggregation condition is similar to the Das Sarma and Tamborenea (DT) model of molecular beam epitaxy [20], but the main difference is that the DT model and its extensions restrict the adatom diffusion to finite distances. Moreover, the model does not consider Ehrlich–Schwoebel (ES) barriers for downhill movement at terrace edges [21] and does not allow uphill movement because adatoms with lateral neighbors are immobile. Indeed, our aim is to consider a model with a minimal set of relevant parameters, in order to search for features that are intrinsically related to the time-varying temperature.

Simulation results presented here were obtained in lattices with L = 512. Some results were also obtained in L = 256 and 1024 in order to confirm that finite-size effects are negligible in the simulated time range. Fixed deposition parameters are ν₀/F = 10¹⁴ and the activation energy E = 0.6 eV. These values are reasonable compared to experimental data for several metals and semiconductors, although the direct applicability is limited due to the differences in lattice structure and aggregation assumptions [22, 4]. Choosing a fixed value of E and changing T does not restrict our conclusions because the physically important parameter is R, which is a function of E/T. In all cases, the temperature is uniform across the film and may vary in time. Some results for fixed temperature are also obtained for comparison.

We will consider temperature values in which R ranges between 10¹ and 10⁵. These values are small if compared to works on submonolayer growth, which frequently consider R up to 10⁹ [4]. However, they are suitable for low-temperature deposition, which is a necessary condition for the applicability of a model with irreversible attachment to steps (in other words, even with a small binding energy to a lateral neighbor, the relative probability of an attached atom to move will be small). For R = 10⁵, results of [19] show that adatoms execute an average of 28 random steps before aggregation, thus diffusion lengths in terraces are much smaller than the lateral size of the films and finite-size effects are negligible.

During the growth of a sample, the adsorption of each new atom takes place in a time interval 1/(FL²). After this process, R/L² steps of randomly chosen free atoms are performed. Since R/L² is usually not integer, and it may be small for large L, we keep its fractional part for determining the number of steps after the next adsorption events. If the temperature changes during the film growth, then the value of R is updated after the deposition of each complete layer (L² atoms), which corresponds to a unit time interval.

The global roughness (or interface width) is defined as

$$\begin{eqnarray} &&W\left (L,t\right )\equiv {\left \lt \overline{\left (h-\overline{h}\right )}^{2}\right \gt }^{1/2}, \end{eqnarray} \tag{ 3 }$$

where the overbars represent spatial averages and the angular brackets represent configurational averages over various samples. For short deposition times (small thicknesses), effects of the finite substrate size L are negligible and W scales as

$$\begin{eqnarray} &&W\sim {t}^{\beta }, \end{eqnarray} \tag{ 4 }$$

where β is called the growth exponent.

For comparison with experimental data, the local roughness w(r,t) is more useful. It is also defined as a root mean square height fluctuation, but the spatial average is limited to a box of size r, with r < L, and the configurational average is performed by gliding this box over the film surface and over various samples. In the regime of thin film growth (negligible effects of substrate size), it follows the Family–Vicsek (FV) dynamic scaling relation [5]

$$\begin{eqnarray} &&w\left (r,t\right )\approx A{r}^{{\alpha }_{l}}g\left (\frac{B t}{{r}^{z}}\right ). \end{eqnarray} \tag{ 5 }$$

Here, α_l is the local roughness exponent (sometimes called Hurst exponent), z is the dynamic exponent, B is a constant, and g is a scaling function such that g(x) ∼ 1 for x ≫ 1 (small box size) and g(x) ∼ x^−α_l/z for x ≪ 1 (large box size, where the local roughness coincides with the global one).

In systems with normal scaling, the amplitude A in equation (5) is a constant (time-independent), thus β = α_l/z and the global roughness exponent α equals the local one. A scaling analysis of w(r,t) of several growth models in one- and two-dimensional substrates is presented in [23]. In systems with AS, the amplitude A scales as

$$\begin{eqnarray} &&A\sim {t}^{\kappa }, \end{eqnarray} \tag{ 6 }$$

with κ characterizing the degree of anomaly [8]. A FV relation can also be defined for the global roughness and involves the exponent α, but it is useful only if very long times/large thicknesses are attained.

Here we restrict the discussion to results on two-dimensional substrates (growth in 2 + 1 dimensions), which is the subject of the present work.

Several models with diffusion of all adatoms and energy barriers proportional to the number of neighbors were previously studied in fixed temperature. Some authors proposed that they had AS [28] or a logarithmic scaling of amplitude A in the FV relation [29]. Wilby et al [30] showed that the growth exponent β was close to the value of the universality class of the fourth order nonlinear growth equation of Villain, Lai and Das Sarma (VLDS) [31, 32]. Recent renormalization studies explained the long crossover to VLDS scaling [33, 34].

The model presented in section 2.1 has previously been studied numerically, showing a FV relation for the global roughness that includes the R-dependence as [19]

$$\begin{eqnarray} &&W=\frac{{L}^{\alpha }}{{R}^{x}}f\left ({R}^{y}t/{L}^{z}\right ), \end{eqnarray} \tag{ 7 }$$

with α ≈ 2/3, z ≈ 10/3, x ≈ 0.5, and y ≈ 1. The exponents α and z agree with those of the VLDS class [3, 35, 36]. The values of exponents x and y were explained by scaling arguments [19].

For R ≥ 10⁵, the surface roughness of this model is very small, even with very large thicknesses (10⁴ or 10⁵ ML) [19]. This is not a regime of dynamic scaling, which is an additional reason for our simulations to be performed with smaller values of R. Full diffusion models with values of R of the same order may show larger roughness, particularly if ES barriers are included [37], but with an extended set of parameters.

A recent work on numerical integration of the VLDS equation showed evidence of AS in 2 + 1 dimensions [38], which raises the question on a possible AS in our lattice model. This feature was not investigated in  [19].

We performed simulations of the model with 10 ≤ R ≤ 1000 for fixed R (fixed temperature). Figure 2 shows the local roughness as a function of box size for three values of that parameter. For R = 10, there is a split in the curves for small r, which suggests AS. However, that split tends to disappear as the thickness increases. For instance, from t = 4000 to 8000 ML, the change in the local roughness for r = 5 is only 7%. For R = 100 and 1000, the split of the curves is almost negligible for t > 10³. For those reasons, we understand that there is no asymptotic AS in our lattice model with fixed R, and an apparent anomaly appears only for small R and for small thicknesses due to some type of scaling correction.

Zoom In Zoom Out Reset image size

Figure 2 also shows that the local roughness scales with α_l ≈ 1/3 for small r (equation (5)). This exponent is significantly below the VLDS value α_l = α ≈ 2/3 [36]. This discrepancy it not a failure of theoretical predictions, but a consequence of scaling corrections also observed in other growth models [23].

First we consider linearly decreasing temperature during the time interval for deposition of t_max = 10⁴ ML, numbered cases (I) and (II), whose conditions are specified in table 1. This value of t_max is reasonable for mesoscopic films, giving thicknesses of a few micrometers for most materials. Figure 3(a) shows the time evolution of the parameter R and of the temperature in those cases.

Zoom In Zoom Out Reset image size

We also consider cases of exponential convergence of the temperature to the final value, with a characteristic time/thickness t_c = 2000 ML, so that

$$\begin{eqnarray} &&T={T}_{\mathrm{F}}+\left ({T}_{\mathrm{I}}-{T}_{\mathrm{F}}\right )\exp \left (-t/{t}_{\mathrm{c}}\right ). \end{eqnarray} \tag{ 10 }$$

The conditions of cases (III) and (IV) considered here are presented in table 1. Deposition of t_max = 10⁴ ML = 5t_c layers is considered, thus the final temperatures of the films are very close to T_F. These cases represent systems exchanging heat by conduction with a colder reservoir (e.g. an initially heated system separated from the surroundings by a conducting wall).

Figure 3(b) shows the corresponding time evolution of the parameter R and temperature. The main difference from the linear decays of cases (I) and (II) is that R rapidly decreases at short times (up to t ∼ t_c) and slowly decreases during most of the deposition time. When t ∼ t_max, the temperature is approximately constant (T_F).

In cases (I)–(IV), E = 0.6 eV in the room temperature range, thus E/k_BT ∼ 23. Maximal temperature changes are close to 10% of the initial ones, from 1 to 10⁴ ML, thus dlnT/dlnt ∼ 0.01. This gives s ∼ 0.2 (equation (9)), which is expected to give a significant change in the roughness scaling (section 3). Depending on the particular form in which the temperature varies, larger or smaller values of s may be obtained in different time ranges, as shown below.

The control of these temperature changes are realistic for thin film growth. For comparison, in film deposition with the substrate subject to a temperature gradient, much larger temperature differences are established [26, 27]. An example is FePt film growth with the temperature varying from 250 to 600 °C along a substrate distance smaller than 1 cm [26].

Figure 4(a) shows the evolution of the global roughness of the films in cases (I) and (II) and the theoretical predictions from equation (8) with the time-dependent R shown in figure 3(a). Those predictions are close to the simulation data when the roughness is larger than 1. The discrepancy for W ≲ 1 is expected because scaling relations are not expected to apply for small roughness. The local slope of the plots in figure 4(a) are effective exponents β_eff.

Zoom In Zoom Out Reset image size

For small thicknesses (t ≤ 10³ ML), the slope β_eff is close to the VLDS value β ≈ 0.2. This corresponds to one tenth of the total growth time and the temperature change is only 3 K, thus it is a growth with nearly constant temperature. For case (I), the roughness is typically smaller than unity, thus it is not a true scaling region where a reliable estimate of β can be measured. For larger thicknesses (t ≥ 10³ ML), W rapidly increases in cases (I) and (II). For 5 × 10³ ML ≤ t ≤ 10⁴ ML, the effective exponents are β_eff ≈ 0.70 and β_eff ≈ 0.82, respectively (see figure 4(a)).

The evolution of the global roughness for cases (III) and (IV) is shown in figure 4(b) and also compared with the extension of equation (8) to time-dependent R (figure 3(b)). Again, W increases slowly for small thicknesses, but shows large effective exponents β_eff for thicknesses above 10³ ML. In case (IV), where R decreases almost two orders of magnitude during the film growth, β_eff > 1/2 is found, similarly to cases (I) and (II). As t approaches t_max, the slopes of the plots in figure 4(b) tend to decrease due to the temperature saturation.

The continuously increasing β_eff, which attains large values at longer times, represents a nontrivial evolution of the roughness. This is solely a consequence of the time-decreasing temperature during deposition, which reduces diffusion lengths and facilitates roughening. The agreement with the theoretical approach of section 3 supports this interpretation.

Some models and experiments show large values of β_eff, but for different reasons. An example is Ag/Ag(100) growth, where mound steepening leads to a crossover from very smooth surfaces (up to 25 ML) to very rough ones (up to 1000 ML) at 300 K [39]. In that case, the presence of ES barriers [21] is responsible for the crossover. Some systems also have β > 1/2 as a consequence of preferential aggregation at surface hills instead of valleys (possibly with surface instability) or shadowing effects [40, 41]. However, in our model, none of those mechanisms are present.

Figure 5 shows the surface topography of films with thicknesses t = 500 and 8000 ML in case (IV) and with fixed temperature T = 270 K. The film grown in fixed temperature is always smooth, but a remarkable roughening is observed in the film grown with decreasing temperature. In the latter, there is no evidence of mound formation, columnar growth, or any other geometrical feature that frequently explains the large values of β_eff. Instead, the nontrivial evolution of the roughness is a consequence of the decreasing temperature.

Zoom In Zoom Out Reset image size

We also recall that the scaling in systems with competitive growth dynamics is very different from figures 4(a) and (b), since they usually show two scaling regions with β_eff ≤ 1/2. Moreover, large β is not necessarily related to large roughness. For instance, very rough deposits are produced by grain deposition models [42] and have β ∼ 0.1, which is much smaller than the value of the asymptotic Kardar–Parisi–Zhang class [43].

For the above reasons, if a nontrivial increase of the roughness is observed in an experimental work, one must consider the possibility that the temperature is changing during the deposition, since this seems to be the simplest mechanism responsible for that feature. Alternatively, the possibility of time-increasing adsorption rates may be considered, as in [11] (where β ≈ 0.5 was obtained), since they also lead to time-decreasing R.

The local roughness in cases (I) and (II) is shown in figure 6(a) as a function of the box size r, for several thicknesses in the range 500 ML < t < 8000 ML. There is a significant split of the curves for small r, in contrast to the small split shown in figure 2 for constant temperature (note that vertical and horizontal scales are the same in figures 2 and 6(a)). This is an AS feature, with the FV relation (equation (5)) having time-dependent amplitude A.

Zoom In Zoom Out Reset image size

Figure 6(b) shows the local roughness as a function of box size for cases (III) and (IV), respectively. These plots also show evidence of AS due to the split of the curves for small r. This trend persists for thicknesses t > 2t_c, where the temperature is changing very slowly.

The AS, corresponding to the continuous increase of local slopes, is another feature intrinsically related to the slow-down of surface diffusion. Indeed, reduced adatom mobility is known to facilitate the formation of steep surface features at small lengthscales.

For small and constant r, the time scaling of the amplitude A (equations (5) and (6)) characterizes the AS. The inset of figure 6(a) shows w(r = 10,t) (proportional to A) as a function of time in cases (I) and (II). It increases faster than the power-law of equation (6), particularly for the longest times. On the other hand, the inset of figure 6(b) shows the scaling of w(r = 10,t) for cases (III) and (IV), which approximately follow power laws with anomaly exponents κ = 0.28 ± 0.01 (III) and 0.47 ± 0.02 (IV). This reinforces the evidence of AS in plausible conditions for experimental work only due to the decreasing temperature.

It is important to stress that the AS found in cases (III) and (IV) is an apparent scaling restricted to the thickness range studied here. At long times, the temperature saturates at T_F, thus the deposit will attain the characteristic features of films grown with fixed temperature. Such films do not have AS, but have normal VLDS scaling (section 2.4). Furthermore, in cases (I) and (II), the temperature decrease is necessarily limited to a finite time range, thus there is no asymptotic scaling. These features contrast with the true AS of the lattice models in [18], but those models were not so closely related to real film growth.

Finally, we observe that the local roughness exponents measured in figures 6(a) and (b) are close to 1/3, similarly to the constant temperature case (section 2.4).