Simulation of conformality of ALD growth inside lateral channels: comparison between a diffusion-reaction model and a ballistic transport-reaction model

Atomic layer deposition (ALD) has found significant use in the coating of high-aspect-ratio (HAR) structures. Approaches to model ALD film conformality in HAR structures can generally be classified into diﬀusion–reaction (DR) models, ballistic transport–reaction (BTR) models and Monte Carlo simulations. This work compares saturation profiles obtained using a DR model and a BTR model. The saturation profiles were compared qualitatively and quantitatively in terms of half-coverage penetration depth, slope at half-coverage penetration depth and adsorption front broadness. The results showed qualitative agreement between the models, except for a section of elevated surface coverage at the end of the structure, ‘trunk’, observed in the BTR model. Quantitatively, the BTR model produced deeper penetration into the structure, lower absolute values of the slope at half-coverage penetration depth and broader adsorption fronts compared to the DR model. These differences affect the values obtained when extracting kinetic parameters from the saturation profiles.


Introduction
2][3][4] Other ALD applications involving high-aspect-ratio (HAR) structures include composite materials for energy storage devices 5 and preparation of heterogenous catalysts. 6Furthermore, HAR test structures can be used to study the reaction kinetics of an ALD process. 1,7aturation profiles plot the surface coverage inside a HAR structure as a function of the distance from the structure entrance (i.e.penetration depth). 8,9Parameters, such as the half-coverage penetration depth and the slope at half-coverage penetration depth can be extracted from the saturation profiles to enable quantitative comparison. 8,9In addition, as long as certain assumptions are valid, the slope at half-coverage penetration depth can be used to determine the lumped sticking coefficient, a kinetic parameter of the reaction. 74][25] Additionally, models based on the Boltzmann equation have been developed. 26,27Fig. 1 illustrates particle transport in a DR model and a BTR model.DR models (Fig. 1a), also sometimes referred to as continuum (or continuous) models, 11,20,[28][29][30] solve particle transport using a diffusion equation. 11,28In BTR models (Fig. 1b), particle fluxes are calculated between discretization segments on the structure surface using view factors derived from the structure geometry. 28The assumption of ballistic transport is valid only in Knudsen diffusion conditions, i.e. when the particles interact with the structure and not with other particles in the gas phase. 28A previous work 30 has compared a DR model and BTR model in the context of low-pressure chemical vapor deposition (LPCVD).
The goal of this work was to compare, both qualitatively and quantitatively, two fundamentally different models, a DR model and a BTR model, to describe the evolution of ALD film conformality in microchannels.The comparison is limited to Knudsen diffusion conditions, because this is where one of the models (BTR model) is valid.We will show that, while the model predictions are qualitatively largely similar, quantitative differences exist, which will lead to different interpretations when extracting kinetic information from the saturation profiles.

Model description
In this work, the DR model 12 (DReaM-ALD 31 ) will be referred to as Model A, while the Markov chain based BTR model 10 (Machball 32 ) will be called Model B. Furthermore, this work will refer to the ALD counter-reactants as Reactant A and Reactant B, with Reactant A being the metal precursor of the process, denoted by subscript 'A'.The parameters describing the inert carrier gas will carry the 'I' subscript (in line with ref. 9 and 31), as opposed to 'B' used in the original description of the diffusion-reaction model. 12A brief description of each model will be given next.More in-depth explanations of the diffusion-reaction model can be found in the original work, 12 as well as more recently elsewhere. 9Model B will be introduced in the context of a single reaction pathway with no possibility of reactant recombination.The full description of Model B is presented in the original work. 10Both models assume immediate thermalization of the reactant molecules, i.e. the molecules enter the structure at surface temperature, as opposed to heating up when they reach the structure surface.

