Coalescence of ultrathin films by atomic layer deposition or chemical vapor deposition: Models of the minimum thickness based on nucleation and growth rates

Fabrication of smooth, pinhole-free thin films is a challenge for manufacturing components of modern electronic and nanoscale devices. Films with these characteristics are desirable as diffusion barriers in metallic interconnects,¹ dielectrics in semiconductor devices,² protective coatings in electrochemical systems,³ and functional surfaces such as plasmonic metamaterials.^4,5 The growth process on a given substrate surface should afford complete substrate coverage (coalescence) within specified limits of film thickness and surface roughness. For many applications, such as diffusion barriers, the film must be rigorously free of bare areas (pinholes) because these defects lead to rapid device failure. If the application requires in-plane conductivity, the roughness will reduce the conductivity compared with a planar deposit containing the same quantity of material. Both outcomes imply that the additional material must be deposited to achieve the desired performance. However, this additional thickness may be incompatible with constraints of the device design, e.g., a diffusion barrier should not consume too much of the width of a deep trench or via. To design a successful growth process, it is crucial to understand how the nucleation and growth kinetics determine the film thickness and roughness at coalescence.

Low-temperature film growth by atomic layer deposition (ALD) or chemical vapor deposition (CVD) leads to an extremely conformal morphology because the growth kinetics are limited by the surface reaction rate. Microscopically, the growth surface reaches a near-saturation of precursor species, either by sequential dosing of coreactants in ALD (Refs. 6 and 7) or by a high ratio of a precursor flux to the surface reaction rate in CVD.^8–11 Conformal growth affords thick films that are smooth compared with a less-conformal process in which the growth rate varies with the flux or angular distribution of arriving precursor. However, even for a very conformal process, the initial stages of growth on a relatively unreactive substrate can cause the film to be rough at low thicknesses. This occurs when the deposit consists of islands of film that grow until coalescence. The onset of deposition is usually described by its effect in time as a “nucleation delay” or a “growth incubation period.”^7,12–17 While it is important to minimize growth time and to reduce process variability, it is the effect of nucleation on morphology that ultimately limits the performance of ultrathin films.

An important consideration is whether the nucleation behavior of the system can be controlled experimentally or whether it is predetermined by the precursor interaction with the given substrate. For example, the onset of HfO₂ deposition is governed by the surface concentration of hydroxyl groups.¹⁸ Plasma pretreatment of substrate surfaces improves nucleation density (which reduces the incubation time before growth can be detected) by removing contaminants and creating dangling bonds on the substrate.^19–21 CVD experiments demonstrate that the density of nuclei can also be enhanced by predosing the substrate with a molecular nucleation promoter.^22–27 When used in conjunction with growth rate inhibitors,¹¹ this can, for example, afford a 3.3 nm thick Co film with root-mean-square (RMS) surface roughness of only 0.35 nm.²² Thus, the utility of the present model is to indicate the target nucleation and growth rates that are needed to achieve a desired ultrathin film.

Numerous microscopic mechanisms may play a significant role in the nucleation and coalescence process. A challenge in modeling, both numerical and analytical, is to reduce the number of unknowns to the minimum set that is supported by the experimental data or by the known physical properties of the system. We show that for ALD and CVD, a useful baseline model can be constructed that considers only the nucleation rate and the isotropic growth rate of islands. It is appropriate, however, to acknowledge the larger set of mechanisms because they may need to be considered in cases where the baseline model proves to be insufficient. These include the distribution of reactive sites on the substrate; the spatial and angular distribution of arriving growth species; the adsorption, diffusion, and attachment kinetics of adspecies on both the bare substrate and the growing film; the surface diffusion of film atoms; and the surface energies that drive wetting behavior and island coarsening. When the substrate or the deposit is crystalline, then all the surface mechanisms will, in general, be a function of the crystallographic planes (and defects) involved.

Relevant models in the literature make varying assumptions. Geometric models use a surface unit cell to define where nuclei form, which in turn defines the nucleation density, and then follow the evolution of roughness, thickness, and coverage as a function of growth cycles (ALD) or time (CVD).^12,28 Monte Carlo simulations describe the evolution of CVD growth surfaces within a specified geometry of crystalline nuclei.^29,30 The effects of film growth on rough substrates have been addressed by assigning a distribution of spherical nuclei over a surface.³¹ A recent analytical model uses the Avrami assumption, which handles the issue of overlap between growing nuclei, to calculate the asymptotic approach of the surface coverage toward unity as a function of the initial density of nuclei or the ongoing nucleation rate.³² 

