Conformality of atomic layer deposition in microchannels: impact of process parameters on the simulated thickness profile

Unparalleled conformality is driving ever new applications for atomic layer deposition (ALD), a thin film growth method based on repeated self-terminating gas–solid reactions. In this work, we re-implemented a diffusion–reaction model from the literature to simulate the propagation of film growth in wide microchannels and used that model to explore trends in both the thickness profile as a function of process parameters and different diffusion regimes. In the model, partial pressure of the ALD reactant was analytically approximated. Simulations were made as a function of kinetic and process parameters such as the temperature, (lumped) sticking coefficient, molar mass of the ALD reactant, reactant’s exposure time and pressure, total pressure, density of the grown material, and growth per cycle (GPC) of the ALD process. Increasing the molar mass and the GPC, for example, resulted in a decreasing penetration depth into the microchannel. The influence of the mass and size of the inert gas molecules on the thickness profile depended on the diffusion regime (free molecular flow vs. transition flow). The modelling was compared to a recent slope method to extract the sticking coefficient. The slope method gave systematically somewhat higher sticking coefficient values compared to the input sticking coefficient values; the potential reasons behind the observed differences are discussed.


Introduction
2][3] ALD enables one to make conformal coatings on almost any desired inorganic substrate including high aspect ratio (HAR) structures such as microelectronics and powder media.Yet, tuning of the process parameters is often required to guarantee conformal coatings in HAR structures.
Several types of feature-scale models have been used to simulate ALD growth in high aspect ratio (HAR) structures [e.g.Fig. 1(a)], as recently reviewed by Cremers et al. 4 Analogously to a recent article on chemical vapor deposition, 5 we classify these ALD models as ballistic line-of-sight (e.g.ref. 6-9), Monte Carlo (e.g.ref. 10-17) and diffusion-reaction models (e.g.ref. 18-23).§ While the heterogeneous gas-solid reactions responsible for ALD growth have been demonstrated in a great range of pressures from atmospheric to ultra-high vacuum, 4,17,24 ALD processes often operate in a low vacuum of roughly 10 2 Pa range. 3Consequently, most feature-scale models for ALD have been developed for low-pressure conditions where the mean-free-path of the molecules l is much higher than the limiting dimension of the feature h (Knudsen number Kn = l/h c 1). 9,25 Here, molecules collide with the feature walls and not with other molecules in the gas phase, and the mass transport regime is referred to with varied names e.g. as (free) molecular flow, Knudsen flow, or Knudsen diffusion. 4,19,26,27Diffusion-reaction models based on Fick's law of diffusion can flexibly be used in free molecular flow (Kn c 1) as well as in transition flow (Kn E 1) and even under continuum flow (Kn { 1) conditions, 25,27,28 as the effective diffusion coefficient D eff (m À2 s À1 ) can be calculated from the gas-phase diffusion coefficient D A (m À2 s À1 ) and the Knudsen diffusion coefficient D Kn (m À2 s À1 ). 17,18,202][13] Irrespective of the theoretical framework, all models reproduce the typical profiles of ideal ALD growth based on self-terminating (i.e., saturating and irreversible) reactions: 3,29 constant film thickness followed by an abrupt decrease to zero, as illustrated in Fig. 1(b).
For ALD process modelling and reactor design, a useful description of the reaction kinetics is essential.Typically, the reaction kinetics of ALD processes are described in a simplified manner assuming irreversible single-site Langmuir adsorption and an associated (lumped) sticking coefficient. 4,7,13,17,20,301][32][33][34] The most straightforward way of analysing the kinetics is by interpreting film termination profiles measured in dedicated test structures or in a cross-flow reactor, where a steep film termination profile is generally associated with high reactivity. 11,12,30,31,35,36With a low sticking coefficient, the process can be in a reaction-limited regime, where no clear termination profile can be identified.Whether film growth is in a reaction-limited or diffusionlimited regime, can be estimated by the Thiele modulus h T (or a 17,20,31 where a = h T 2 ) which gives the ratio of the reaction rate to the diffusion rate. 27In the diffusion-limited regime, the value of Thiele modulus is much higher than one. 27,370][41] Furthermore, a slope method has been developed by Arts et al. 31 to be used in conjunction with microscopic LHAR test channels, where the sticking coefficient c can be calculated from the slope [at a surface coverage y (À) of 1/2] of the Type 1 normalized thickness profile 40 [Fig.1(b)] through a simple square root relationship. 31Here, a Type 1 normalized thickness profile 40 refers to the normalized amount of growth (one for a saturated surface, or y = 1) as the vertical axis and the dimensionless distance (distance divided by channel height) as the horizontal axis.This slope method was derived empirically from the diffusion-reaction model of Yanguas-Gil and Elam 20 at free molecular flow conditions.
In this work, we have re-implemented the diffusion-reaction model by Ylilammi et al. 18 and used it to simulate the evolution of ALD growth at various scenarios of kinetic and process parameters and diffusion regimes.We first describe the assumptions and equations behind the Ylilammi et al. 18 diffusion-equation model.We then demonstrate how process parameters influence the ALD thickness profile by varying individual parameters at various diffusion regimes and channel filling levels.Finally, we compare the simulations of the Ylilammi et al. 18 model with the Arts et al. 31 slope method and discuss the likely reasons for the observed slight differences.