Model A: diffusion-reaction model
The DR model by Ylilammi et al. 12 assumes constant concentration of Reactant A along the height and width of the channel.Thus, the one-dimensional diffusion equation can be used to model reactant molecule transport: 12 where p A is the partial pressure of Reactant A (Pa), t is time (s), D eff is the effective diffusion coefficient (m 2 s À1 ), x is the distance from the channel entrance (m), g is the net adsorption rate (m À2 s À1 ), h is the hydraulic diameter of the channel (m), R is the molar gas constant (J K À1 mol À1 ), T is the temperature (K) and N 0 is Avogadro's constant (mol À1 ).The model uses the Bosanquet relation to calculate the effective diffusion coefficient: where D A is the gas-phase diffusion coefficient of Reactant A (m 2 s À1 ) and D Kn is the Knudsen diffusion coefficient (m 2 s À1 ). 33,36In the model, analytical approximations of eqn (1) are derived for two regions inside the channel: a region of linearly decreasing pressure at the beginning of the channel and a region of exponential decrease at the adsorption front. 12or the linear region at the beginning of the channel, it is assumed that the surface becomes instantly saturated (g E 0) and the flow of gas remains constant (i.e.qp A /qt E 0). 12 Thus, eqn (1) reduces to 12 D eff @ 2 p A ðx; tÞ @x 2 % 0; (3) which resolves to 12 p A ðx; Here, p A0 is the partial pressure of Reactant A at the channel entrance (Pa) (i.e.p A (0,t)).In turn, x s is the distance the gas front has travelled at time t (m), approximated from the diffusion length: 12,28 x s ¼ ffiffiffiffiffi ffi Dt p ; ( where D is the apparent longitudal diffusion constant (m 2 s À1 ). 12 The linear approximation in eqn ( 4) is valid up to 12 Here, c is the sticking coefficient (À) and Q is collision rate at unit pressure (m À2 s À1 Pa À1 ): 12 where M A is the molar mass of Reactant A (kg mol À1 ).In the region of exponential decrease in pressure, the partial pressure of Reactant A is given by 12 where 12 This journal is © the Owner Societies 2023 The Langmuir model of adsorption is used to determine how the fractional surface coverage y (À) evolves over time: 12 where q is the adsorption capacity at saturation (nm À2 ).The value of q can be estimated using the experimentally determined saturation growth per cycle gpc sat : 12 where b film is the number of metal atoms present in one formula unit of the deposited film (À), b A is the number of metal atoms in one Reactant A molecule (À), r is the mass density of the film (kg m À3 ) and M film is the molar mass of the deposited film (kg mol À1 ).The net adsorption rate is defined through the rate of adsorption and desorption: 12 Here, the rate of adsorption is: 12 and the desorption rate is: 12 where P d is the desorption probability (s À1 ).Inserting eqn ( 12)-( 14) in eqn (10), we obtain 12 dyðx; tÞ dt ¼ cQp A ðx; tÞ q À cQp A ðx; tÞ q þ P d yðx; tÞ: Alternatively, eqn (15) can be expressed as In Model A, the analytical approximation of the partial pressure of Reactant A (eqn (4) and ( 8)) is inserted in the adsorption model (eqn (15)), which is then solved numerically.In the original work, 12 the numerical solution was obtained using the 4th order Runge-Kutta method.However, the implementation of Model A used in this work, DReaM-ALD, 31 utilized MATLAB's built-in ODE23 differential equation solver instead.Both solvers have been shown to produce similar saturation profiles. 9

Model B: ballistic Markov chain model
The ballistic Markov chain model developed by Yanguas-Gil and Elam 10,32 (Model B) uses a stochastic method to determine reactant molecule transport inside the structure.In ballistic models, reactant molecule transport is expressed in terms of particle fluxes between discretization segments. 10,37The particle flux reaching segment i can be expressed as: 10,37 The particle flux reaching segment i per unit surface area per unit time f i (m À2 s À1 ), multiplied by the segment's surface area S i (m 2 ), is equal to the sum of particles arriving from the other segments in the structure (q ij (1 À b j )S j f j ) and from the structure entrance (q 0i S 0 f 0 ).The particle flux originating from segment j is determined by multiplying the particle flux per surface area arriving at segment j, f j , by the surface area of segment j, S j , and the view factor between segments j and i, q ji (À), as well as the probability that the particle did not react at segment j, (1 À b j ), where b j is the local reaction probability at segment j (À).Similarly, the particle flux originating from the entrance of the structure is calculated by multiplying the particle flux per entry area, f 0 (m À2 s À1 ), by the entry area, S 0 (m 2 ), and by the view factor from the structure entrance to segment i, q 0i (À).The view factors indicate the probability of transmission from a specific segment to another and are determined from the structure geometry. 10,37n principle, the model can be used with structures of any three-dimensional geometry, provided their surfaces can be discretized into segments and the view factors between these segments can be determined.The channel examined in this work was constructed by dividing the channel side walls (upper and lower walls in Fig. 1) into vertical discretization segments, with the entrance and end wall (left and right side in Fig. 1, respectively) represented by a single segment each. 10,32nce the structure geometry has been defined, the model solves reactant molecule transport as a Markov chain. 10In a Markov chain, the probability of an event (i.e. a gas molecule interacting with a specific part of the structure) is only dependent on the outcome of the previous event (i.e.where the molecule previously collided with the structure). 10,38These events can be seen as transitions between states, where each state represents the conditions a molecule is in. 10For example, a molecule could be located at segment i and transition to segment j by colliding with the structure wall. 10Other state transitions include reacting at the segment or exiting the structure. 10In Model B, the probabilities between these state transitions (dg) are defined in a transition probability matrix P dg (À), consisting of the following probabilities (excluding recombination): 10 Here, P ij is the transmission probability matrix (À), P i0 is the reemission probability matrix (À) and P iF i is the reaction probability matrix (À).The relationship between the molecule state and the transition probabilities can be summarized as: 10 p g ðmÞ ¼ Here, the probability of a reactant molecule existing in state g after m collisions, p g (m) (À), is determined by the probability of the reactant molecule existing previously in state d, p d (m À 1), multiplied by the probability of the reactant molecule transitioning from state d to g, P dg .In the case of irreversible adsorption with a single reaction pathway, the evolution of surface coverage depends on the incident flux f i s 0 , multiplied by the local reaction probability b i (y i ): 10 Here, s 0 is the average area of a surface reaction site (m 2 ), which is the inverse of the adsorption capacity at saturation q in Model A: As the adsorption is saturating and irreversible, the local reaction probability at each discretization segment b i (À) is dependent on the specific discretization segment's current surface coverage y i (À): 10 where b 0 is the bare-surface reaction probability (À), similar to c in eqn ( 13).
If we define P i as the probability of a reactant molecule adsorbing in segment i (À), a balance equation can be formulated: 10 Combining eqn ( 20) and ( 23), we obtain: 10 This equation is solved numerically with an adaptive time step, which is based on the change in surface coverage after each step. 10Simulation details