Here, we model the initial stages of film growth using two complementary approaches: (1) in a Monte Carlo model, nucleation is spatially random on bare areas of the substrate and can be instantaneous or random in time; (2) in an analytic model, Poisson statistics determine the probability of finding bare areas before coalescence. Both approaches assume hemispherical nuclei that increase in radius at a constant rate, i.e., film growth is conformal down to nanoscale dimensions. These assumptions are consistent with AFM data for the nucleation and growth of many amorphous films by ALD or low-temperature CVD.^8,22,33,34

The importance of nucleation is illustrated in Fig. 1, which was generated by the Monte Carlo simulation described below. In this example, nuclei are formed only at the start of the growth process and are located at random positions on the substrate surface [Fig. 1(a)]. As deposition proceeds, the nuclei grow into islands that progressively merge to form a continuous film [Fig. 1(b)]. We have halted the simulation just prior to complete coalescence: One tiny area of the substrate remains uncoated, which corresponds to a position on the substrate where, due to the statistics of the random location, no nuclei exist within the region indicated by the dashed white circle in Fig. 1(a). Device technologies require that there be no holes; hence, the minimum film thickness must be incrementally larger than in this simulation, such that the film fully covers the substrate. As we show, if the areal density of nuclei is low, then both the minimum thickness and the surface roughness will be large at the point when pinholes are eliminated; conversely, when the areal density of the nuclei is large, the film is smooth and pinhole-free at a small thickness.

The model employs a minimal number of physical parameters with the aim of revealing the role of the nucleation behavior without considering the additional mechanisms cited above. See the supplementary material⁴⁰ for a copy of the Python code. It has three functions to account for film growth on a planar substrate: one generates nuclei, one grows existing nuclei to produce a height distribution, and one determines the fractional surface coverage. Each of these functions occurs during a discrete time step. Growth proceeds from the nuclei as islands, i.e., the Volmer–Weber mode.³⁵ The fundamental assumptions are as follows:

1. Nuclei are hemispherical. This corresponds to an amorphous morphology and the absence of wetting behavior on the substrate surface.

2. Growth of nuclei (the formation of islands) occurs as a constant increase in the radius in all directions. This corresponds to perfectly conformal attachment of adspecies and the absence of favored growth directions.

3. Nuclei have an initial radius of 0 and begin to grow on the time step t after they are created, e.g., a nucleus created at t = 1 has no measurable radius until t = 2.

4. Nuclei are randomly placed on the surface of the substrate. The nucleation probability is the same on every bare site and zero on sites covered by the film.

5. When two islands intersect by growth, the height at every point of intersection is the larger value of the two islands had they each grown without intersection.

6. Islands do not migrate on the surface, and they do not coarsen upon contact.

The Monte Carlo implementation is as follows. A random number generator with an input nucleation probability determines whether nucleation occurs on each site within a 1000 × 1000 matrix with periodic boundary conditions. Nucleation can have two distinct probabilities: during the first step, which we call initial nucleation, and continuously in parallel with film growth, which we call ongoing nucleation. Initial nucleation occurs at t = 0, and ongoing nucleation occurs at or after t = 1. The matrix stores the film height z over the area of the substrate, where the indices of the matrix elements represent the in-plane x and y coordinates. The coordinates are multiplied by a scaling constant to afford physical dimensions. In our scaling, each surface site is a square 0.25 nm on a side; hence, the simulation substrate is a square 250 nm on a side. Once nuclei are created, the growth function updates the height on all coordinates that fall within the radius of the island. Island overlap is accounted for by choosing the highest value of all reference files that cover each point on the matrix. This treatment of overlap avoids double counting the growth rate and eliminates height discontinuities.

Several approaches are used to reduce computational time. The height of every point on a hemisphere of a given radius is calculated in advance and used (as a table lookup) to update the height profile of each island. In the case of ongoing nucleation, new nuclei are created at every time interval (step), but the net surface heights are only recalculated every four steps; this does not affect the outcome because of the treatment of overlap. After the surface coverage exceeds 95%, the heights are updated, and the surface coverage is calculated during each step. The model stops when 100% coverage is reached, but it can easily be changed to continue after coalescence. Data can be exported after a given number of steps or after a specified surface coverage is reached; however, the discrete growth steps used in the model cause the coverage to vary from the target value by a small amount. The code outputs a file containing the height data every time the height matrix is updated and a file containing the surface coverage values for all the output height matrices. The RMS roughness of the surface heights is calculated using a separate method.