Description of the Ylilammi et al. 18 diffusion-reaction model
For full details behind the model derivation, please see the article by Ylilammi et al. 18 Here, core concepts are presented that are used in the current implementation.In some cases, somewhat expanded explanations are provided, to help the reader follow the model and connect it to other models on ALD.

Basic ALD process and geometry assumptions
The Ylilammi et al. 18 model was built to describe a typical ALD process, based on the use of inert gas for transporting the reactant from the source to the growth surface. 3he considered high-aspect-ratio (HAR) geometry is a wide microchannel similar to the one in Fig. 1(a), where the height H (m) of the microchannel is the limiting dimension.The width W (m) of the microchannel is orders of magnitude larger than H, and the length L (m) is considered infinite (the channel end effects are not considered).While the model has been constructed with lateral HAR (LHAR) structures in mind, 39,40 it is Typically, an ALD process has at least two reactants, often called Reactant A and Reactant B. 3 In this model, one of the two reactions is assumed to limit the extent of film growth in the microchannel; typically, this is assumed to be the reaction of Reactant A. The partial pressure of Reactant A is denoted as p A (Pa).Reactant A is brought to the microchannel entrance at a partial pressure p A0 . 18It is expected that the partial pressure of Reactant A behaves like a step function: during the reactant pulse, the pressure at the microchannel entrance is p A0 , and otherwise it is zero. 3,42An inert carrier gas, denoted here with ''I'' (instead of the notation ''B'' used in the Ylilammi et al. 18 article, to avoid confusion with Reactant B 3 ), is used to aid the transport of Reactant A from the source to the surface.The inert gas has the same partial pressure p I (Pa) inside and outside of the microchannel.Inside the microchannel, p A decreases, as Reactant A is consumed in the adsorption (i.e., ALD) process.The time t from the beginning (t = 0) until the end of the exposure of Reactant A (t = t 1 ) is considered; purge is excluded from this model.Thus, of the four typical steps in an ALD sequence, 3 the current simulation concerns Step 1 (or Step 3) only.

Mass transport by diffusion and partial pressure of Reactant A
In the chosen geometry, the partial pressure of Reactant A can be considered constant in the y-and z-direction. 18The onedimensional diffusion equation for the partial pressure of Reactant A p A (eqn (10) of Ylilammi et al. 18 ) is: The second term on the right side of this equation is called the adsorption loss term.R is the gas constant (8.314461J K À1 mol À1 ), T is the absolute temperature (K), and N 0 is Avogadro's constant (mol À1 ).The effective diffusion coefficient D eff considers both the gas-phase collisions and the channel wall collisions through the gas-phase diffusion coefficient D A (m 2 s À1 ), as well as the Knudsen diffusion coefficient D Kn (m 2 s À1 ), in the Bosanquet relation (eqn (6) of Ylilammi et al. 18 ): In eqn (1), g (m À2 s À1 ) is the net adsorption rate of molecules from the gas phase to the surface.h (m) is the hydraulic diameter of the microchannel (eqn (5) of Ylilammi et al. 18 ): The gas-phase diffusion coefficient depends on the average speed of the Reactant A molecules % v A (m s À1 ) and the collision rate of the Reactant A molecules in a mixture of A and B, z A (s À1 ) (eqn (3), Ylilammi et al. 18 ): The average speed of molecules A (i.e. the thermal velocity) is, from the kinetic theory of gases, obtained as (eqn (2) of Ylilammi et al. 18 ): The collision frequency of molecules A in a mixture of A and inert gas I is, from the kinetic theory of gases, obtained as (eqn (1), Ylilammi et al. 18 ): Here, d A and d I are the (hard-sphere model) diameters (m) of molecules A and the inert gas, respectively, and M I is the molar mass of the inert gas (g mol À1 ).¶ The diameters can be estimated for example from the gas-phase viscosity (eqn (7) of Ylilammi et al. 18 ) or the liquid phase density (eqn (8) of Ylilammi et al. 18 ).The Knudsen diffusion coefficient D Kn depends on the microchannel's hydraulic diameter h, the temperature T, and the molar mass of Reactant A, M A (kg mol À1 ) (eqn (4) of Ylilammi et al. 18 ): Instead of solving the differential equation for the partial pressure of Reactant A (eqn (10) of Ylilammi et al. 18 ) numerically, Ylilammi et al. 18 derived an approximate analytic solution of the diffusion equation, which is implemented in this work.With the approximate solution, the partial pressure of Reactant A p A can be analytically calculated for any position x and time t.In the part of the profile where the net adsorption rate g is approximately zero (x o x t , where t stands for ''transition''), the partial pressure p A decreases linearly with x (eqn (18) of Ylilammi et al. 18 ) (Fig. S1, ESI †): and beyond the point x = x t , the decrease is exponential (eqn (24) of Ylilammi et al. 18 ): This journal is © the Owner Societies 2022 In eqn ( 8) and (9), x s , where the linearly extrapolated partial pressure p A is zero, is obtained from (eqn (19), Ylilammi et al. 18 ): Here, D is the apparent longitudinal diffusion coefficient (m 2 s À1 ), which is obtained from (eqn (23), Ylilammi et al. 18 ): Here, q is the adsorption density of the metal M atoms in the growth of the film of the M y Z x material (m À2 ) (i.e., the growth per cycle in ALD, expressed as the areal number density), which can be calculated from the thickness-based growth per cycle (GPC) of the ALD process gpc sat (eqn (9), Ylilammi et al. 18 ): Here, b film is the number of metal atoms in a formula unit of the growing film (e.g., 2 for Al 2 O 3 ), b A is the number of metal atoms in a Reactant A molecule (e.g., 1 for trimethylaluminium), r (kg m À3 ) is the mass density of the ALD film material (composition denoted here as M y Z x ), gpc sat (m) is the ALD GPC (corresponding to saturated reactions) in thickness units, M (kg mol À1 ) is the molar mass of one formula unit of the growing film (M y Z x ), and N 0 is Avogadro's constant (mol À1 ).The transition point x t from eqn ( 8) and ( 9) occurs at (eqn (28), Ylilammi et al. 18 ): Here, c is the sticking probability of Reactant A in collision with the microchannel wall (0 r c r 1, unitless).Q is the wall-collision rate at unit pressure (m À2 s À1 Pa À1 ), calculated from (eqn (14) of Ylilammi et al. 18 ): In this model implementation, the gas-phase diffusion coefficient D A is updated for all positions and times in each cycle, as D A depends on the partial pressure of Reactant A p A (x,t).The apparent longitudinal diffusion coefficient D and the effective diffusion coefficient D eff are also updated accordingly, as they are influenced by the gas-phase diffusion coefficient D A .An illustration of how the partial pressure of Reactant A decreases inside the microchannel is shown in Fig. 2(a).Simulation conditions similar to those of Ylilammi et al. 18 were used.Fig. 2(a) shows a similar trend as Fig. 2 in Ylilammi et al., 18 suggesting that the model was correctly reimplemented.