Saturation profile terminology and analysis
The saturation profiles presented in this work follow the classification proposed in ref. 8.The amount of film grown is expressed as surface coverage y, indicated on the vertical axis.On the horizontal axis, the distance from the channel entrance is indicated either as dimensionless distance x ˜= x/H (À) (Type I normalized saturation profile), or normalized distance x = x/L (À) (Type II normalized saturation profile).Here, H is the channel height (m) and L is the channel length (m).
The half-coverage penetration depth x ˜y=0.5 was determined using two points from each saturation profile: the first point from the channel entrance where y r 0.5, and the point prior to it.These points were interpolated to y = 0.5 to obtain the half-coverage penetration depth, as well as the slope at halfcoverage penetration depth.The method has been described elsewhere 9 in detail.Absolute values were calculated for the slope to ease comparison, as the slope is typically negative.
As a further point of comparison, the broadness of the adsorption front, x ˜y=0.1 À x ˜y=0.9 , was measured.The adsorption front broadness was defined as the dimensionless distance between the first point where y o 0.9 and the first point where y o 0.1, measured from the channel entrance.Saturation profiles, where the surface coverage was below one at the channel entrance, or where the surface coverage did not decrease below 0.1, were excluded from the comparison.

Selection of parameters
The trimethylaluminum (TMA)-water ALD process to grow aluminum oxide 8,39 was used as a starting point for the selection of the simulation parameters, similarly as in ref. 9 and 12.The molar mass of Reactant A M A was selected to be 0.1 kg mol À1 , which is in the same order of magnitude as TMA.The molar mass of the inert carrier gas M I was selected based on the commonly used nitrogen as 0.028 kg mol À1 .Similarly, the molecular diameters d A = 0.6 nm and d I = 0.4 nm were chosen based on approximations 12 of TMA and nitrogen gas, respectively.The molar mass of the deposited film M film was 0.05 kg mol À1 and the film density r was 3500 kg m À3 , which are in the same order of magnitude as the respective values for aluminum oxide. 40The number of metal atoms in a reactant molecule b A and molecular unit of film b film were set to one.The saturated adsorption capacity q was 4 nm À2 , which is in the same order of magnitude as indicated by the growth-per-cycle observed in the TMA-water process. 39,41An equivalent value of 0.25 nm 2 for the average surface site area s 0 was determined using eqn (21).The temperature T was selected to be in a range typical for ALD, 11 with a baseline value of 573.15 K.
A channel height H of 0.5 mm was selected based on the dimensions commonly used in PillarHallt lateral HAR test structures. 8,12,42,43The channel width W was arbitrarily set to 0.01 m, to be orders of magnitude larger than the channel height.The length of the channel L was selected to be 500 mm, which gives a sufficiently high aspect ratio (AR = L/H) of 1000 to obtain visible adsorption fronts for most of the simulations.
A proper resolution for the simulations was determined through preliminary simulations, where the number of discretization points f p (Model A) or discretization segments f s (Model B) was varied.To account for varying channel lengths, the resolution was defined relative to the channel aspect ratio.For example, the aspect ratio of a channel of L = 500 mm and H = 0.5 mm is (500 mm)/ (0.5 mm) = 1000 (À).When the relative number of discretization points f p is four, the number of discretization points used is 4 Â 1000 = 4000 (À).A summary of the studied resolutions is given in Table 1.
A range of values was selected for the sticking coefficient c to cover the upper and lower extremes (0.00001 to 1), with a baseline value of 0.001.Similarly, the exposure duration t end ranged through multiple orders of magnitude: 0.01 to 100 s, with a baseline value of 1 s.
As Model B is only applicable in Knudsen diffusion conditions (Kn c 1), the initial partial pressure of Reactant A p A0 had to be sufficiently low.Therefore, the initial partial pressure of This journal is © the Owner Societies 2023 are summarized in Table 2.The selected simulation parameters resulted in Knudsen numbers ranging from B99 to B62000 and the Thiele modulus ranged from B2.7 to B870.The Knudsen number and Thiele modulus were calculated as described in ref. 9.