Poisson point statistics can be used to derive an analytical relationship between nucleation density or ongoing nucleation rate, film growth rate, and the maximum film thickness at coalescence. This approach predicts, on a mathematical basis, the same power law relationships that are found in the Monte Carlo simulations (Sec. IV).

Poisson statistics describe a search space in which a process occurs infrequently with an average probability. Here, the space is an area of the substrate and the process is the formation of nuclei, expressed as N (area⁻¹) or N′ (area⁻¹ time⁻¹). The probability P that the process occurs exactly x times is

where μ is the average expected number in the search space. Variables used in this model and in the previous model are summarized in Table I.

Consider a film just prior to complete coalescence, such that only one spot of the substrate surface remains uncovered because no nucleus formed close to it. For a growth rate r and growth time t, the nearest nuclei must be separated from the spot by a radius (R) incrementally larger than R = r ⋅ t [Fig. 1(a)]. That radius defines a circular search space of area (πR²) in which no nucleation took place. Note that the angular position, the number of nuclei, and the detailed shape of the bare spot do not affect the result.

The Poisson probability of no nucleation, x = 0, is

For nuclei that form at the start of growth, the average number

According to Eq. (3), the probability of finding a bare spot approaches zero asymptotically with increasing R, whereas in a real film, complete coalescence (zero bare spots) is reached above a certain thickness. This apparent discrepancy occurs for three reasons. First, the nucleation may not be random but associated with chemically active sites that are never absent over unusually large areas. Second, in a real growth process, the bare spots may be filled in by other mechanisms, such as short-range diffusion of adspecies and preferential attachment in concave regions, which are not included in the present models. Third, the Poisson formula assumes an infinite substrate area, for which there is always a nonzero probability of finding a spot where all surrounding nuclei are unusually far away. A lower bound is the probability P that corresponds to one bare area divided by the finite area of the substrate under consideration; e.g., for one bare area over the 1000 × 1000 Monte Carlo matrix, P = 10⁻⁶.

Equation (3) shows that, for a chosen P, the argument of the exponential is constant, i.e., the maximum height R at coalescence varies as the inverse square root of the areal density of nuclei N. This is the expected result based on geometry. A useful form emerges when we rewrite P as

Taking the natural logarithm and rearranging,

For example, α = 6 (a probability P = 10⁻⁶ of a bare spot) requires a density of nuclei N = 14/πR². To achieve coalescence at a maximum film thickness of R = 3 nm, N must be 0.5 nuclei nm⁻², which is a high but not unphysical value. This vividly illustrates the challenge of reaching full coalescence in a very thin film if growth proceeds by random nucleation and spherical growth in the absence of other smoothing mechanisms.

For a fixed density of nuclei, the growth rate r determines how rapidly coalescence is reached but does not affect the morphology under the given assumptions. This implies that reducing the growth rate, e.g., via coflow of an inhibitor, is not needed. However, many available CVD precursors have relatively high reactivity that leads to subconformal growth; in this case, the peaks of the islands will grow faster than the diameters, which increases roughness and delays coalescence. The addition of an inhibitor affords highly conformal growth that promotes smoothness. For this reason, our CVD process for ultrasmooth cobalt films employs both a nucleation promotor to increase N and an inhibitor to decrease r.³⁶ 

If nucleation occurs at a rate N′ (area⁻¹ time⁻¹), then the average expected number of nuclei within the search space, μ, is calculated from an integral over the growth time. It is convenient to define a dummy variable, τ, as the growth time for any nucleus, i.e., a nucleus that forms at the start of the process can grow for the full duration t, whereas a nucleus that forms later can only grow for a shorter time until the process ends. Therefore,

It is useful to rewrite the growth time t in terms of R, the maximum radius (height) of islands at coalescence, which are those for which nucleation occurred at the start of the process. Then,

The equivalent to Eq. (5) is

The maximum height R at coalescence, therefore, varies as (N′/r)^−1/3. Small R requires large values of N′/r, the ratio of the nucleation rate to the growth rate, but the cube root makes the dependence slow. The reason that ongoing nucleation has only a modest effect on the result is that the bare area available for nucleation decreases continuously during the process. On a log-log plot,

Note that the chosen value of α offsets the curve vertically, but it does not change the slope of −1/3.