Langmuir adsorption model and surface coverage
The model is built on the assumption of reversible single-site Langmuir adsorption describing the gas-solid reaction step in ALD: 18 A + * " A*. ( Here, A is the reactant molecule, * is a surface site, and A* denotes a molecule adsorbed on a site.In the Langmuir adsorption model, a surface consisting of a checkerboard can be imagined: all sites are equal, and the adsorbed species are assumed to not interact with each other.If an elementary reaction was assumed, eqn (15) would correspond on the assumption of an associative adsorption mechanism. 3However, it is acknowledged that the actual surface reactions are more complex, 18 and in the model, eqn (15) does not describe an elementary reaction but rather a lumped reaction.The fraction of occupied adsorption sites is called the surface coverage and is denoted with y (0 r y r 1).The fraction of unoccupied or vacant adsorption sites is 1 À y.The rate of adsorption per unit surface area f ads (m À2 s À1 ) is proportional to the fraction of vacant sites, the probability that a collision leads to adsorption c, and the frequency of collisions p A , either as (eqn (11), Ylilammi et al. 18 ): or through the use of the concept of the (gas-phase) collision rate at unit pressure Q: The rate of desorption f des (m À2 s À1 ) depends on the surface concentration of the adsorbed species (yq, m À2 ) and the desorption probability in unit time P d (s À1 ) (eqn (12), Ylilammi et al. 18 ): The net adsorption rate g (m À2 s À1 ) is (eqn (15), Ylilammi et al. 18 ): At equilibrium, the net adsorption rate would be zero, the surface coverage would have reached the equilibrium value y eq , and the equilibrium constant can be defined as (eqn ( 13), Ylilammi et al. 18 ): During adsorption, the ALD reactions are generally not at equilibrium, and the surface coverage y is a function of x and time t.In the model, the surface coverage is solved numerically from the rate equation describing the rate of change of the surface coverage with time (eqn (31), Ylilammi et al. 18 ): The solution requires the partial pressure of Reactant A as a function of position and time, for which eqn ( 8) and ( 9) are used.An illustration of the surface coverage inside the microchannel is shown in Fig. 2(b).This figure shows a similar trend to Fig. 3 by Ylilammi et al., 18 indicating that the model has been correctly re-implemented.

Effect of cycles on film thickness and parameters such as narrowing of the channel
For each cycle, the surface coverage profile is calculated separately, as in eqn (21).The thickness increment caused by the surface coverage is (eqn (37), Ylilammi et al. 18 ): In calculating the thickness profile s(x, N), the thickness increments caused by the N cycles are summed up: In the Ylilammi et al. model, 18 a simplification is made to assume that the free height of the microchannel H is decreased in each cycle by twice the GPC value, as the film grows both on the top and bottom of the microchannel (eqn (35), Ylilammi et al. 18 ): The constant free channel height simplification increases the computational speed. 18The consequence is that the surface coverage for an individual cycle decreases somewhat too steeply in Region III of the thickness profile [see Fig. 1(b)].Ylilammi et al. 18 estimated that the assumption is valid when the film is thin compared to the height of the microchannel H and when the film does not grow much beyond the half-thickness penetration depth x 50% (x p in the Ylilammi et al. 18 model). 18n illustration of the simulated thickness profiles after 1000 cycles in microchannels with various heights is shown in Fig. S2 (ESI †).Here, a similar trend is observed as that in Fig. 4 by Ylilammi et al., 18 confirming that the model has been re-implemented properly.In this implementation, the Knudsen diffusion coefficient D Kn is updated from cycle to cycle as the free height of the microchannel H(N) is updated in each cycle (eqn ( 7)).