Effect of the simulation resolution
The number of discretization points (Model A) and discretization segments (Model B) was varied to determine a suitable resolution for the simulations.Fig. 2 shows how the resolution affects the resulting saturation profile in channels with different lengths.In Model A, increasing the number of discretization points merely increases the saturation profile roundness.In contrast, the saturation profiles generated using Model B show a change in shape and penetration depth when the number of discretization segments is changed: as the number of discretization segments is increased, the surface coverage at the end of the channel is reduced and the penetration depth decreases.Additionally, the adsorption front becomes more steep.Once the relative number of discretization segments is increased beyond one (i.e. one discretization segment per dimensionless distance unit of channel length), the saturation profiles stabilize.For the comparison of the saturation profiles produced by the models, an intentionally high resolution was selected to ensure consistent results.The relative number of discretization points was selected to be two for the simulations performed using Model A. For Model B, the relative number of discretization segments was selected to be four.

Comparison of saturation profiles
The series of simulated saturation profiles is presented in Fig. 3, while the half-coverage penetration depth and the absolute value of the slope at half-coverage penetration depth extracted from each saturation profile are shown in Fig. 4. When the temperature was increased (Fig. 3a and 4a), the penetration depth of the film decreased, while the shape of the saturation profile remained similar.An increase in initial partial pressure of Reactant A (Fig. 3b and 4b) caused the film to penetrate deeper into the structure and, in Model B, led to an increase in coverage at the end of the channel.Increasing the molar mass of Reactant A (Fig. 3c and 4c) decreased the penetration depth, with the decrease in penetration depth tapering off as the molar mass was increased.The higher the aspect ratio of the channel (Fig. 3d), the lower the penetration depth was, relative to the total length of the channel.Longer exposure durations (Fig. 3e and 4e) led to deeper penetration depths in both models, compared to shorter exposure durations.In Model B, the channel was covered with a uniform film at an exposure duration of t = 10 s, which was not sufficient to coat the entire channel in Model A. In both models, increasing the sticking coefficient (Fig. 3f and 4f) resulted in an increase in the absolute value of the slope at half-coverage penetration depth.Decreasing the adsorption capacity at saturation (Model A) or increasing the average surface site area (Model B) (Fig. 3g and 4g) generated deeper penetration and increased the coverage at the end of the channel in Model B. The increase in penetration depth with decreasing adsorption capacity at saturation observed in Model A is consistent with a previous work. 9Finally, varying the initial partial pressure of Reactant A and the exposure duration, while keeping the exposure constant (Fig. 3h), had no impact on the resulting saturation profile in either model.Fig. 5 shows the adsorption front broadness measured from the series of saturation profiles.Increasing the temperature (Fig. 5a) had no effect in Model A, while the adsorption front broadness decreased slightly in Model B. An increase in initial partial pressure of Reactant A (Fig. 5b) led to no change in  Model A, while an increase in adsorption front broadness was observed in Model B. When the molar mass of Reactant A was increased (Fig. 5c), the adsorption front broadness was not affected in Model A. In Model B, the adsorption front broadness decreased and appeared to plateau.The aspect ratio of the channel (Fig. 5d) did not affect the adsorption front broadness in either model.The variation in exposure duration (Fig. 5e) produced few saturation profiles, from which the adsorption front broadness could be determined.Between the two valid saturation profiles obtained from Model A, no change in adsorption front broadness could be observed.Increasing the sticking coefficient (Fig. 5f) caused the adsorption front broadness to decrease sharply and then plateau in both models.
When the average surface site area was increased (Model B) (Fig. 5g), i.e. the adsorption capacity at saturation was decreased (Model A), the adsorption front broadness remained unaffected in Model A, while it increased in Model B. The series of constant exposure by variation of the initial partial pressure of Reactant A and the exposure duration (Fig. 5h) produced identical saturation profiles in both models.Fig. 6 compares the half-coverage penetration depth, the absolute value of the slope at half-coverage penetration depth and the adsorption front broadness observed throughout the series of simulations.Model B yielded systematically larger penetration depths (Fig. 6a) and broader adsorption fronts (Fig. 6c) than Model A. The highest deviation in half-coverage penetration depth between the models was observed when the average surface site area was increased (Fig. 6a and 4g), while the models were closest at low partial pressures of Reactant A (Fig. 6a and 4b).Varying most of the simulation parameters had little impact on the slope at halfcoverage penetration depth, with the exception of the sticking coefficient, which led to an increase in the absolute value of the slope as it was increased (Fig. 6b and 4f).Similarly, the broadness of the adsorption front was influenced most by the sticking coefficient in both models (Fig. 6c and 5f).In Model A, the other varied parameters had no impact on adsorption front broadness (Fig. 6c).In contrast, Model B showed variation in adsorption front broadness in all cases, especially when the partial pressure of Reactant A and average surface site area were varied (Fig. 6c, 5b and 5g).

