Atomistic Simulations of Plasma-Enhanced Atomic Layer Deposition

Plasma-enhanced atomic layer deposition (PEALD) is a widely used, powerful layer-by-layer coating technology. Here, we present an atomistic simulation scheme for PEALD processes, combining the Monte Carlo deposition algorithm and structure relaxation using molecular dynamics. In contrast to previous implementations, our approach employs a real, atomistic model of the precursor. This allows us to account for steric hindrance and overlap restrictions at the surface corresponding to the real precursor deposition step. In addition, our scheme takes various process parameters into account, employing predefined probabilities for precursor products at each Monte Carlo deposition step. The new simulation protocol was applied to investigate PEALD synthesis of SiO2 thin films using the bis-diethylaminosilane precursor. It revealed that increasing the probability for precursor binding to one surface oxygen atom favors amorphous layer growth, a large number of –OH impurities, and the formation of voids. In contrast, a higher probability for precursor binding to two surface oxygen atoms leads to dense SiO2 film growth and a reduction of –OH impurities. Increasing the probability for the formation of doubly bonded precursor sites is therefore the key factor for the formation of dense SiO2 PEALD thin films with reduced amounts of voids and –OH impurities.


Introduction
Atomic layer deposition (ALD) is a coating technology based on self-terminating reactions of gaseous precursor molecules with surface functional groups of the substrate. As shown in Figure 1, a precursor vapor (usually an organometallic compound) is introduced in the chemical reactor. It reacts with the functional groups at the surface of the substrate, until the surface is saturated (first half-cycle). Excess precursor and reaction byproducts are purged by inert gases. Next, a second precursor is pulsed in the reactor, such as an oxidizing agent like H 2 O, ozone, O 2 , or O 2 plasma, which reacts with the remaining ligands of the precursor material (second half-cycle). The reactor is purged again. The ALD process is continued by repeating the above pulsing and purging steps, which constitute the ALD cycle, until the desired film thickness is reached. These reactions depend only on the abundancy and variety of the surface functional groups for a specific precursor, and not on the shape of the substrate. Hence, the ALD technology leads to extremely conformal coatings on highly curved [1][2][3][4][5] or nanostructured substrates. Additionally, the composition of the thin films can be controlled with atomic precision by alternating the organometallic compound for binary or more complex mixtures [1,5]. Plasma-enhanced ALD (PEALD), which utilizes oxygen plasma, is applied in order to increase the reactivity of the oxidizing agent [6]. This enables, for example, room temperature deposition of SiO 2 thin films, which is possible only at 200 • C in typical thermal processes using O 3 [7]. However, recently developed highly reactive precursors allow for thermal ALD processes at temperatures approaching room temperature [8,9]. The O radical species present in the plasma allow for additional tailoring of promising trend in the development of simulation methods for nanomaterials is the use of so-called reactive force fields, such as the ReaxFF approach of van Duin and co-workers [27]. This has been applied to a wide range of materials including amorphous and crystalline SiO 2 [28,29] and Al 2 O 3 [30], yielding very good agreement with experimental data. Although parameterization of ReaxFF is relatively complex, it is very powerful since it allows for an accurate description of the chemistry of nanostructures with a computational cost far below that of quantum mechanical methods. In our recent study, we demonstrated that the structure of amorphous SiO 2 could be accurately described at the ReaxFF level by employing relatively small semi-amorphous periodic models [31]. A simulation cell containing only 64 SiO 2 units yielded a structure that properly reproduced experimental structural parameters, mechanical properties, and IR spectra of bulk silica, using subsequent refinement at the DFT level.
In this article, we present an atomistic simulation scheme for a PEALD process that combines the Monte Carlo (MC) deposition algorithm and structure relaxation using MD. Our scheme is based on the one proposed by Hu et al. [15]. However, in contrast to their implementation that directly deposits hydroxylated products of the second ALD half-cycle, we use a real, atomistic model of the precursor deposited in the first half-cycle. In addition, for the MC deposition step, our approach introduces site occupation probabilities for precursor reaction products. The performance of our simulation scheme is demonstrated for PEALD synthesis of SiO 2 thin films, using a bis-diethylaminosilane (BDEAS) precursor with the structural formula H 2 Si(NEt 2 ) 2 , Et = C 2 H 5 . We show that the resulting structure of the SiO 2 coating strongly depends on the occupation probabilities for precursor products. Densification of the SiO 2 film and reduction of -OH impurities are observed for increasing occupation probability of doubly coordinated surface sites. and promising trend in the development of simulation methods for nanomaterials is the use of so-called reactive force fields, such as the ReaxFF approach of van Duin and co-workers [27]. This has been applied to a wide range of materials including amorphous and crystalline SiO2 [28,29] and Al2O3 [30], yielding very good agreement with experimental data. Although parameterization of ReaxFF is relatively complex, it is very powerful since it allows for an accurate description of the chemistry of nanostructures with a computational cost far below that of quantum mechanical methods. In our recent study, we demonstrated that the structure of amorphous SiO2 could be accurately described at the ReaxFF level by employing relatively small semi-amorphous periodic models [31]. A simulation cell containing only 64 SiO2 units yielded a structure that properly reproduced experimental structural parameters, mechanical properties, and IR spectra of bulk silica, using subsequent refinement at the DFT level.
In this article, we present an atomistic simulation scheme for a PEALD process that combines the Monte Carlo (MC) deposition algorithm and structure relaxation using MD. Our scheme is based on the one proposed by Hu et al. [15]. However, in contrast to their implementation that directly deposits hydroxylated products of the second ALD half-cycle, we use a real, atomistic model of the precursor deposited in the first half-cycle. In addition, for the MC deposition step, our approach introduces site occupation probabilities for precursor reaction products. The performance of our simulation scheme is demonstrated for PEALD synthesis of SiO2 thin films, using a bis-diethylaminosilane (BDEAS) precursor with the structural formula H2Si(NEt2)2, Et = C2H5. We show that the resulting structure of the SiO2 coating strongly depends on the occupation probabilities for precursor products. Densification of the SiO2 film and reduction of -OH impurities are observed for increasing occupation probability of doubly coordinated surface sites.  Figure 2 shows the implemented protocol for simulation of the first PEALD deposition half-cycle. In contrast to the original scheme of Hu et al. [15], it used the MC method to simulate the full precursor deposition process. In the first step of each precursor deposition half-cycle, the surface -OH groups were enumerated. Then, during each MC step, a free surface -OH group was randomly selected. The surface was scanned for neighboring free -OH groups and, depending on their number and precursor type, a MC move was performed that creates one of possible precursor products. This  Figure 2 shows the implemented protocol for simulation of the first PEALD deposition half-cycle. In contrast to the original scheme of Hu et al. [15], it used the MC method to simulate the full precursor deposition process. In the first step of each precursor deposition half-cycle, the surface -OH groups were enumerated. Then, during each MC step, a free surface -OH group was randomly selected. The surface was scanned for neighboring free -OH groups and, depending on their number and precursor type, a MC move was performed that creates one of possible precursor products. This part is also different than in the original scheme of Hu et al. [15], where direct deposition of the product of the second half-cycle (i.e., -Al(OH) 2 ) was performed. Steric hindrance and overlap restrictions at the surface corresponding to the real precursor deposition step were accounted for only in an approximate way, by defining an exclusion zone around each successfully deposited -Al(OH) 2 group. In contrast, our approach employs a real, atomistic model of the precursor. After each attachment step, the structure was relaxed, either by local structure optimization or by short simulated annealing using molecular dynamics (MD). The energy of the chemisorption was calculated and used to calculate the acceptance criterion employing the usual Metropolis-Hastings rule [32]. The procedure was repeated until all enumerated surface -OH groups were either reacted or excluded from MC moves due to steric hindrance and overlap with chemisorbed precursor molecules. Alternatively, a predefined concentration of surface -OH groups can be left unreacted. The final precursor monolayer structure was relaxed and equilibrated using MD. The second ALD half-cycle was modeled by simply replacing all organic surface groups with -OH groups, followed by relaxation and equilibration using MD. This corresponds to full conversion of the products of the precursor pulse to surface hydroxyl species. Such an approach is particularly well suited for simulating ALD growth of systems with primary reactions that have relatively low activation energies [14]. It has been shown that, unlike thermal ALD, PEALD can be efficiently performed even at room temperature [6,7], justifying the assumption of low activation energies. It has also been shown by calculations that oxygen plasma can quantitatively oxidize organic surface precursor species [17].