Model implementation in MATLAB s s
In this work, ALD thickness profiles in LHAR with different conditions were simulated by implementing the Ylilammi et al. 18 diffusion-reaction model (Model A).The sticking coefficients used for the simulation were compared with those backextracted by the Arts et al. 31 slope method, which is based on the Yanguas-Gil and Elam 20 diffusion-reaction model (Model B).
To discuss the reasons behind the observed differences, the partial pressure of Reactant A and surface coverage in LHAR were simulated also with Model B.
3.1.1.Ylilammi et al. 18 model (Model A).The resulting set of equations for surface coverage y, reactant partial pressure p A , and film thickness s along the microchannel was solved using the software MATLAB s .For discretisation of the geometric domain along the microchannel (x-axis), an equidistant array was used.Based on this, the temporal evolution of the surface coverage y(x,t) (eqn ( 21)) was solved numerically by using MATLAB's ODE23 ordinary differential equation solver with a relative tolerance of 10 À3 and an absolute tolerance of 10 À5 .A simplified flowchart of the algorithm is shown in Fig. 3.The temporal evolution of the reactant partial pressure (p A ) was calculated using eqn (8) and (9).To perform these calculations, it was assumed that p A was zero along the entire microchannel when t = 0.The transport properties (i.e.D eff and D) required for obtaining x s and x t (eqn (10) and ( 13), respectively) were computed for each element along the microchannel using the reactant partial pressure from the previous time step.8  The surface coverage profile y(x,t) was computed by solving eqn (21) until the target pulse time was reached.The film thickness profile s(x) was obtained from eqn (22) while the thickness increment and the updated microchannel height were calculated from eqn (23) and (24), respectively.The This journal is © the Owner Societies 2022 previous procedure was repeated until the defined number of cycles (N) was achieved.
3.1.2.Yanguas-Gil and Elam 20 model (Model B).The slope method reported by Arts et al., 31 which is used to back-extract the value of the (lumped) sticking coefficient, is based on the diffusion-reaction model reported by Yanguas-Gil and Elam. 20his model is similar to Model A, 18 but with two main differences.First, Model B 20 does not use a desorption term to calculate the evolution in y, that is, the adsorption is irreversible.Second, Model B calculates the partial pressure directly from eqn (1), by numerically solving the coupled equations for p A (eqn (1)) and y (eqn ( 21)) simultaneously, while Model A 18 uses an approximate solution (eqn ( 8) and ( 9)).
In this work, and in the work of Arts et al. 31 reporting on the slope method, 31 the coupled equations of the Yanguas-Gil and Elam model 20 were solved assuming free molecular flow (i.e., D eff = D Kn ), using MATLAB's pdepe solver. 43This function solves a system of parabolic and elliptic partial differential equations with one spatial parameter (here, the distance x) and one time parameter t.For the implementation of Model B, the symmetry of the problem was set to 0, corresponding to slab geometry, and default tolerance values of 10 À3 relative tolerance and 10 À6 absolute tolerance were used.Analogous to the implementation of Model A, the parameters x and t were discretised using constant spacing.Finally, the initial conditions p A (x,0) = 0 and y(x,0) = 0 were used, in combination with the boundary conditions p A (0, t) = p A0 and @p A @x ðL; tÞ