Figure 2(a) shows the result of the Monte Carlo simulations that are terminated at a surface coverage of 99%, together with analytical results for the same coverage, α = 2, in the case of a fixed density of initial nuclei. The maximum film thickness decreases as the square root of the density of initial nuclei (a slope of −0.5 on the log-log plot), as predicted by Eq. (5). This relationship between the nucleation density and the time (or thickness) to reach a fixed surface coverage is identical using the Poisson formalism or using the Avrami assumption.³² This is the case because the same mathematical form emerges [compare Eq. (3) in this article and Eq. (8) in Ref. 32] despite the very different assumptions used to construct these models. Results for the maximum thickness and coverage in the Monte Carlo simulation agree very well with the analytical results; this is expected because these methods utilize the same assumptions for nucleation density and growth rate on the substrate surface. The Monte Carlo simulation shows that the RMS roughness scales linearly with the maximum film thickness, as expected from geometric scaling; the roughness at 99% coverage is ∼1/7 of the maximum thickness.

For the case of ongoing nucleation, the important parameter is the ratio (N′/r). Figure 2(b) shows the results for the Monte Carlo simulation and the analytical model at the same terminations shown in Fig. 2(a) (∼99% coverage and α = 2) for a growth rate of 0.25 nm min⁻¹. The maximum film thickness decreases as the inverse cube root of the ongoing nucleation rate (a slope of -1/3), as predicted in Eq. (8). Similar to the case of initial nucleation [Fig. 2(a)], the RMS roughness scales linearly with the maximum film thickness. In both cases, small fluctuations in the Monte Carlo predictions occur due to the finite size of the model and because the simulation halts on a time step where the coverage is slightly different from the target value of 99%.

Figure 3 shows the Monte Carlo output at a coverage of ∼98.5% for two different ongoing nucleation rates; we display this result, prior to 100% coverage, to emphasize the roughness and bare areas of the substrate. The simulation parameters for these height maps, including the duration of growth and maximum thickness, are listed in Table II. The differences are striking: the simulation with the higher nucleation probability of 1 × 10⁻⁴ site⁻¹ step⁻¹ [Fig. 3(a)] has a RMS roughness of 1.6 nm, which is much smoother than the roughness of 3.6 nm obtained from a nucleation probability of 1 × 10⁻⁵ [Fig. 3(b)]. The uncoated areas are higher in number but smaller in area for the higher nucleation probability. In addition, many additional simulation steps, i.e., a higher maximum thickness, are required for the film with the lower nucleation rate to reach 100% coverage. According to the analytical model, a decrease in the ongoing nucleation rate by a factor of 10 will increase the maximum film thickness by a factor of 10^1/3 (=2.15) for any value of α. This factor fits the trend for maximum thickness in Table II (where the ratio of growth steps for t₁₀₀ is 106/50 = 2.12), which is also the trend for roughness and growth time as explained above. The simulated films are similar in appearance to SEM images and AFM height maps of precoalescence films in other nucleation studies.^8,22

At the onset of nucleation, the steep sides on growing islands cause a rapid increase in roughness; smoothing takes over at higher coverages as islands merge to form a continuous film. This is a well-known phenomenon in both ALD and CVD that motivated earlier modeling efforts.^12,28 Figure 4 plots these metrics from the height maps of simulated films as a function of the number of simulation steps for an ongoing nucleation probability of 10⁻⁴ site⁻¹ step⁻¹ and a growth rate of 0.25 nm step⁻¹. Growth per cycle (GPC) is calculated as the difference in average thickness divided by the number of growth steps between measurements; GPC peaks at 0.33 nm step⁻¹. GPC exceeds the local growth rate when the surface area of the islanded microstructure, and thus the attachment rate of adspecies, is larger than the projected area of the substrate. Determining growth per cycle through Monte Carlo simulation affords a simple route to the complicated analytical problem of finding the increase in volume per step from a random arrangement of overlapping nuclei.

Roughness increases as the height of nuclei above the bare substrate increases, but it begins to decrease significantly once a large portion (or all) of the bare substrate is coated. The peak roughness at t = 26 steps occurs at 85% coverage. Up to and after 100% coverage is reached, the roughness continuously decreases because the curvature of hemispherical nuclei decreases (the radii increase). The effect of smaller nuclei will also be lost as they are overgrown by larger nuclei. The power spectral density of surface heights would show a reduction in short-range roughness. What the present simulation will not reproduce is a simultaneous increase in long-range roughness that can occur experimentally when the growth rate depends on the flux and the incident flux is slightly shadowed by surface topography.³⁷ 

We recently studied the nucleation behavior of HfB₂ films grown by CVD on SiO₂ and Al₂O₃ surfaces.²² Using the single-source precursor Hf(BH₄)₄ at 220 °C, these surfaces afford selective deposition: nucleation is very slow on SiO₂ but rapid on Al₂O₃. The resulting films after a long growth time are discrete islands on SiO₂ and smooth, continuous films on Al₂O_3. Film growth was interrupted at different points during the nucleation process in order to measure the density of nuclei and film roughness by AFM; the areal density of deposited Hf atoms was measured by Rutherford backscattering spectrometry and converted to an equivalent bulk thickness, which can be compared directly to the average height over the matrix in the Monte Carlo simulation.³⁸ 