Discussion
This section comments on the general trends observed in the results and provides reasoning for the resolution-dependent behaviour of Model B. Additionally, the inclusion of the end wall in the structure geometry and its effect on the shape of the simulated saturation profile is discussed.

Similarities and differences between the models' output
Model A and Model B produced qualitatively similar trends in terms of half-coverage penetration depth when the simulation  2. Table 1 shows the absolute number of discretization points or segments used in the simulations.
This journal is © the Owner Societies 2023 parameters were varied (Fig. 3 and 4).For example, an increase in temperature resulted in a slight decrease in half-coverage penetration depth (Fig. 3a and 4a).However, Model B showed systematically deeper penetration into the structure than Model A (Fig. 4 and 6).The relative difference in half-coverage penetration depth ranged from B38% higher in Model B compared to Model A (p A = 1 Pa, Fig. 4b and c = 0.00001, Fig. 4f) to B57% higher in Model B compared to Model A (s 0 = 1 nm 2 , Fig. 4g).On average, the half-coverage penetration depth was B48% higher in Model B compared to Model A. This systematic difference is consistent with an earlier comparison 30 of a DR model and BTR model, which operated in the context of LPCVD and studied the relationship between the film thickness at the channel end and the outer, exposed surface (i.e. the step coverage).In the earlier comparison, 30 the BTR model gave a systematically higher step coverage compared to the DR model.
On average, the absolute value of the slope at half-coverage penetration depth was B2.0 times higher in Model A compared to Model B. In the series of simulated saturation profiles, the most agreement between Model A and Model B, in terms of slope, was at p A = 1 Pa (Fig. 4b), where the absolute value of the slope was B1.6 times higher in Model A compared to Model B. The models agreed the least at c = 1 (Fig. 4f), where the absolute value of the slope at half-coverage penetration depth was B4.2 times higher in Model A compared to Model B. Overall, varying the sticking coefficient had the most extreme effect on the slope at half-coverage penetration depth for both models.When the sticking coefficient was increased (Fig. 3f), the absolute value of the slope at halfcoverage penetration depth increased (Fig. 4f), with the saturation profiles produced by both models pivoting around a fixed point.Varying the sticking coefficient (Fig. 4f), produced a broader range in the absolute value of the slope in Model A (B0.0031 to B0.30) compared to Model B (B0.0051 to B0.071).
The broadness of the adsorption front remained mostly constant throughout the simulations in Model A, while Model B produced more variation in broadness (Fig. 5).The adsorption front broadness in Model B was, on average, B2.8 times the respective value given by Model A (Fig. 5).The series where the Fig. 3 Series of saturation profiles that shows the effect of (a) temperature T (K), (b) initial partial pressure of Reactant A p A0 (Pa), (c) molar mass of Reactant A M A (kg mol À1 ), (d) channel aspect ratio L/H (À), (e) exposure duration t end (s), (f) sticking coefficient c (À), (g) surface site area s 0 (nm 2 ) and adsorption capacity at saturation q (nm À2 ), and (h) constant dose with varying initial partial pressure of Reactant A p A0 and exposure duration t end .In each panel, the left graph shows results for the diffusion-reaction model (Model A), while the right graph shows the saturation profiles obtained using the ballistic transport-reaction model (Model B).See Table 2 for a summary of the simulation parameters, the corresponding symbols used in each model, as well as the variation of the adsorption capacity at saturation q in panel (g).Note that panel (d) has normalized distance on the horizontal axis, as opposed to dimensionless distance in the other panels.2. Saturation profiles, where the surface coverage at the channel entrance was less than one, or where the surface coverage was equal to one throughout the channel, were excluded from parameter extraction.Fig. 5 Comparison of the adsorption front broadness measured from each series of simulations as a function of the varied simulation parameter, indicated on the horizontal axis of each graph.The adsorption front broadness was determined by measuring the distance between the first point from the channel entrance where the surface coverage y o 0.9 and the first point where y o 0.1.If the surface coverage was below one at the beginning of the channel or did not fall below 0.1 in the rest of the channel, the saturation profile was excluded from the broadness measurements.The simulation parameters are summarized in Table 2.
This journal is © the Owner Societies 2023 sticking coefficient was varied produced the greatest difference in adsorption front broadness between Model A and Model B: when c = 1, the adsorption front broadness in Model B was equal to B19 times the respective value in Model A (Fig. 5f).If the sticking coefficient series is excluded from the examined saturation profiles, the adsorption front broadness in Model B is, on average, B2.0 times higher than the value obtained from Model A. As varying the sticking coefficient also had the highest impact on the absolute value of the slope at half-coverage penetration depth, the results suggest that the measured slope and the broadness of the adsorption front correlate to some degree.
While the saturation profiles produced by both models were generally similar in shape, Model B produced a section of increasing coverage at the closed end of the channel, which we call a trunk in this work.The shape of the saturation profiles was especially sensitive to variation of the sticking coefficient, with the adsorption front turning steeper and approaching the shape of a step function as the sticking coefficient was increased (Fig. 3f).In Model B, increasing the sticking coefficient also increased the height of the trunk.Furthermore, in Model B, the whole trunk section (i.e. from the saturation profile's minimum coverage to the end of the channel) was elevated as the adsorption front approached the end of the channel.

Resolution and end wall effect in Model B
In Model B, the number of discretization segments was found to affect the saturation profiles.Increasing the number of discretization segments decreased the penetration depth of the profiles, up to a ratio of approximately one discretization segment per aspect ratio unit, beyond which the profiles stabilized (Fig. 2b and d).
This change in penetration depth is explained by how Model B treats particle transport as probabilistic transitions between discretized segments of the channel.In Model B, the probability distribution of the potential target discretization segments for a given molecule to travel to is determined by the view factors between the departure segment and each target discretization segment. 10,37These view factors follow a cosine distribution with the maximum probability of transmission being along the surface normal n (À) of the departure segment. 10,37ig. 7 illustrates how the angle between the surface normal and the connecting line between neighboring discretization segments O (À) (i.e. the normal angle) behaves as the channel aspect ratio and discretization resolution are varied.As the aspect ratio of the channel is increased (Fig. 7b), or the number of segments is decreased, the normal angle increases.An increase in the normal angle, in turn, decreases the view factor between the segments, causing the view factor distribution of a given channel wall segment to narrow.Conversely, increasing the resolution (Fig. 7c) broadens the view factor distribution, which leads to an increase in the number of potential target discretization segments for the molecule to travel to.A visual representation of the view factors can be found in the ESI † (Section 2).Furthermore, when the resolution is low, a molecule advancing the length of a single discretization segment has a more pronounced effect on the saturation profile, compared to when the resolution is high, as the resolution directly affects the length of each discretization segment on the channel wall.
With a higher resolution and a higher number of potential target discretization segments, the molecule has to interact with a higher number of segments to penetrate the structure to a given depth, compared to a lower resolution with fewer potential target discretization segments from each channel wall segment.Given that each collision with a bare or partially covered segment entails the possibility of adsorption, a high number of segments will curb the probability of a molecule reaching a given penetration depth.A high number of collisions with bare or partially covered segments (i.e. a broad view factor distribution) also produces narrower adsorption fronts, as the propagation of gas molecules in the channel is more dependent on the local degree of surface coverage, compared to a case with few collisions (i.e. a narrow view factor distribution).
Model B produced a region of increasing surface coverage at the channel end, a trunk, which was not observed in Model A. This is explained by the fact that the end wall of the channel is included in the modelled structure in Model B, while in Model A, it is not.As the molecules approach the closed end of the channel, the normal angle towards the end wall decreases, which increases its probability as a destination segment.The view factor of the end wall is further increased by its larger solid angle, 37,44 in comparison to the channel side wall segments. 10s a result of the relatively high view factor of the end wall, reactant molecules are likely to collide with it, and either adsorb directly on the end wall or continue to the discretization segments in its vincinity.As the molecules are not able to propagate past the end wall, they cumulate in the channel end area, forming the trunk shape.A trunk shape has previously been produced by a Monte Carlo model, 45 in the context of plasma-assisted ALD.However, to the best of our knowledge, the trunk has not been reported experimentally yet.Further examination of the trunk in channel and hole structures is presented in the ESI † (Section 1).

Conclusions
This work provided a qualitative and quantitative comparison of a diffusion-reaction model (DR model, Model A) and a ballistic transport-reaction model (BTR model, Model B) for ALD conformality simulation.The models were used to generate series of saturation profiles based on identical simulation parameters, with one simulation parameter varied in each series.The parameter values were chosen so that the Knudsen number was Kn c 1, i.e. free molecular flow, as this is where the BTR model is valid.These saturation profiles were compared qualitatively, and quantitatively in terms of half-coverage penetration depth, slope at half-coverage penetration depth and adsorption front broadness.
Qualitatively, the DR model and the BTR model resulted in similar trends of ALD film conformality in each series of parameter variation.For example, in both models, increasing the partial pressure of Reactant A caused the film to penetrate further into the structure, while an increase in the Reactant A molar mass decreased penetration depth.Increasing the temperature caused the half-coverage penetration depth to decrease.Furthermore, an increase in the sticking coefficient resulted in steeper adsorption fronts in both models.However, the BTR model produced an elevated section at the channel end, i.e. a trunk, which was not present in the DR model.This lack of a trunk in the DR model is explained by the fact that, in contrast to the BTR model, the DR model does not account for the end wall of the channel.The trunk observed in the BTR model is likely the result of Reactant A molecules accumulating at the channel end.
Quantitatively, the half-coverage penetration depth, as well as the slope at half-coverage penetration depth, were extracted from the saturation profiles and compared.The half-coverage penetration depth was systematically higher in the BTR model---on average, by B48% compared to the DR model.In contrast, the absolute value of the slope at half-coverage penetration depth was consistently higher in the DR model, reaching B1.6 to B4.2 times the respective value in the BTR model.These gentler slopes Fig. 7 Schematics depicting the normal angle in microchannels of varying aspect ratios discretized with varying resolutions, as well as the resulting saturation profiles.Three situations are shown: (a) reference case, (b) case with higher aspect ratio and the same number of segments as in the reference case and (c) case with higher aspect ratio and the same resolution as in the reference case.The left column shows schematics depicting how the channel aspect ratio and discretization resolution affect the normal angle O between neighboring discretization segments (i and i + 1, highlighted in orange).The blue arrow marks the trajectory of the molecule, while the surface normal n is marked by a black arrow.In the right column, the resulting saturation profiles (orange dotted line in panels (b) and (c)) are illustrated in comparison to the reference case (grey solid line).If the aspect ratio is increased without increasing the number of discretization segments (panel (b)), the resolution decreases and the normal angle increases, compared to the reference case (panel (a)): O 1 o O 2 .As a result, the resulting saturation profile (orange dotted line in panel (b)) does not align with the reference case (grey solid line).If the aspect ratio is increased and the resolution is maintained (panel (c)), the normal angle remains similar to the reference case (panel (a)).Here, the resulting saturation profile (orange dotted line in panel (c)) aligns with the reference case (grey solid line).The saturation profiles were produced using the baseline parameters shown in Table 2, except (a) f s = 0.250, (b) f s = 0.125, L = 1000 mm and (c) f s = 0.250, L = 1000 mm.An intentionally low resolution was used to make the difference in saturation profiles more clearly visible.Note that, as the channel aspect ratio is varied, the length of the saturation profile in dimensionless distance units varies as well.Therefore, the trunk is not visible in the longer saturation profiles (orange dotted line in panels (b) and (c)).
This journal is © the Owner Societies 2023 observed in the BTR model are connected to the broader adsorption fronts it produces.The adsorption fronts obtained from the BTR model were generally approximately twice as broad, compared to the DR model, if the sticking coefficient is constant.When the sticking coefficient was varied, the adsorption front broadness observed in the BTR model exceeded the DR model by a factor of up to B19.
In the BTR model, the number of discretization segments was found to significantly influence the shape and penetration depth of the saturation profile, depending on the aspect ratio of the channel.When the resolution (i.e. the number of discretization segments per channel aspect ratio unit) was increased, the penetration depth decreased.Beyond a resolution of one discretization segment per aspect ratio unit, the saturation profiles stabilized.This behaviour is explained by the geometry involved in the calculation of the view factors between discretization segments, which determine molecule transport inside the channel.The authors recommend to allocate at least two discretization segments per channel aspect ratio unit in future simulations with the BTR model.
This work showed that the choice of model has a significant impact on the quantitative characteristics obtained in conformality simulations.When fitting experimental data, the backextracted kinetic parameters will be influenced by the selected model.For example, when the steepness of the adsorption front is used to determine the sticking coefficient of the process, 7 the two models used in this work would lead to different results.

Fig. 1
Fig. 1 Simplified schematic depicting the core principles of reactant molecule transport in (a) a diffusion-reaction (DR) model (Model A) and (b) a ballistic transport-reaction (BTR) model (ModelB).In the DR model, the partial pressure of Reactant A p A (x,t) (indicated by a blue gradient) is determined by the one-dimensional diffusion equation.In the BTR model, the structure is discretized along the channel wall into segments i,jA(0,1,. ..N), with the first segment being the entrance of the channel 0, and the last segment N being the channel end wall (example segments marked in orange).Reactant molecule transport is modelled as transitions between particle states (for example, from discretization segment i to N) and the transition probabilities are given by a probability matrix P dg , which is based on the structure geometry and local reaction probability.

a
The following simulation parameters are denoted by different symbols compared to when Model A was first presented:12 p I was denoted as p B , M I was denoted as M B , M film was denoted as M, t end was denoted as t p and d I was denoted as d B .b In Model B, the channel height H and length L are considered only through the aspect ratio AR = L/H.c L was varied in the preliminary simulations to find a suitable resolution as indicated in Table1.dThe adsorption capacity at saturation q and the average surface site area s 0 are inversely related (eqn(21)).e In this work, f p and f s are defined as the ratio of the number of discretization points or segments, respectively, to the aspect ratio of the channel.As such, they should not be confused with the usage in the Machball model,10 where f denotes particle flux per unit surface area per time.This journal is © the Owner Societies 2023 Phys.Chem.Chem.Phys., 2023, 25, 22952-22964 | 22957

Fig. 2
Fig. 2 Series of saturation profiles obtained with varying numbers of (a) and (c) discretization points (Model A) and (b) and (d) discretization segments (Model B).The legends indicate the resolution as the number of discretization points or segments relative to the channel length in dimensionless distance units.In panel (a), a set of insets shows the shape of the saturation profile in the highlighted area for two relative numbers of discretization points: (upper inset) 0.0625 and (lower inset) 1.00.The channel length L was (a) and (b) 500 mm and (c) and (d) 1000 mm.The other simulation parameters were as in the baseline values indicated in Table2.Table1shows the absolute number of discretization points or segments used in the simulations.

Fig. 4
Fig. 4 Comparison of the parameters extracted from each series of simulations as a function of the varied simulation parameter.Each panel shows how (left) the half-coverage penetration depth x ˜y=0.5 and (right) the absolute value of the slope at half-coverage penetration depth (Dy/Dx ˜)y=0.5 changes as the simulation parameter indicated on the horizontal axis is varied.The simulation parameters are summarized in Table2.Saturation profiles, where the surface coverage at the channel entrance was less than one, or where the surface coverage was equal to one throughout the channel, were excluded from parameter extraction.

Fig. 6
Fig. 6 Parity plots comparing the (a) half-coverage penetration depth x ˜y=0.5 , (b) absolute value of slope at half-coverage penetration depth (Dy/Dx ˜)y=0.5 and (c) adsorption front broadness x ˜y=0.1 À x ˜y=0.9 observed in Model A and Model B. The legend shows which simulation parameter was varied for each plotted value.On the grey dashed line, the results from both models would be identical.The simulation parameters are summarized in Table2.

Table 1
Relative number of discretization points or segments and the corresponding number of discretization points or segments in channels with different lengths

Table 2
Values of the varied and constant simulation parameters, with the baseline value of each varied parameter marked in bold a t c 0.01, 0.1, 1, 10, 100 s VI c b 0 1, 0.1, 0.01, 0.001, 0.0001, 0.
This journal is © the Owner Societies 2023 Phys.Chem.Chem.Phys., 2023, 25, 22952-22964 | 22963 y i Fractional surface coverage at segment i (À) x Normalized distance (À) r Mass density of film (kg m À3 ) f i Particle flux reaching segment i per unit surface area per unit time (m À2 s À1 )