Simulation details
Simulations were made with MATLAB s scripts by varying an individual parameter while keeping other parameters constant.
To extract the half-thickness penetration depth, the script chose the first point (x i ,y i ), where x i is equal to or smaller than the half-thickness penetration depth, and then chose another discretisation point (x iÀ1 ,y iÀ1 ), which was one point before (x i ,y i ).Once the two discretisation points were chosen, the half-thickness penetration depth and the slope at the halfthickness penetration depth were interpolated linearly between the two discretisation points (see Fig. S3, ESI †).The total number of discretisation points were selected so that the difference between those two discretisation points in the y-axis is less than or equal to 3% of the whole range.
To compare the simulations made with the Ylilammi et al. model 18 and the Yanguas-Gil and Elam model, 20 and to backextract the sticking coefficient by the slope method, 31 we chose conditions with Kn Z 100 9,25 and Thiele modulus h T 4 1. 27,37 The Knudsen number was calculated as (eqn (1), Cremers et al. 4 ) where l (m) is the mean free path, and h (m) is the hydraulic diameter of the microchannel (eqn (3)).
The mean free path l was calculated as (eqn (3), Cremers et al., 4 and eqn (5.21), 3 Chapman and Cowling 44 ): where k B (m2 kg s À2 K À1 ) is the Boltzmann constant, T (K) is the temperature, p I (Pa) is the partial pressure of the inert gas I, p A0 is the partial pressure of Reactant A at the microchannel entrance (0,t), m I (kg) is the mass of the inert gas molecule I, and s A,I (m 2 ) is the collision cross section between Reactant A and the inert gas I.The collision cross section between molecules i and j is calculated using the following equation (eqn (4), Cremers et al. 4 ): where d i (m) is the hard-sphere diameter of molecule i.**For a first-order reaction with respect to the gas phase species on a LHAR channel geometry, the Thiele modulus 27,37 h T is given by The excess number g, which refers to the amount of Reactant A existing per adsorption site in the LHAR structure, 18 is calculated by using the following equation (eqn (6) Yanguas-Gil and Elam 20 ): where V (m 3 ) and S (m 2 ) are the volume and surface area of the HAR structure, respectively, q (m À2 ) is the adsorption density, and n A is the particle concentration (number density) of Reactant A (m À3 ) at the microchannel entrance (0,t).We simulated thickness profiles at conditions where the excess number g { 1 (e.g. in the baseline condition, was ca.4.5 Â 10 À4 ).Such conditions (g { 1) are required for the slope method to be valid. 31he sticking coefficient was back-extracted with the Arts et al. 31 slope method derived from the Yanguas-Gil and Elam model 20 as follows: where y (À) is the surface coverage and x ˜(À) is the dimensionless distance.In this work, the surface coverage y was extracted from a Type 1 normalised thickness profile expressed as the normalised thickness s/(Ngpc sat ) against the dimensionless distance x ˜.To back-extract the sticking coefficient the total number of discretisation points was selected so that the difference between the two discretisation points chosen in the y-axis was below 1% of the whole range.

ALD in microchannels: general trends with the baseline process
Trends in the evolution of conformality in the microchannels were initially investigated by defining a baseline process with parameters inspired by the experimental TMA-water process. 18,39,40,45For these simulations, the microchannel height H was chosen to be 500 nm, as typically used in microscopic PillarHall TM LHAR structures. 18,31,40,46,47The temperature was chosen to be 250 1C, which is in the typical temperature range of the TMA-water process. 32The partial pressure of Reactant A and the inert gas were chosen as 100 Pa and 500 Pa, respectively, and the Reactant A pulse length was chosen as 0.1 s; these conditions are similar to earlier reported experimental conditions. 40,48The adsorption density of the surface was set to 4 nm À2 , which is in the range observed for the TMA-water process 32,49,50 and in the range typical for ALD. 42The molar mass of Reactant A was set to an arbitrary value of 100 g mol À1 , while that of the purge gas was typical for nitrogen (28 g mol À1 ).
The diameters of Reactant A and the inert gas were 600 and 374 pm, respectively, as in the TMA-water simulation. 18,40The choice of H = 500 nm, combined with the pressure range used, resulted in the Knudsen number Kn for the baseline conditions being 7.6, which is in the transition flow regime (0.1 r Kn r 10). 25The varied parameters are presented in Table 1.
The thickness profiles simulated with varied parameters are shown in Fig. 4(a)-(i).Panels (a) and (b), respectively, show that an increase in the partial pressure of Reactant A, p A , and the pulse time of Reactant A, t 1 , both significantly increase the penetration depth of the film.This result is as expected: the product of p A and t 1 is the dose (Pa s) that defines the penetration depth at free molecular flow conditions (halfthickness penetration depth / ffiffiffiffiffiffiffiffiffi p A t 1 p ). 20,26,31 Panel (c) shows the effect of varying the molar mass of Reactant A, M A .The penetration depth of the film is higher when the molecules are lighter.This is consistent with the fact that the diffusion of light molecules is faster than that of heavy molecules (eqn ( 4) and ( 7)).Panel (d) illustrates the effect of varying the mass density of the grown M y Z x material, r.The lower the density, the higher the grown thickness (note that the penetration depth is not affected).With a lower density, one unit of M y Z x takes up a larger space, so a constant adsorption density q leads to a larger film volume and therefore a larger thickness (eqn (12)).(Note: the thickness-based GPC is not constant in such a case. 3) Panel (e) illustrates the effect of the adsorption density, q (i.e.GPC expressed as the areal number density) on the thickness profile.With other parameters constant, the adsorption density has a strong influence on the growth.9,32 The higher the gpc sat , the higher the film thickness in the saturated region but the lower the penetration depth.This observation is consistent with and explains recent experimental findings where a higher gpc sat resulted in a lower penetration depth. 40,45,51Panel (f) shows how varying the desorption probability, P d , affects the simulation (the Ylilammi et al. 18 model allows reversible reactions).High values of desorption probability affect the shape of the thickness profile especially in Region II, before the Region III of fast decrease (for the regions, see Fig. 1).Panel (g) illustrates the effect of varying the sticking coefficient of Reactant A, c.The sticking coefficient strongly affects the shape of the resulting thickness profile, as already known from earlier simulations made for the diffusionlimited regime. 12,30,31Varying the process temperature T and the inert gas pressure p I has a minor effect on the penetration depth, as seen from panels (h) and (i), respectively.
Earlier works have shown the importance of the sticking coefficient, 11,12,30 as well as the components defining the reactant dose 18,22,25 -i.e., reaction time and reactant pressureon the characteristics of the thickness profile.Simulations made in this work for a typical baseline process resembling the archetypical TMA-water ALD process demonstrated that the process parameters such as the molar mass of the reactant, the adsorption density (derived from the GPC), and the mass density of the film also influence the detailed features of the thickness profile.
This journal is © the Owner Societies 2022  When a film grows into a microchannel, in each ALD cycle, the channel gets narrower from both sides by twice the value of the GPC (eqn ( 24)).The film thickness that completely fills the microchannel is thus half of the microchannel height H.Although an experimental ALD thickness profile can be measured after any number of cycles, the expected shape will depend on the number of ALD cycles, as shown in ref. 40.
How much can a microchannel be filled so that a ''fingerprint'' ALD thickness profile can be measured, whose shape and characteristics are not yet affected by the already grown film?From this fingerprint thickness profile, it is possible to extract the sticking coefficient with the simple slope method. 31reviously, a preferable filling of less than 10% was proposed. 40ere, the effect of the channel filling on the resulting thickness profile was simulated, using the same baseline conditions as in the previous section (Table 1), and varying either the microchannel height H or the number of cycles N.
The results of the simulation series illustrating the effect of channel filling are shown in Fig. 5.In the thickness profiles of panels (a) and (b), the expected features are observed: with a larger microchannel height, the penetration depth increases, and with an increasing number of cycles, the film thickness increases.The scaled thickness profiles of panels (c) and (d) reveal finer trends.With a constant film thickness and varied microchannel height (c), the half-thickness penetration depth x ˜50% (À) first increases with increasing microchannel height, and then starts to decrease [panel (e)].With a constant microchannel height and a varying film thickness of panel (d), the scaled thickness profiles simulated for the smallest cycle numbers (5 to 20) approximately overlap, but already for channel filling of a few percent, the penetration depth starts to decrease with channel filling [panel (e)].Numerical information regarding the half-thickness penetration depth x ˜50% and the slope at this point is presented in Table S1 (ESI †).
The slope at half-thickness penetration depth is shown for both simulation series in Fig. 5(f).For the series where the number of cycles was varied, the slope settles to a constant value (ca.À0.0029) for the smallest amounts of channel filling, as expected.The case where the channel height was varied, shows a different trend: with decreasing channel filling, after a knee, the absolute value of the slope increases again.Table S1 (ESI †) shows the Knudsen number Kn calculated in each case.The knee point occurs at Kn = 8 which is in the transition flow regime. 25Because of the increasing channel height, the Knudsen number decreases with decreasing channel filling (see Fig. S5, ESI †).The reason for the somewhat unexpected trend of the increasing slope with increasing channel height (and decreasing channel filling) was the transition from free molecular flow towards transition flow (Kn o 10), where gas-phase collisions make the diffusion coefficient smaller.
From the simulations made to explore the effect of channel filling, the following can be concluded.(i) To be in the region where the thickness profile is independent of the number of cycles, the channel filling should not exceed a few percent.(ii) The flow regime affects the thickness profile, including the numerical characteristics of the half-thickness penetration depth and the slope at half-thickness penetration depth.To measure a fingerprint thickness profile for an ALD process, the flow condition must be free molecular flow (Kn c 1).To check whether this is the case, the mean free path of the molecules should be calculated (eqn (26)) and compared to the limiting dimension of the feature (eqn ( 25)).

Comparison of thickness profile trends at free molecular flow and transition flow regimes
The simulations in the previous sections revealed that (i) the thickness profile characteristics depend on the flow regime and (ii) the thickness of the grown film affects the characteristics of the thickness profile already from a filling of a few percent.
To compare the trends of the ALD thickness profile in different flow regimes in a well-defined way, we varied individual process parameters to make ALD thickness profiles in free molecular flow (Kn c 1) and transition flow (Kn E 1) conditions. 25,27,28A comparison was made with a single cycle, so that the channel filling does not influence the trends of the thickness profile.Both the scaled thickness profile and the Type 1 normalized thickness profile were used as a basis for comparison.The scaled thickness profile is the most informative thickness profile for an ALD process, and the Type 1 normalized thickness profile is the basis of the slope method. 31he thickness profiles are presented in Fig. S6-S9 (ESI †).For each case, the trends in the half-thickness penetration depth x ˜50% (À) and the absolute value of the slope at x ˜50% were analysed.The numerical values are shown in Fig. S10-S15 (ESI †).
The qualitative thickness profile trends in the free molecular flow and transition flow regimes are summarised in Table 2.In the free molecular flow, the half-thickness penetration depth and the absolute value of the slope remained constant with varying channel heights, as expected.In the transition flow, the penetration depth decreased, and the absolute value of the slope increased with increasing channel heights, most likely resulting from gas-phase collisions.An increase in the reactant partial pressure and pulse time highly increased the penetration depth in both free molecular flow and transition flow, as expected.The desorption probability did not affect the penetration depth and the absolute value of the slope in either flow regime.The penetration depth decreased slightly with the increasing process temperature in both flow regimes.Some process parameters affected the trends of the thickness profile differently in different flow regimes.The molar mass of Reactant A did not affect the absolute value of the slope in free molecular flow while the absolute value of the slope slightly increased with increasing molar mass in the transition flow.The inert gas influenced the thickness profile differently in free molecular flow and transition flow.The inert gas parameters did not affect the penetration depth or the slope in the free molecular flow regime, as expected.In the transition flow regime, the half-thickness penetration depth slightly decreased with increasing pressure and decreasing molar mass of the inert gas.The absolute value of the slope slightly increased with the increasing pressure and decreasing molar mass of inert gas.The absolute value of the slope in the transition flow regime slightly increased with increasing reactant size, while the reactant size did not affect the thickness profile in the free molecular flow regime.
In general, the scaled thickness profile showed the same trends as the Type 1 normalized thickness profile.However, there were two exceptions.With increasing adsorption density q, the absolute value of the slope of the scaled thickness profile markedly increased, while that of the Type 1 normalized thickness profile remained constant.With increasing density of the grown film r, the absolute value of the slope of the scaled thickness profile slightly decreased, while that of the Type 1 normalized thickness profile remained constant.

Comparison of the simulations with different models
Sticking coefficients used for simulations of Model A 18 were compared to those back-extracted from their thickness profiles with the slope method 31 (eqn ( 30)).Different scenarios with varying process temperatures, molar masses, and sticking coefficients were tested, with parameters defined so that the mass transport was always in the free molecular flow regime (Kn Z 100), and the excess number g was {1, 36 as it is where the slope method 31 is valid.Fig. 6(a) shows that when Kn Z 100, the effective diffusion coefficient D eff becomes practically identical to the Knudsen diffusion coefficient D Kn .Note that the comparison was made at conditions where the number of ALD cycles was one.If a larger number of cycles were used and part of the microchannel got filled by the growing film, the slope and penetration depth would have decreased.This channel filling would affect the extracted sticking coefficient and, thus, the comparison.3 lists the parameter values used for simulations and the backextracted sticking coefficients.Fig. S16 (ESI †) shows Type 1 normalized thickness profiles used for the back extraction of (lumped) sticking coefficients.While the order of magnitude is the same, the back-extracted sticking coefficients are systematically ca.25% higher than the set values (see Table 3).Therefore, it seems that the slope method 31 can be used to backextract sticking coefficients from thickness profiles simulated with the current implementation of the Ylilammi et al. 18 model by simply applying a correction factor.
To analyse possible sources of differences in the sticking coefficient values, we compared the partial pressure of Reactant A along the microchannels simulated with the Ylilammi et al. 18 model (Model A) and the Yanguas-Gil and Elam 20 model (Model B); Model B forms the basis of the slope method. 31ig. 7 shows the surface coverage and partial pressure simulated with the above two models against the dimensionless distance., which is the adsorption front where the thickness rapidly decreases. 40A difference compared to Model B is expected, since the Ylilammi et al. 18 model introduced a simplified analytical approximation to the partial pressure (eqn ( 8) and ( 9)) (Fig. S1 and S17, ESI †).We conclude that the more rapid drop of pressure p A at the adsorption front simulated by Model A caused a higher absolute slope value extracted at halfthickness penetration depth, and thus a slightly higher backextracted sticking coefficient compared to the set value.Despite this limitation, Model A is still useful to predict the effect of various process conditions on thickness profile, as shown in this work.

Conclusions and outlook
This work re-implemented the Ylilammi et al. 18 diffusion-reaction model for ALD conformality analysis through thickness profile simulation and used that model to explore trends in the thickness profile inside wide microchannels in different diffusion regimes encountered in reality.
A series of simulations were made to explore the effect on thickness profile characteristics at free molecular flow and transition flow conditions of kinetic and process parameters, such as temperature, (lumped) sticking coefficient, molar mass of the ALD reactant, the reactant's exposure time and pressure, total pressure, density of the grown material, and GPC of the ALD process.Increasing the molar mass and the GPC, for example, resulted in a decreasing penetration depth into the LHAR channel.Trends with the parameter changed depending on the flow regime.To obtain an ALD measurable or a ''fingerprint'' characteristic for a specific ALD process the following conditions should be met: (i) free molecular flow should be the governing mass transport regime, (ii) the channel filling should remain below 5%, and (iii) the scaled thickness profile should be presented, with the dimensionless distance on the   horizontal axis and the thickness divided by cycles on the vertical axis.From such a fingerprint thickness profile, the characteristic GPC is evident, and the kinetic information can be extracted by various means.
The simulations were compared with the recent slope method by back-extracting the sticking coefficient from the ALD thickness profiles at free molecular flow conditions.The slope method gave systematically somewhat higher sticking coefficient values than input values.The difference is most likely related to how, to speed up simulations, the partial pressure of Reactant A inside the channel is analytically approximated in the re-implemented model.
For reactor modelling, kinetic information of real ALD processes is needed.Recent advances have made it possible to measure experimental thickness profiles, which contain the necessary kinetic information, without the need of time-consuming post-preparation of HAR samples.Several theoretical models have been developed to extract (lumped) sticking coefficient parameters from such experimental data.This work has shown that (i) to obtain experimental data for kinetic experiments, detailed knowledge of the experimental conditions, especially pressure, is important, to choose a suitable model for the parameter extraction (most models are based on free molecular flow assumption).Furthermore, (ii) there are differences between the models.The same data fitted with different models may give different results for the extracted fundamental kinetic growth parameters, as the details of the model implementation may affect the results.
For speedy development of the fundamental understanding of ALD processes, and to compare models with each other and with data, it would be advantageous if the scientific ALD community could publish experimental thickness profiles as Open Data and models as Open Code.First such initiatives have already been made: an Open Data community has been initiated in Zenodo.org, 52and the first ALD simulation code has been published in Github. 53The simulation codes of this work are to be published accordingly.

Fig. 1
Fig. 1 Illustration of the ALD film in wide microchannel structures: (a) side and top views of the microchannel with length L and height H, containing a film with thickness s at the channel entrance (illustration intentionally not to scale).(b) Illustration of the different regions (I-IV) of a thickness profile superimposed on a thickness profile with different axes shown.The classifications of thickness profiles and the regions are as in ref. 40, except that here we use the term thickness profile instead of saturation profile.

Fig. 3
Fig. 3 Simplified algorithm for the simulation of thickness profiles with the re-implemented Ylilammi et al. 18 diffusion-reaction model.

Fig. 4
Fig. 4 Illustration of the effect of varying individual parameters on the thickness profile in microchannels, simulated with the Ylilammi et al. model 18 re-implemented in this work.The parameter values used in the simulation are presented in Table 1.Simulations with the baseline values are shown as a solid blue line.The effect of the (a) initial partial pressure of the Reactant A, p A0 , (b) pulse length of Reactant A, t 1 , (c) molecular mass of Reactant A, M A , (d) film density, r, (e) adsorption density, q (i.e.GPC expressed as areal number density), (f) desorption probability, P d , (g) (lumped) sticking coefficient, c, (h) ALD process temperature, T, and (i) inert gas pressure, p I .Note that the image of panel (b) has a larger distance in the horizontal axis than the other images.

Fig. 5
Fig. 5 Illustration of the channel filling effect on the ALD thickness profile in wide microchannels.The parameters used in the simulations are given in Table 1.Baseline simulation results are marked in blue.(a) Thickness profiles simulated with a constant number of cycles of 250 and a varied channel height.(b) Thickness profiles simulated with a constant channel height of 500 nm and a varied number of cycles.(c) The scaled thickness profile from the data of panel (a).(d) The scaled thickness profile from the data of panel (b).(e) The half-thickness penetration depth of the scaled thickness profile as a function of the channel filling fraction, for the data presented in panels (c) and (d).(f) The slope at half-thickness penetration depth as a function of the channel filling fraction (1 À 2s/H), from the data presented in panels (a/c) and (b/d).

Table 2
Summary of the qualitative effects of varying specific parameters on the thickness profile, characterised by the half-thickness penetration depth and the slope at half-thickness penetration depth.The trends are reported separately for different diffusion regimes: free molecular flow (Kn c 1) and transition flow (Kn E 1).Indicators: m increases slightly, mm increases markedly, and mmm increases strongly, -no change, k decreases slightly, kk decreases markedly with increasing parameter values a Simulation parameter (increases) Fraction of reactant pressure of total pressure (p A0 /p) This journal is © the Owner Societies 2022 Phys.Chem.Chem.Phys., 2022, 24, 8645-8660 | 8655

Fig. 6
panels (b)-(d) show the sticking coefficients backextracted by the slope method 31 compared to set values.Table

Fig. 7 (
b) shows observable differences in partial pressures and

a
Parameters used, if not otherwise stated: W = 10 mm, H = 0.1 mm, T = 523.15K, N = 1, t 1 = 2 s, p A0 = 10 Pa, M A = 0.1 kg mol À1 , d A = 6.0 Â 10 À10 m, M I = 0.028 kg mol À1 , d I = 3.74 Â 10 À10 m, p I = 50 Pa, q = 4 nm À2 , r = 3500 kg m À3 , M = 0.050 kg mol À1 , P d = 10 À5 s À1, and c = 0.01.To satisfy the criteria of a difference between the two discretisation points in the y-axis below 1% of the whole range, 5000 discretisation points were used in the simulations with varied molar mass of Reactant A and process temperature while in the simulation with varied sticking coefficient 30 000 discretisation points were used.

Table 3
18icking coefficient values back-extracted (c ext ) by the slope method31against the ones (c) used in simulations implementing the Ylilammi et al. model18with varying process temperatures, molar masses of Reactant A, and sticking coefficients a

b
Number of metal atoms in a reactant molecule in the Ylilammi et al. 18 model (À) c Sticking coefficient (À) c ext Sticking coefficient back-extracted with the slope method 31 (À) D Apparent longitudinal diffusion coefficient in the Ylilammi et al. 18 model (m 2 s À1 ) D A Gas-phase diffusion coefficient of Reactant A (m 2 s À1 ) D eff Effective diffusion coefficient (m 2 s À1 ) D Kn Knudsen diffusion coefficient (m 2 s À1 ) d A Hard-sphere diameter of molecule A (m) d I Hard-sphere diameter of the inert gas molecule (m) f ads Adsorption rate (m À2 s À1 ) f des Desorption rate (m À2 s À1 ) g Net adsorption rate (m À2 s À1 ) gpc sat Saturation growth per cycle, thickness-based, in the Ylilammi et al. 18 model (m) At Partial pressure of Reactant A at x t (Pa) P d Desorption probability in unit time in the Ylilammi et al. 18 model (s À1 ) q Adsorption density of metal M atoms in the growth of film of the M y Z x material (m À2 ) (i.e.GPC expressed as areal number density) Q Collision rate with surface at unit pressure in the Ylilammi et al. 18 model (m À2 s À1 Pa À1 )