For growth on Al₂O₃, the simulated thickness and roughness are in reasonable agreement with the experimental data [Fig. 5(a)]. However, the growth rate of islands, as evaluated by the tallest ones in the distribution, is twice as fast during the initial period, 0–4 min, than in steady-state (data not shown). The model agreement can be improved by doubling the growth rate over the first 4 min. The significant point is that another mechanism must be involved, for example, the lateral diffusion of an adsorbed precursor from the bare surface to the growing islands, which allows small islands to grow faster. At growth times longer than 4 min, the maximum height of nuclei on this substrate and on SiO₂ increases linearly with time, which suggests that this effect does not govern the majority of the nucleation process.

On SiO₂, the simulated thickness and roughness considerably exceed the experimental results [Fig. 5(b)]. The discrepancy indicates that HfB₂ growth on SiO₂ must involve other mechanisms. AFM data acquired at different times (not shown) indicate that the nucleation probability increases during the initial period, 0–25 min, before stabilizing at ∼30 nuclei μm⁻² min⁻¹. This behavior may imply that the presence of growing islands promotes the nucleation of additional islands; this trend is the opposite of what would be expected with a denuded zone in which island growth lowers the density of the mobile precursor on the surface. A working hypothesis is that reaction by-products, such as BH_x, are released from the growing islands and promote nucleation on the bare SiO₂ surface. We previously advanced this hypothesis to explain the observed nonrandom pair correlation function for HfB₂ islands on the SiO₂ surface.³⁹ 

We present two modeling approaches that capture the effect of the film nucleation rate and the film growth rate on the evolution of coverage, total thickness, and roughness during the onset of growth. (1) A Monte Carlo simulation generates a height profile of simulated films as a function of the number of simulation cycles. The height maps are useful for showing the development of film height, roughness, and the distribution of bare areas on the substrate, and they simplify the calculation of these parameters over the course of growth. (2) Poisson point statistics are used to evaluate the last site to be coated on the substrate, which is the site with the largest radius away from any nuclei. This analytical model affords the same power law relationship between maximum thickness and film nucleation as the Monte Carlo simulation.

Growth inhibitors and nucleation promoters are essential in designing film processes that afford a smooth, pinhole-free morphology at low thickness. When an experimental system of interest does not intrinsically afford a high density of nuclei or a high nucleation rate relative to the film growth rate, then kinetic modifications are required. The present models can be used to determine the rates that are needed. Modifications consist of increasing the nucleation density, e.g., by predosing the surface with a promoter or decreasing the growth rate by adding a coflow of inhibitor. In combination, these techniques can afford ultrasmooth and ultrathin films.

D.K.L. and K.L.C. contributed equally to this work.

The authors acknowledge project funding from the National Science foundation (NSF) through No. CMMI 1825938 (J.R.A.) with a Research Experience for Undergraduates (REU) supplement (D.K.L.), Grant No. CER 1914769 (K.L.C.), and the REU program of the Materials Research Science and Engineering Center at the University of Illinois at Urbana-Champaign, No. DMR 1720633.

The data that support the findings of this study are available within the article and its supplementary material.⁴⁰ 

Diana K. LaFollette is currently a first year PhD student in the Materials Science and Engineering Department at the Georgia Institute of Technology. The research described here occurred while she was participating in the Materials Research Science and Engineering Center REU program at University of Illinois at Urbana-Champaign. While completing her Bachelor's degree in Chemical Engineering at the University of Southern California in 2021, she was awarded the National Science Foundation Graduate Research Fellowship. She is now studying the structure-property relationships of lead-halide perovskites for solar cell applications using advanced characterization methods with Dr. Correa-Baena.

Kinsey L. Canova is in her fifth year pursuing a Ph.D. at the University of Illinois at Urbana-Champaign in Urbana, Illinois, USA. She received her Bachelor's degree in Materials Science and Engineering at the Georgia Institute of Technology in 2017. Her Ph.D. research focuses on the kinetic controls available in low temperature, surface-reaction limited chemical vapor deposition of thin films, and she is particularly interested in how these effects can be used to engineer growth processes that enhance smoothness, filling in recessed features, and control over relative deposition rates in alloy growth. The applications for these process designs include integrated electronic device manufacturing and protective coatings.