PEALD Simulation Procedure
The selection of created chemisorbed precursor products depends on predefined occupation probabilities which remain constant during the simulated thin film growth. For example, considering a precursor yielding two possible chemisorbed precursor products P1 or P2 requires the corresponding occupation probabilities p P1 , p P2 ∈ [0, 1] , with This approach can be easily generalized to a set of N products {P1, P2, . . . , PN} with occupation probabilities p P1 , p P2 , . . . , p PN ∈ [0, 1] fulfilling the condition Materials 2019, 12, 2605 4 of 11 part is also different than in the original scheme of Hu et al. [15], where direct deposition of the product of the second half-cycle (i.e., -Al(OH)2) was performed. Steric hindrance and overlap restrictions at the surface corresponding to the real precursor deposition step were accounted for only in an approximate way, by defining an exclusion zone around each successfully deposited -Al(OH)2 group. In contrast, our approach employs a real, atomistic model of the precursor. After each attachment step, the structure was relaxed, either by local structure optimization or by short simulated annealing using molecular dynamics (MD). The energy of the chemisorption was calculated and used to calculate the acceptance criterion employing the usual Metropolis-Hastings rule [32]. The procedure was repeated until all enumerated surface -OH groups were either reacted or excluded from MC moves due to steric hindrance and overlap with chemisorbed precursor molecules. Alternatively, a predefined concentration of surface -OH groups can be left unreacted. The final precursor monolayer structure was relaxed and equilibrated using MD. The second ALD half-cycle was modeled by simply replacing all organic surface groups with -OH groups, followed by relaxation and equilibration using MD. This corresponds to full conversion of the products of the precursor pulse to surface hydroxyl species. Such an approach is particularly well suited for simulating ALD growth of systems with primary reactions that have relatively low activation energies [14]. It has been shown that, unlike thermal ALD, PEALD can be efficiently performed even at room temperature [6,7], justifying the assumption of low activation energies. It has also been shown by calculations that oxygen plasma can quantitatively oxidize organic surface precursor species [17]. The selection of created chemisorbed precursor products depends on predefined occupation probabilities which remain constant during the simulated thin film growth. For example, considering a precursor yielding two possible chemisorbed precursor products P1 or P2 requires the corresponding occupation probabilities { , This approach can be easily generalized to a set of products {P1, P2, …, PN} with occupation probabilities { , , … , ∈ [0,1]} fulfilling the condition + + ⋯ + = 1 (2)

Implementation
The implementation of our PEALD simulation protocol used the Python programming language. It was designed in a modular way employing the Atomic Simulation Environment (ASE) [33], which was used for storing and manipulating structure models as well as for the reading and writing of all necessary input and output files. Structure relaxation and simulated annealing steps were accomplished employing the General Utility Lattice Program (GULP) [34]. The Python library ASE provided a built-in GULP interface which enabled single-energy calculations, structure optimizations, and MD simulations.

Test Application-PEALD of SiO 2
As a test, the PEALD simulation protocol was applied to investigate the synthesis of SiO 2 thin films using a bis-diethylaminosilane (BDEAS) precursor with the structural formula H 2 Si(NEt 2 ) 2 , Et = C 2 H 5 . BDEAS belongs to the class of aminosilanes widely used in ALD, and it has been successfully applied for improving surface characteristics in PEALD of SiO 2 [10,35]. It has been shown in [17] that Si precursors with amino ligands can dramatically lower the activation energies for reactions of aminosilane with surface hydroxyl groups, making BDEAS-based PEALD a very good case for testing our simulation protocol.

Model Systems
Hydroxylated SiO 2 substrate and deposited precursor products are modeled using two types of cluster models of two neighboring -OH groups and a two-dimensional (2D) periodic surface. The implementation of our PEALD simulation protocol used the Python programming language. It was designed in a modular way employing the Atomic Simulation Environment (ASE) [33], which was used for storing and manipulating structure models as well as for the reading and writing of all necessary input and output files. Structure relaxation and simulated annealing steps were accomplished employing the General Utility Lattice Program (GULP) [34]. The Python library ASE provided a built-in GULP interface which enabled single-energy calculations, structure optimizations, and MD simulations.

Test Application-PEALD of SiO2
As a test, the PEALD simulation protocol was applied to investigate the synthesis of SiO2 thin films using a bis-diethylaminosilane (BDEAS) precursor with the structural formula H2Si(NEt2)2, Et = C2H5. BDEAS belongs to the class of aminosilanes widely used in ALD, and it has been successfully applied for improving surface characteristics in PEALD of SiO2 [10,35]. It has been shown in [17] that Si precursors with amino ligands can dramatically lower the activation energies for reactions of aminosilane with surface hydroxyl groups, making BDEAS-based PEALD a very good case for testing our simulation protocol.

Model Systems
Hydroxylated SiO2 substrate and deposited precursor products are modeled using two types of cluster models of two neighboring -OH groups and a two-dimensional (2D) periodic surface. Figure  3a,b shows both cluster models. The smaller, flexible cluster model has the composition Si2O7H6, while the bigger and more rigid cage cluster model has the composition Si8O14H8. The 2D periodic substrate model shown in Figure 3c

Methods
All DFT calculations employed our implementation of Kohn-Sham DFT for molecular and periodic systems [36][37][38][39][40] within the TURBOMOLE program package [41,42]. The Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [43], triple zeta valence plus polarization (def2-TZVP) [44] basis sets, and Grimme dispersion correction (DFT-D3) [45,46] were used. MD and MC simulations employed the ReaxFF reactive force field [27] using GULP-provided combined force field parameters [27,47,48]. Figure 4 shows the possible elementary surface reactions for PEALD synthesis of SiO2 thin films using a BDEAS precursor. The first half-cycle starts with the deposition of BDEAS on the hydroxylated SiO2 surface and its reaction with -OH groups. There are two possible reactions (1 and 2a in Figure 4) that yield the organic precursor bonded to one or two surface oxygen atoms, denoted

Methods
All DFT calculations employed our implementation of Kohn-Sham DFT for molecular and periodic systems [36][37][38][39][40] within the TURBOMOLE program package [41,42]. The Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [43], triple zeta valence plus polarization (def2-TZVP) [44] basis sets, and Grimme dispersion correction (DFT-D3) [45,46] were used. MD and MC simulations employed the ReaxFF reactive force field [27] using GULP-provided combined force field parameters [27,47,48]. Figure 4 shows the possible elementary surface reactions for PEALD synthesis of SiO 2 thin films using a BDEAS precursor. The first half-cycle starts with the deposition of BDEAS on the hydroxylated SiO 2 surface and its reaction with -OH groups. There are two possible reactions (1 and 2a in Figure 4) that yield the organic precursor bonded to one or two surface oxygen atoms, denoted as P1 and P2, respectively. Reaction 2b converts P1 to P2. After the plasma pulse, P1 and P2 lead to a surface Si atom with three or two -OH groups, denoted as S1 and S2, respectively. Alternatively, S1 can be converted to site S2 by reaction 3, which releases H 2 O. We have concentrated on reactions 1, 2a, 2b, and 3, since calculations indicate that oxidation of the organic precursor by oxygen radical species does not alter bonds between the precursor Si atom and the surface [17]. as P1 and P2, respectively. Reaction 2b converts P1 to P2. After the plasma pulse, P1 and P2 lead to a surface Si atom with three or two -OH groups, denoted as S1 and S2, respectively. Alternatively, S1 can be converted to site S2 by reaction 3, which releases H2O. We have concentrated on reactions 1, 2a, 2b, and 3, since calculations indicate that oxidation of the organic precursor by oxygen radical species does not alter bonds between the precursor Si atom and the surface [17].

Reaction Energies
As the first step, we investigated the energetics of the elementary reactions using DFT and the cluster models of two neighboring OH groups on a hydroxylated SiO2 surface, shown in Figure 3a,b. The aim of the cluster studies was to obtain insights into the intrinsic energetics of the reactions without the spatial and steric constraints present at silica surfaces. All three reactions (1, 2a, and 2b) of the precursor with the cluster model in Figure 3a are exothermic, with reaction energies ∆ r of −58.3, −64.3, −6.0 kJ/mol, respectively. Reaction 3, converting S1 to S2, is slightly endothermic, with an ∆ r of 37.5 kJ/mol. For the bigger cluster model in Figure 3b, the calculated reaction energies ∆ r of reactions 1, 2a, 2b, and 3 are -67.0, -30.5, +36.5, and +63.0 kJ/mol, respectively. In the next step, we investigated the energetics of the same reactions on a hydroxylated -quartz (0001) surface ( Figure  3c). The calculated reaction energies ∆ r were −63.0, −18.4, +44.6, and +60.9 kJ/mol for reactions 1, 2a, 2b, and 3, respectively. The value for exothermic reaction 1, involving only one surface OH group and yielding site P1, is very similar for all three models shown in Figure 3. Reaction energies for reaction 2a increase with increasing rigidity of the structure, while remaining exothermic. In contrast, reaction 2b is slightly exothermic for the small, flexible cluster model, but becomes endothermic for the structurally rigid models of the bigger cluster and surface. Reaction 3 is endothermic for all models, but reaction energies increase with increasing rigidity of the structure. Thus, reactions 1 and 2a are always energetically favorable. In contrast, reactions 2b and 3, transforming P1 to P2 and S1 to S2, respectively, are energetically unfavorable in a rigid atomic environment. This can be attributed to constraints of the surface on the one hand, since bringing two surface oxygen atoms close enough to create sites P2 and S2 is connected with energy penalty. On the other hand, this demonstrates the high stability of the Si-O bonds, especially in the case of S1.

PEALD Simulations
The cluster and periodic calculations demonstrate that the actual amount of precursor deposition products P1 and P2 as well as the plasma products S1 and S2 formed at each PEALD cycle will be influenced by the actual surface structure and the availability of neighboring -OH Elementary surface reactions during PEALD deposition of SiO 2 film using a bis-diethylaminosilane (BDEAS) precursor. P1 and P2 are products of the precursor deposition. S1 and S2 are products of the plasma pulse.

Reaction Energies
As the first step, we investigated the energetics of the elementary reactions using DFT and the cluster models of two neighboring OH groups on a hydroxylated SiO 2 surface, shown in Figure 3a,b. The aim of the cluster studies was to obtain insights into the intrinsic energetics of the reactions without the spatial and steric constraints present at silica surfaces. All three reactions (1, 2a, and 2b) of the precursor with the cluster model in Figure 3a are exothermic, with reaction energies ∆E r of −58.3, −64.3, −6.0 kJ/mol, respectively. Reaction 3, converting S1 to S2, is slightly endothermic, with an ∆E r of 37.5 kJ/mol. For the bigger cluster model in Figure 3b, the calculated reaction energies ∆E r of reactions 1, 2a, 2b, and 3 are -67.0, -30.5, +36.5, and +63.0 kJ/mol, respectively. In the next step, we investigated the energetics of the same reactions on a hydroxylated α-quartz (0001) surface (Figure 3c). The calculated reaction energies ∆E r were −63.0, −18.4, +44.6, and +60.9 kJ/mol for reactions 1, 2a, 2b, and 3, respectively. The value for exothermic reaction 1, involving only one surface OH group and yielding site P1, is very similar for all three models shown in Figure 3. Reaction energies for reaction 2a increase with increasing rigidity of the structure, while remaining exothermic. In contrast, reaction 2b is slightly exothermic for the small, flexible cluster model, but becomes endothermic for the structurally rigid models of the bigger cluster and surface. Reaction 3 is endothermic for all models, but reaction energies increase with increasing rigidity of the structure. Thus, reactions 1 and 2a are always energetically favorable. In contrast, reactions 2b and 3, transforming P1 to P2 and S1 to S2, respectively, are energetically unfavorable in a rigid atomic environment. This can be attributed to constraints of the surface on the one hand, since bringing two surface oxygen atoms close enough to create sites P2 and S2 is connected with energy penalty. On the other hand, this demonstrates the high stability of the Si-O bonds, especially in the case of S1.

PEALD Simulations
The cluster and periodic calculations demonstrate that the actual amount of precursor deposition products P1 and P2 as well as the plasma products S1 and S2 formed at each PEALD cycle will be influenced by the actual surface structure and the availability of neighboring -OH groups. In addition, temperature, pressure, and concentration of precursor species, and plasma energy are expected to play an important role. Accurate simulations including all these parameters would be a very complex task. As described earlier, our scheme takes the process parameters into account in a simplified way, employing predefined occupation probabilities at each deposition MC step. The probabilities for the creation of P1 (precursor bonded to one surface O atom) and P2 (precursor bonded to two surface O atoms) are p P1 and p P2 = 1 − p P1 , respectively. The value p P2 = 1 means that each Monte Carlo step attempts to create P2, and P1 is created only if P2 cannot be formed due to missing neighboring surface -OH groups or steric constraints. In contrast, p P2 = 0 means that no P2 is directly created during MC moves. As shown in Figure 5, these probabilities have significant impacts on the structures of PEALD layers. The simulated structures of SiO 2 deposited for probabilities p P2 = 0 (p P1 = 1), p P2 = 0.5 (p P1 = 0.5), and p P2 = 1 (p P1 = 0) clearly demonstrate the densification of the growing SiO 2 layer with increasing p P2 . For p P2 = 0, each precursor molecule is allowed to react with only one surface -OH group, yielding P1. The plasma pulse converts P1 to S1, with three hydroxyl groups connected to the deposited silicon atom. This increases the number of hydroxyl groups compared to the initial state, thus favoring the creation of voids, -OH impurities, and amorphous growth of deposited films. In contrast, for p P2 = 1, after the first ALD half-cycle, the maximum possible number of P2 precursor sites is created. The plasma pulse converts P2 to S2, with two hydroxyl groups connected to the deposited silicon atom. In this case, the number of hydroxyl groups remains unchanged after the deposition step compared to the initial surface state, and, in combination with compact Si-O-Si bridges, it results in a dense film growth. Figure 6 shows the average growth per cycle (GPC) and mass densities of deposited SiO 2 thin films obtained by 30 simulations of seven PEALD cycles for each of p P2 = 0.0, 0.25, 0.5, 1.0. Virtually identical results were obtained for MD simulations at temperatures ranging from 300 to 1000 K at a time scale of up to 10 ps. GPC and mass densities were basically converged after seven simulated PEALD cycles. Decreasing GPC and increasing density were observed for increasing p P2 , confirming the observed densification of thin films shown in Figure 5. The very low values for GPC and density in the case of p P2 = 0 are mainly due to the strong steric effects when only creation of singly coordinated P1 and S1 is allowed and creation of P2 and S2 is completely excluded. This confirms again the key role of doubly coordinated surface sites for the densification of SiO 2 PEALD thin films. groups. In addition, temperature, pressure, and concentration of precursor species, and plasma energy are expected to play an important role. Accurate simulations including all these parameters would be a very complex task. As described earlier, our scheme takes the process parameters into account in a simplified way, employing predefined occupation probabilities at each deposition MC step. The probabilities for the creation of P1 (precursor bonded to one surface O atom) and P2 (precursor bonded to two surface O atoms) are and = 1 − , respectively. The value = 1 means that each Monte Carlo step attempts to create P2, and P1 is created only if P2 cannot be formed due to missing neighboring surface -OH groups or steric constraints. In contrast, = 0 means that no P2 is directly created during MC moves. As shown in Figure 5, these probabilities have significant impacts on the structures of PEALD layers. The simulated structures of SiO2 deposited for probabilities = 0 ( = 1), = 0.5 ( = 0.5), and = 1 ( = 0) clearly demonstrate the densification of the growing SiO2 layer with increasing . For = 0, each precursor molecule is allowed to react with only one surface -OH group, yielding P1. The plasma pulse converts P1 to S1, with three hydroxyl groups connected to the deposited silicon atom. This increases the number of hydroxyl groups compared to the initial state, thus favoring the creation of voids, -OH impurities, and amorphous growth of deposited films. In contrast, for = 1, after the first ALD half-cycle, the maximum possible number of P2 precursor sites is created. The plasma pulse converts P2 to S2, with two hydroxyl groups connected to the deposited silicon atom. In this case, the number of hydroxyl groups remains unchanged after the deposition step compared to the initial surface state, and, in combination with compact Si-O-Si bridges, it results in a dense film growth. Figure 6 shows the average growth per cycle (GPC) and mass densities of deposited SiO2 thin films obtained by 30 simulations of seven PEALD cycles for each of = 0.0, 0.25, 0.5, 1.0. Virtually identical results were obtained for MD simulations at temperatures ranging from 300 to 1000 K at a time scale of up to 10 ps. GPC and mass densities were basically converged after seven simulated PEALD cycles. Decreasing GPC and increasing density were observed for increasing , confirming the observed densification of thin films shown in Figure 5. The very low values for GPC and density in the case of = 0 are mainly due to the strong steric effects when only creation of singly coordinated P1 and S1 is allowed and creation of P2 and S2 is completely excluded. This confirms again the key role of doubly coordinated surface sites for the densification of SiO2 PEALD thin films.   The actual reaction conditions were accounted for only in an approximate way in our occupation probability model. For instance, possible precursor desorption due to the substrate heating was not directly simulated. Possible influences from precursor temperature and concentration were not directly simulated, since our MC scheme deposits one precursor molecule after the other, disregarding the actual reaction path. However, we expect that our model can describe plasma energy effects, since the presence of plasma ions in the second half-cycle may affect reactions 2b and 3 (Figure 4), and therefore the occupation probabilities of doubly coordinated products P2 and S2.
Our simulation scheme has been tested for BDEAS as a bis-aminosilane precursor containing two amino ligands. For mono-aminosilane precursors, we expect simulated structures resembling the case of = 1 in Figure 5, but with denser films and less steric hindrance effects, since in this case there are no remaining precursor amino ligands after the first deposition half-cycle. Aminosilane precursors with three or four amino ligands can lead, in principle, to even higher surface densification, since they can increase the number of compact Si-O-Si bridges. However, the higher number of amino ligands favors steric hindrance effects, and triply and quadruply occupied surface sites strongly reduce the availability of surface -OH groups, thus inhibiting the overall film growth. In summary, the size and number of precursor amino ligands are expected to be determining factors for the steric effects and availability of -OH groups in our MC simulation scheme.

Conclusions
In this study, we presented an atomistic simulation scheme for PEALD that combines the Monte Carlo deposition algorithm and structure relaxation using molecular dynamics. In contrast to previous implementations, our approach employed a real, atomistic model of the precursor. This allowed us to account for steric hindrance and overlap restrictions at the surface corresponding to the real precursor deposition step. In addition, our scheme took various process parameters into account in a simplified way, employing predefined occupation probabilities for precursor products at each Monte Carlo deposition step. The new simulation protocol was applied to investigate the PEALD synthesis of SiO2 thin films using a bis-diethylaminosilane precursor. Initial analysis of the precursor products with hydroxylated SiO2 surfaces employing DFT calculations demonstrated a strong dependence of reaction energies on the surface structure, in particular for the precursor bonded to two surface oxygen atoms. PEALD simulations employing our new method revealed that precursor binding to one surface oxygen atom favors amorphous layer growth, a large number of -OH impurities, and the formation of voids. In contrast, precursor binding to two surface oxygen atoms leads to dense SiO2 film growth and reduction of -OH impurities. Increasing the probability for the formation of doubly bonded precursor sites is therefore the key factor for the formation of dense SiO2 PEALD thin films with reduced amounts of voids and -OH impurities. The actual reaction conditions were accounted for only in an approximate way in our occupation probability model. For instance, possible precursor desorption due to the substrate heating was not directly simulated. Possible influences from precursor temperature and concentration were not directly simulated, since our MC scheme deposits one precursor molecule after the other, disregarding the actual reaction path. However, we expect that our model can describe plasma energy effects, since the presence of plasma ions in the second half-cycle may affect reactions 2b and 3 (Figure 4), and therefore the occupation probabilities of doubly coordinated products P2 and S2.
Our simulation scheme has been tested for BDEAS as a bis-aminosilane precursor containing two amino ligands. For mono-aminosilane precursors, we expect simulated structures resembling the case of p P1 = 1 in Figure 5, but with denser films and less steric hindrance effects, since in this case there are no remaining precursor amino ligands after the first deposition half-cycle. Aminosilane precursors with three or four amino ligands can lead, in principle, to even higher surface densification, since they can increase the number of compact Si-O-Si bridges. However, the higher number of amino ligands favors steric hindrance effects, and triply and quadruply occupied surface sites strongly reduce the availability of surface -OH groups, thus inhibiting the overall film growth. In summary, the size and number of precursor amino ligands are expected to be determining factors for the steric effects and availability of -OH groups in our MC simulation scheme.

Conclusions
In this study, we presented an atomistic simulation scheme for PEALD that combines the Monte Carlo deposition algorithm and structure relaxation using molecular dynamics. In contrast to previous implementations, our approach employed a real, atomistic model of the precursor. This allowed us to account for steric hindrance and overlap restrictions at the surface corresponding to the real precursor deposition step. In addition, our scheme took various process parameters into account in a simplified way, employing predefined occupation probabilities for precursor products at each Monte Carlo deposition step. The new simulation protocol was applied to investigate the PEALD synthesis of SiO 2 thin films using a bis-diethylaminosilane precursor. Initial analysis of the precursor products with hydroxylated SiO 2 surfaces employing DFT calculations demonstrated a strong dependence of reaction energies on the surface structure, in particular for the precursor bonded to two surface oxygen atoms. PEALD simulations employing our new method revealed that precursor binding to one surface oxygen atom favors amorphous layer growth, a large number of -OH impurities, and the formation of voids. In contrast, precursor binding to two surface oxygen atoms leads to dense SiO 2 film growth and reduction of -OH impurities. Increasing the probability for the formation of doubly bonded precursor sites is therefore the key factor for the formation of dense SiO 2 PEALD thin films with reduced amounts of voids and -OH impurities. Funding: This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in the framework of the Priority Programme 1959 "Manipulation of matter controlled by electric and magnetic field: Towards novel synthesis and processing routes of inorganic materials" within the project SI 938/8-1.

Conflicts of Interest:
The authors declare no conflict of interest.