# Minimum Energy Atomic Deposition: A novel, efficient atomistic simulation method for thin film growth

Shivraj Karewar, Germain Clavier, Marc G.D. Geers, Olaf van der Sluis,
Johan P M Hoefnagels

► To cite this version:

Shivraj Karewar, Germain Clavier, Marc G.D. Geers, Olaf van der Sluis, Johan P M Hoefnagels. Minimum Energy Atomic Deposition: A novel, efficient atomistic simulation method for thin film growth. Surface and Coatings Technology, 2024, 494, pp.131462. 10.1016/j.surfcoat.2024.131462. hal-04859314

HAL Id: hal-04859314

https://hal.science/hal-04859314v1

Submitted on 30 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Full length article

# Minimum Energy Atomic Deposition: A novel, efficient atomistic simulation method for thin film growth

Shivraj Karewar<sup>a,b</sup>, Germain Clavier<sup>a,c</sup>, Marc G.D. Geers<sup>a</sup>, Olaf van der Sluis<sup>a</sup>, Johan P.M. Hoefnagels<sup>a,*</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

**Keywords:**

Thin film deposition

Atomistic simulations

Molecular dynamics

Residual stress

MEAD deposition

## ABSTRACT

Thin-film growth is an area of research concerned with complex phenomena happening at atomic scales. Therefore, molecular simulation has been an important tool to confront experimental results to theoretical assumptions. However, the traditional thin film growth simulation methods, i.e., Molecular Dynamics (MD) and kinetic Monte-Carlo (kMC) and combinations thereof, suffer from limitations inherent to their design, i.e., limitations in system size and simulation time for MD and predetermined reaction rates and reaction sites for kMC. Consequently, it is practically impossible to simulate the evolution of polycrystalline growth resulting in ~100 nm thick films with realistic stress fields and defect structures, such as grain boundaries, stacking faults, etc. In this work, we propose a versatile and efficient atomistic simulation method (Minimum Energy Atomic Deposition) which works by direct insertion of atoms at points of minimal potential energy through efficient scanning of candidate positions and rapid relaxation of the system. This method allows simulating ≥100 nm film thickness while mimicking experimental growth rates and high crystallinity and low-defect concentration and enables in-depth studies of atomic growth mechanisms, the evolution of crystal defects, and residual stress build-up. We demonstrate the efficiency and versatility of the method through the deposition of Al on Si, Al on Al, and Si on Si. The simulation results are systematically compared with experimental observations of thin-film deposition, yielding consistent observations. The method has been implemented in open-source LAMMPS software, making it easily accessible to the research community.

## 1. Introduction

Understanding the atomic mechanisms during thin-film growth and the resulting evolution of internal defects and intrinsic residual stresses constitutes a major scientific and industrially relevant research area, as it governs the performance, reliability, and durability of a wide range of microsystems (micro- and nano-electromechanical systems (MEMS, NEMS)) [1], battery electrodes, fuel cells, optoelectronic devices (LEDs), display technology (OLED), sensors and detectors (radiation detectors, biosensors), quantum dots, etc [2–5]. These applications rely on stacks of thin films with thicknesses ranging from tens to hundreds of nanometers, each with strict requirements in terms of high crystallinity and low concentrations of defects, such as grain boundaries (GBs), stacking faults (SFs), dislocation, voids, etc., therefore, much effort is put in controlling and reducing the (remaining) defect densities. Moreover, the built-up of intrinsic residual stresses in each layer of the thin film stack needs to be carefully controlled, to ensure mechanical

reliability or to enhance device functionality (e.g., micromachined ultra-sound transducers for in-body imaging that can be optimized by tuning the residual stress of the active layer [6]).

Intrinsic residual stresses result from the complex interplay between the atomic deposition processes (e.g., adsorption, diffusion, and desorption of reactive species coming to the surface), as well as from the interplay between defects in the bulk, all making the residual stress profile notoriously hard to predict. Therefore, revealing underlying mechanisms is key to optimizing the deposition process and improving the efficiency and performance of these micro-electronic devices. However, an atomistic simulation method that can simulate the growth of ~100 nm thick, high-crystallinity, low-defect films with realistic stress fields and defect structures is still being pursued. Such a method would enable in-depth studies of relevant atomistic growth mechanisms, the evolution of crystal defects, and residual stress build-up. In contrast, existing atomistic simulation methods have limitations that severely hamper such investigations.

<footnote><sup>a</sup> Eindhoven University of Technology, Dep. Mechanical Engineering, Eindhoven, The Netherlands</footnote>

<footnote><sup>b</sup> Advanced Materials Research Center, Technology Innovation Institute, Abu Dhabi, United Arab Emirates</footnote>

<footnote><sup>c</sup> CIMAP, Normandie Univ, ENSICAEN, UNICAEN, CEA, CNRS, CIMAP, 14000, Caen, Normandie, France</footnote>

<footnote>* Corresponding author.
E-mail address: J.P.M.Hoefnagels@tue.nl (J.P.M. Hoefnagels).</footnote>
In deterministic techniques, such as Molecular Dynamics (MD) sim-
ulations of thin film growth through direct atom sputtering, the in-
trinsic timescale (~fs) and length scale (~Å) are very short, and the
deposition rates are a factor of 10⁸ to 10¹⁰ higher than typical ex-
perimental growth rates, in order to achieve a relevant deposition
thickness in a reasonable computational time [7–23]. The short time
scale in MD simulations does not take into account the long time scale
phenomena and slower relaxation events, such as the adsorption of
radicals at the surface and the subsequent long-range surface diffusion
of these adatoms, with some exceptions [24]. The result is that most
adatoms never reach their low-energy equilibrium positions that lie
within the experimental adatom surface diffusion length, in contrast to
experimental conditions. Alternatively, MD simulations can be set up to
deposit high crystal-packing or crystalline morphologies [16,17], but
computational feasibility limits these simulations to ultra-short total
simulation times (up to a few hundred nanoseconds) and to maximally
tens of thousands of deposited atoms (for a computing time of sev-
eral days on a supercomputer). However, the defect structures and
stress-thickness profiles in the thin films used in the micro-electronic
applications, mentioned above, manifest themselves in the range of
tens to hundreds of nanometers film thickness. Therefore, simulation
of the evolution and interactions of the various defects and the built-
up of intrinsic residual stresses during deposition requires simulation
of high-crystallinity, low-defect growth of film volumes in the order
of ~10000 nm³ (e.g., 50 nm × 50 nm × 50 nm) or above, i.e., tens
of millions of deposited atoms, as this is the lower bound for many
defect types considering a minimal degree of statistics. This combina-
tion is inaccessible to MD simulations within realistic computational
resources.

In contrast, in On-lattice kinetic Monte Carlo (kMC) simulations of thin film growth, the atoms are directly placed at predetermined surface sites, while the evolution of these adatoms through the transition states is computed until their final deposition sites are reached, although the kinematics of these transition processes is not explicitly simulated but implicitly determined by predetermined reaction rates. As a result, kMC enables the deposition of micrometer-thick films, but it does require prior knowledge of all possible reaction rates and atomic deposition configurations [25–30]. Importantly, the atoms are not allowed to relax during deposition, resulting in artificially high atomic stresses. MC simulations still use symmetries relations to optimize simulation with regard to specific problems [31]. To integrate the advantages of MD and kMC, several combinations of MD + MC algorithms [32,33] or time-stamped force-bias Monte Carlo (tfMC) [33–36] were proposed. However, for all these methods, simulated film thicknesses are limited to tens of nanometers in thickness as well as in width, preventing a detailed investigation of the most relevant defects and defect mechanisms, such as the mutual interaction of dislocations, SFs, and/or grain/phase/twin boundaries, as well as their correlation with the origin of intrinsic residual stresses.

We propose a new method in which the atoms are directly deposited
(and relaxed) at the minimum energy positions (MEPs) of the surface
under deposition, completely surpassing the computationally expensive
evolution of the system through the transition states, hence the name
Minimum Energy Atomic Deposition (MEAD). This approach depends
on the atoms’ local equilibrium energies and is independent of the
adatom kinetics. The MEAD method is motivated by the fact that thin
films with high crystallinity and low-defect concentrations at near-
equilibrium growth conditions and high Young’s modulus and yield
strength are preferred for the range of applications mentioned above,
whether serving electrical, optical, and/or mechanical functionality.
Note, however, that these high-quality films still have many types
of defects and a complex residual stress profile. Often the industry
aims to (further) reduce this defect density/distribution and control
the residual stress profile [1,37]. MEAD simulations can aid in this
process by simulating under which conditions the grain shapes/defect
densities/distributions/residual stress profiles are changing and in what

manner exactly. In the industrial processes, this is achieved by slow
growth rates at high deposition temperatures to maximize the surface
diffusion length of adatoms, to let (almost) all of them find the MEPs
and deposit there. To this end, the MEAD algorithm parameters can
be tuned to simulate thin film growth where the grain shapes/defect
densities/distributions/residual stress profiles (approximately) agree
with the experiments, after which the parameters can be further tuned
to optimize the film quality, providing guidelines for experimental-
ists to improve their films. Finally, note that the existing simulation
methods – MD and off-lattice kMC – already well capture the far-from-
equilibrium growth conditions, but they struggle to accurately simulate
near-equilibrium deposition conditions, especially for simulation do-
mains that are large enough to study the interactions of different kinds
of defects and the built-up of the residual stress over the full film
thickness as used in the application, requiring volumes of ~10000 nm³
or more as mentioned above. Therefore, to complement the existing
methods, the MEAD method is focused on simulating growth conditions
that result in high-crystallinity thin films with low defect concentration
and/or large film thickness, with the aim to further optimize these
growth conditions.

As mentioned, the core idea behind the MEAD method is to detect
and place atoms directly at the MEPs. The thus simulated thin-film
growth mimics experimentally observed near-equilibrium deposition
rate, while high computational efficiency is achieved by fast and ac-
curate determination of the MEPs of the surface under deposition, as
will be demonstrated below. The aim of this work is not to replace
the existing algorithms, such as MD or kMC, but instead to improve
upon the existing methods in certain aspects by integrating them into
the MEAD method. One can readily combine the MEAD algorithm with
the existing techniques such as minimization, tfMC, and MD equilibra-
tion, as detailed later in the manuscript, thereby also incorporating
the capabilities of these techniques. This can be done by tuning the
parameters of the MEAD algorithm to grow structures that are not
known in advance under various growth conditions to yield valuable
insight into the growth mechanism(s) that can support the experiments.

The manuscript is organized as follows. First, the MEAD method,
with its basic conceptual and computational steps, is introduced (Sec-
tion 2). Subsequently, some of the many possibilities of the MEAD
method are successfully showcased for the industrially relevant cases
of Al and Si atoms on a Si substrate (Section 3). Finally, the results
and perspectives offered by the MEAD method are discussed, as well
as current limitations and routes to overcome them in future work
(Section 4).

## 2. The MEAD method

In MEAD, the atoms are directly deposited at the MEPs. To this end,
it is important to efficiently compute the potential energy surface (PES)
of the deposition surface, as experienced by isolated atoms arriving at
the surface, assuming their kinetic energy is low enough to not disturb
the surface structure in a significant way. We will explore two ways
to compute the PES: (i) an accurate but too-slow method to compute
the full PES, which will serve as a reference for the MEAD method and
is presented first, and (ii) the fast, approximate method used in MEAD
that is accompanied by an additional energy minimization step to find
the MEPs, as explained thereafter.

### 2.1. *Accurate, but slow computation of the PES*

The PES can be computed accurately by introducing a ghost atom
above the surface that moves downward at a very small step size (here
taken 1 × 10⁻² nm) passing through the surface, while the position with
the lowest potential energy is stored. This process is repeated for a
regular grid in the XY plane with a very small lateral step size (here
again taken 1 × 10⁻² nm) [38]. Subsequently, a continuous PE surface is
constructed from the selected low-energy positions using the Gaussian
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** (a–d) The potential energy surfaces (PES) for an Al atom on a (001) Al surface for four different potentials (ADP, EAM, Morse, and Tersoff), color-coded by the color bars. The distance between two MEPs is indicated for ADP potential. (e–f) The potential energy (PE) as a function of the Z coordinate of the Al atoms on the PES for different potentials, the topmost atom on the Al substrate surface is at the height of z = 20.24 Å. (g–h) The morphologies of the Al-on-Al thin films for Morse and Tersoff potentials, respectively, color-coded as per a-CNA (green—FCC, blue—BCC, red—HCP, gray—unidentified), the inset shows the defect structures observed in the deposited systems.</figcaption></figure>

Density Surface algorithm for surface reconstruction [39]. The accuracy of the PES can be further increased by decreasing the step size.

Fig. 1(a–d) shows such a PES for an Al atom on top of a (001) Face Centered Cubic (FCC) Al crystal surface, for four interatomic potentials, namely Angular Dependent Potential (ADP) [40], Embedded Atom Method (EAM) [8], Morse [41], and Tersoff [42]. The PES surface shows the MEPs in a dark blue color, corresponding to an energy of ~ -2.5 eV, where new atoms should end up after a sufficiently long diffusion time. In this simple case, the MEPs are periodically located on the surface of this crystalline substrate. As the atomic forcefields are simulated through continuous and derivable functions, it is assumed that the potential energy smoothly decreases towards the MEPs in their surroundings. This was confirmed via computation of PES for the four potentials, which all show a smooth PES (see Fig. 1), and was also found for all cases tested in this work (Al and Si as the depositing atoms on Al and Si substrates).

It is useful to study the details of the PES near their minima and compare the differences for different interatomic potentials, as this explains why some of them are more suitable for thin-film deposition simulations using the MEAD method than others. To this end, Fig. 1(e,f) shows the distribution of the potential energy as a function of the Z coordinates above the deposition surface, for the four interatomic potentials. For the ADP, EAM, and MORSE potentials, the energy level of locations near the MEPs is highly correlated with their height. In

contrast, the environment of minimum energy positions for the Tersoff potential are very flat, i.e., there exist regions of minimum energy, see Fig. 1(d) and 1(f). These plateaus of equal energy, instead of a clearly identifiable single minimum energy location, make it difficult to identify the optimum minimum energy position for deposition. The effect of the interatomic potential on the deposition simulations will be discussed by running MEAD on all four interatomic potentials (see Section 4).

2.2. *Fast, approximate PES computation in MEAD*

One could compute a high-resolution PES, similar to Fig. 1(a), deposit an atom at a MEP, and repeat this process. However, this is (i) computationally very expensive and (ii) of insufficient accuracy because even a very small step size of the ghost atom of 1 × 10⁻² nm typically still results in too large steps in energy (see Fig. 1) resulting in high stresses and because relaxation of the neighboring atoms can move the precise location of the MEP away from the computed one. Therefore, it will always be needed to add a system energy minimization step (or thermal equilibration step) to move the atoms that are close to their computed MEP precisely to the real (relaxed) MEP locations. But with such an added energy minimization step it is unnecessary to compute the PES with high accuracy, i.e., fast computation of an approximate PES suffices, as will be shown below. Therefore, for the MEAD method, we propose the following 3-step iterative approach:
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Lower-left: A full Al(001) surface with PES as experienced by a single Al atom. The PES is shown here as the (underlying) continuous surface of varying color according to the energy scale bar given at the center of the figure. The color bar at the center-right corresponds to the magnitude of the displacement arrows in the top-row figures. (a, b, c) Three viewpoints of a regular periodic grid of phantom atoms (shown in pink) are used to probe the PES above the substrate surface. Phantom atoms at different locations in space converge to the same position after energy minimization, as shown by the arrows pointing to the blue central atom, located at the same position as the minimum in the underlying PES as indicated by the green vertical arrow. The transparent phantom atoms, only shown in (a), indicate positions that converge to other minimum energy positions. (d) The grid of phantom atoms, colored to their potential energy, after they have been slightly displaced in random directions, to prevent probing the exact same locations at each deposition step. (e) Only the atoms with the lowest energies are kept in the system, and these move to their MEPs upon energy minimization, as shown by the gray arrows.</figcaption></figure>

* Step (1): fast, approximate PES computation by simultaneously prob-
ing the unrelaxed energies around the deposition surface for a
massive number of *non-interacting* phantom atoms.

* Step (2): selection of the phantom atoms with the lowest energies, i.e. those that are closest to a MEP.

* Step (3): turning them into actual atoms that interact with the surface and moving them to their real MEPs (using energy minimization or thermal equilibration).

The method is illustrated for Al deposition on an Al(001) substrate (Al-on-Al) in Fig. 2 and for Al deposition on a Si(001) substrate (Al-on-Si) in Fig. 3.

**Step (1):** The 3D geometry of the deposition surface is first constructed, in low resolution, using the alpha-probe method [43], and a periodic grid of non-interacting phantom atoms with a grid spacing *a* are placed in a volume surrounding the deposition surface. The height of the grid ranges from -0.5 Å to +2.5 Å relative to the local Z coordinate of the deposition surface, and it spans the whole XY plane (see Figs. 2(a, b, c) and 3(a)). These phantom atoms only interact with the substrate and previously deposited atoms, but not with each other. This is required since we are using models that take into account three-body interactions, and as such, the interaction of phantom atoms with one another and with the surface would not only shift the potential energy but change its value with a dependence on the grid geometry without any physical meaning. The individual interaction energy *E*ᵢ of each phantom atom *i* is computed directly, yielding a potential energy value at all phantom atom positions simultaneously, thus identifying all proximate MEPs locations in one molecular statics computational step. High accuracy is not needed, as atoms within a fairly large distance to a MEP still move to the MEP during energy minimization, as demonstrated by the arrows in Fig. 2(a, b, c). Finally, a random displacement of maximum magnitude *a*/2 in the X and Y directions and *a*/3 along the Z-axis is added to each phantom atom during each deposition loop (see Figs. 2(d) and 3(b)) to avoid probing the same positions between successive deposition loops.

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Illustration of the steps of the MEAD algorithm for the Al-on-Si system: (a) A periodic grid of non-interacting Al phantom atoms (pink) above the Si(001) substrate surface (brown atoms). (b) A random displacement is applied to each phantom atom in the grid to prevent probing the exact same locations at each deposition step. The potential energy of the phantom atoms is calculated at this step, (c) Removal of the high-energy phantom atoms, keeping the phantom atoms with a PE below a certain threshold, as shown by their blue color. The cluster of atoms in the black circles is closer than what would be physically possible. Only the atom with the lowest PE in the cluster is kept in the system, the rest is discarded. (d) The system's energy is then minimized, resulting in displaced atoms with a darker blue color, while the orange arrows show the displacements during minimization.</figcaption></figure>
**Step (2):** The phantom atoms with the lowest potential energies are kept in the simulation system, while the others are discarded (Figs. 2(e) and 3(c)). The selection of these low-energy phantom atoms is made here through a deterministic criterion, $|\frac{E'-E_{Min}}{E_{Min}}| < \lambda$, where $E_{Min}$ is the lowest of all phantom atom energies computed in the deposition loop, while $\lambda$ is an energy cutoff parameter that controls the allowed deviation from $E_{Min}$. Since the phantom atoms do not mutually interact, the selected atoms can be non-physically close to each other, as shown by the black circles in Fig. 3(c). Hence, for atoms closer than a cutoff radius $r_c$, only the one with the lowest potential energy is kept. $r_c$ is chosen just higher than the first neighbor distance of the crystal structure of the pure material.

Note that a more advanced selection criterion can easily be selected. For instance, to simulate low-temperature non-equilibrium deposition resulting in low-density films, a non-deterministic selection criterion can be designed that captures the stochastic nature of the deposition process, including the constraint that each depositing atom/radical can only reach MEPs that lie within the experimental adatom surface diffusion length from its random impact location. However, this is outside of the scope of the present article and shall be detailed in future work.

**Step (3):** For the remaining isolated phantom atoms, the interatomic interaction is enabled. Subsequently, either the system's energy is minimized, as demonstrated first in Sections 2 and 3, or a thermal equilibration step is performed, as demonstrated in Section 4. During the energy minimization or thermal equilibration step, computed here with respectively the FIRE algorithm [44] and time-stamped force-bias Monte Carlo (tfMC) [33–36], the box dimensions are kept fixed. Figs. 2(e) and 3(d) show the result of an energy minimization step, in which the selected atoms move to their closest MEPs. The surface atoms are also able to reorganize, however, since the surface is already close to the energy minimum and the newly inserted atoms are close to the MEPs, the reorganization for all atoms is typically minor.

This three-step deposition loop is repeated until the desired deposition thickness is achieved. Depending on the choice of the MEAD parameters ($\lambda$ and $a$) and the size of the simulation system, it takes approximately 10 to 15 MEAD iterations to deposit one monolayer of thickness, i.e., about ~7 to ~10% fraction of monolayer is deposited during each loop. One deposition loop takes approximately 120 seconds on 1 node with 16 processors, depending on the simulation system size and interatomic potential. The computational time scales roughly linearly with the number of atoms due to the use of standard MD algorithms and the use of a chargeless model.

### 3. Proof of principle: Large-scale deposition simulations

The MEAD method is first demonstrated on the homoepitaxial growth of Al (Al-on-Al), for which a clean crystalline structure is expected. Next, Si growth on a silicon wafer (Si-on-Si), is studied, which is interesting as it does not always result in homoepitaxial growth due to the Si substrate surface reconstruction. Then, we present the industrially relevant case of Al deposition on the same silicon wafer (Al-on-Si), which is studied in terms of achievable film thicknesses, defect structures, and associated residual stress built up. Finally, an improvement to the MEAD method is introduced to take temperature effects into account. Film growth evolution movies of the three deposition cases are given in the supplementary information.

The MEAD algorithm only has two user-defined parameters of significant influence, i.e., the grid spacing and the energy cut-off $\lambda$, which were chosen as, respectively, $a = 0.07$ nm and $\lambda = 0.15$, based on a study of their effect that is discussed in the next Section 4.1. Additionally, the choice of the interatomic potential – here the ADP potential [40] is consistently employed based on a comparison in the next Section 4.2 – and the choice between energy minimization and thermal equilibration in step (3) are important, as will be shown below. Simulations are carried out using the LAMMPS code [45,46] and analyzed using OVITO tools adaptive common neighbor analysis (a-CNA) and identify diamond structure (IDS) [47–49].

Regarding the substrates, for all simulations and the supplementary information videos, the following simulation parameters were used: The 24.2×24.2×2.42 nm³ Al(001) substrate (of 93 600 atoms) and 24.5×24.5 × 3.29 nm³ Si(001) substrate (of 101 250 atoms) were equilibrated at 298 K and 673 K, respectively, to obtain a stress-free configuration using an isothermal-isostress (Nσxx,yy,T) ensemble with Nosé-Hoover integrator thermostat and anisotropic barostats at 0 MPa. Coupling constants were set to 100 fs and 1000 fs with a timestep value $dt = 1$ fs. Their energy was minimized, keeping the box dimensions fixed. Thin film growth occurs in the Z-direction where a large vacuum was left in the simulation box so the surfaces can freely relax their stress along the Z-axis, and with periodic boundary conditions in the XY plane. The FIRE algorithm [44] was employed for energy minimization of the substrates with relative energy variation of $10^{-8}$ and force tolerance of $10^{-8}$ eV Å⁻¹.

#### 3.1. Al-on-Al: homoepitaxial growth

Fig. 4(a) shows the layer-by-layer atomic deposition mechanism for Al deposition on Al substrate. As can be seen, a layer-by-layer deposition mechanism is clearly identifiable. Such behavior is coherent with experimentally observed deposition from CVD or Molecular Beam Epitaxy (MBE) homoepitaxial growth that happens with limited chemical reactions (see, e.g., Ref. [50]). This shows that the MEAD method can simulate the experimentally observed deposition phenomena and also keep the crystalline deposition morphology. The formation of defects is also observed (see bottom insets of this figure) such as SFs, deposition of individual adatoms, terrace plane formation, atomic cluster formation, presence of vacancies, and kink formation and propagation. This illustrates the capability of MEAD to simulate the relevant defect mechanisms, which are not explicitly included in its algorithm.

#### 3.2. Si-on-Si: deposition of $3.87 \times 10^6$ Si atoms

Now, let us focus on the Si-on-Si deposition example. In this case, the same “ADP” potential and set of parameters were used as that for the Al-on-Al case, but homoepitaxial growth was not observed. This is not surprising when considering the following factors. First, the Si(001) substrate surface is not crystalline but consists of dimers at the surface, see the inset of Fig. 4(b). As such, the ADP potential correctly reproduces experimental behavior, as Si atoms tend to evolve as dimers on the surface [51]. The consequence is that the depositing atoms experience MEPs that do not correspond to crystalline structure sites, contrary to the Al-on-Al case. Second, this Si surface reconstruction is coupled to a high energy barrier for the atoms to leave their coupled states [52]. Therefore, the pure Si homoepitaxial deposition is carried out at a much higher temperature [51,53–58]. However, it can be observed (using IDS analysis on the top inset of Fig. 4(b)) that small Si crystal nuclei are being formed in the bulk, with 1st and 2nd crystal shells formed around a substantial number of atoms. This shows that the energy minimization step is effective in bringing the atoms to their local crystalline configuration minima in the bulk. But, starting from the initial amorphous layer (induced by the silicon substrate dimers), the relative position of the local minima is typically too far removed from the ideal crystalline configuration such that the amorphous configuration is passed onto the next layer, a process that is repeated for each layer. Even the small silicon crystallites that do form are not sufficient to pull the system into large-scale crystalline growth in subsequent layers, not even after more than 30 nm film thickness, see Fig. 4(b).

The Si crystalline nuclei in the system indicate the possibility of crystalline growth if the conditions allow, i.e., at a higher temperature. Indeed, epitaxial growth of Si-on-Si by MBE occurs at very slow deposition rates and at a high temperature above 750 K in order
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Large-scale deposition of test cases with zoomed views of various phenomena observed. (a) An example of Al deposition on a (001) Al substrate system showing zoomed views during the deposition process in the bottom insets. The atoms are color-coded as per a-CNA where green denotes FCC, red denotes HCP, gray denotes either unidentified atoms or GB atoms or atoms at the surface, and blue denotes BCC. (b) An amorphous deposition morphology of $3.87 \times 10^6$ Si atoms on Si substrate. The two insets on top of this figure show magnified views of the Si(001)(2x1) surface reconstruction with dimers and crystalline nuclei formation within the bulk of the deposited film. (c,d) Film growth at 1200 K results in a crystalline deposition morphology for (c) a Si-on-Si(100) substrate and (d) a Si-on-Si(111) substrate, despite the Si surface reconstructions. (b–d) The atoms are color-coded as per IDS analysis where dark blue denotes cubic diamond, light blue denotes cubic diamond (1st or 2nd nearest neighbor), and yellow indicates atoms having a hexagonal diamond structure locally. (a–d) The atoms below the dashed lines constitute the substrate.</figcaption></figure>

to provide enough energy for full atomic reconfiguration. We have
tested this for Si-on-Si(100) and Si-on-Si(111) substrate surfaces by
combining MEAD with tfMC at a much higher temperature (1200 K)
and with an increased number of tfMC iterations (5000 instead of
1000) to (partially) account for the almost infinite atomic equilibration
times at the very slow MBE deposition rates, and a slightly reduced
energy cutoff parameter value ($\lambda$ = 0.10 instead of 0.15) to (partially)
account for the thermally-induced high adatom surface diffusion length
in MBE resulting in almost all adatoms finding a very low energy
deposition site. Under these conditions that approach MBE conditions,

indeed, crystalline growth occurs for Si-on-Si(100) and Si-on-Si(111);
see Figs. 4(c–d). Note that the gray atoms at the surface are caused
by the continuous formation of ‘dimers’, disrupting the first few lay-
ers, yet, these atoms converge to their bulk crystalline positions in
a few MEAD deposition cycles. Interestingly, no defects are detected
at the substrate-bulk interface or in the crystalline bulk, therefore
these growth conditions seem to correspond to the high end of the
temperature range of MBE conditions, i.e., around 1100 K, used in
the industry to achieve high-quality crystalline growth with almost
perfect epitaxial layer quality and very low defect densities [59]. The
remaining difference in temperature (between 1200 K and 1100 K) may
solely be caused by the accuracy of the Si interatomic potential at high
temperatures. The structural characteristics of the substrate and thin
film at 1200 K for the Si(100) and Si(111) surfaces are shown in the
supplementary information Fig. S1, including views of, respectively,
the (2x1) and (7x7) surface reconstruction, the epitaxial substrate-thin
film interfaces, and the surface-roughness-induced amorphous layers at
the surface of the depositing thin films. These results show that the
MEAD algorithm is able to capture the transition from amorphous to
crystalline Si growth, in agreement with experiments, demonstrating
the ability to simulate various experimental conditions.

### 3.3. Al-on-Si: deposition of 7.78 × 10⁶ Al atoms

The main proof-of-principle simulation case of Al-on-Si helps in understanding the role of the surface in the structure formation. An example of such an Al deposition with large thickness (213 nm) and a high number of deposited atoms (7.78 × 10⁶ atoms) is shown in Fig. 5(a), while a second example with 60 nm thickness is shown in 5(b). As can be seen, the deposition structure transitions first from amorphous growth at the interface to nano-polycrystalline growth with many planar defects, such as SFs, twin boundaries (TBs), and GBs, and afterwards to single crystalline growth without defects above ~100 nm. In this case, these defects do not originate from an incorrect selection of the MEAD algorithm parameters, but they are generated because of different reasons, including the dimers at the initial Si substrate surface, the lattice mismatch between the Si substrate and Al FCC structure, and the simultaneous formation of nano-grains with different orientations necessarily resulting in GBs. The fact that the defects are not artificially generated by the MEAD method is evidenced by the fact that, for the same MEAD algorithm parameters and same Al interatomic potential, (i) the Al-on-Al system in Fig. 4(a) shows single crystal epitaxial growth, while (ii) for the Al-on-Si system in Fig. 5(a) the grain size quickly increases with increasing film thickness, i.e., the well-known process of competitive grain growth, eventually converging to single crystal growth (due to the limited lateral dimensions of the deposition system). The convergence to epitaxial crystalline deposition clearly shows the ability of the MEAD method to deposit thin films with high crystallinity and low-defect concentration at near-equilibrium growth conditions.

Nevertheless, the interface with the substrate still shows a large number of defects and many atoms in the amorphous phase before convergence. Therefore, the displacement of the deposited atoms before and after minimization is investigated in more detail. Interestingly, the initially deposited atoms do end up near computed MEPs of the PES. This can be seen in Fig. 6(a), which shows the displacement magnitude corresponding to the unrelaxed and relaxed energies of Al atoms during the first deposition loop. It can be seen that the relationship between the energy and the displacement during relaxation is approximately linear. The average atomic displacement during minimization is only 0.3 Å. In subsequent loops, as the depositing surface becomes irregular, this linear relationship is lost. This can be seen in Fig. 6(b), which shows the same quantities as Fig. 6(a) but after deposition of 40 nm of film. Larger displacements of up to 1.5 Å are observed because of the presence of planar defects within the deposited Al thin film at this stage (see defects such as SFs, TBs, and GBs in Fig. 5). This difference can be seen directly by comparing Fig. 6(c) and 6(d), which show the evolution of the system before and after the minimization for the two cases. Still, the typical small atom displacements and fast convergence (at a few hundred iterations) demonstrate that the identified locations of approximate MEPs are close to the actual MEPs on the PES, resulting in limited reorganization during the energy minimization.

Experimentally, deposition of high-crystallinity, low-defect films is
typically carried out at high temperatures to allow for adatoms to
diffuse on the deposition surface. So far, in our proof-of-principle sim-
ulations, only configurations of minimal energy have been taken into

<figure><img src="image_8.png" /><figcaption>Fig. 5. (a) A 213 nm-thick film of 7.78 × 10⁶ Al atoms deposited on a Si(001) wafer. (b) Another 60 nm-thick film to better show the defect structures and topography. (c) The system shown in (b) but after an additional thermal equilibrium at 400 K. Various marked defects (d) and the corresponding Von Mises stress field (e) for the cross-section marked in (c). In (a–d) atoms are identified by a-CNA as: FCC (green), HCP (red) for SFs or TBs, or unidentified structures (gray) for GBs, non-crystalline regions, or surface atoms.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. The unrelaxed and relaxed energies, for Al deposition on Si substrate, plotted against the relative displacement before and after minimization, respectively, for atoms deposited during the (a) first deposition loop and (b) after 40 nm deposition. Black arrows show examples of energy reduction magnitudes during minimization. Zoom-in at the atoms, for Al deposition on Si substrate, before (pink atoms) and after (PE color scale bar) the energy minimization during the (c) first deposition loop and (d) after 40 nm deposition. The arrows show the displacements, scaled times 3 for better visualization, and color-coded with the displacement magnitude.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** The evolution of the (a) number of deposited Al atoms and (b) fraction of defect atoms as a function of MEAD steps on a (001) Al substrate surface for several $\lambda$ values at $a = 0.7$ Å, respectively, (c) the visualization of the deposition morphology for the same parameters, (d–e) the evolution of deposited Al atoms as a function of MEAD steps and the fraction of non-crystalline defect atoms as a function of an increase in deposited atoms, respectively, for different values of $a$ parameter (expressed in Å units) at a fixed $\lambda = 0.15$, and (f) the corresponding visualization of the deposition morphology.</figcaption></figure>

account without any effect of the temperature. To investigate the effect of temperature, we have carried out post-process finite temperature MD simulations of the deposited system to investigate how thermal treatments, often used in industry, affect the structure, especially at the interface. To this end, 1 ns-long thermal equilibration simulations were carried out using Nosé-Hoover thermostat with timestep $dt = 1$ fs and coupling constant $\tau_T = 100 \cdot$ dt. An initial random velocity distribution was used with average kinetic energy so that the instantaneous temperature at the start is 400 K and the momentum is canceled in all directions. The effect of such treatment is shown in the difference between Fig. 5(b) and 5(c). As can be seen, after the additional thermal equilibrium, the interface contains a lower density of defects, while the atoms are reorganized into bigger crystal grains in the whole system. A second 1 ns simulation with similar parameters was carried out, with the system configuration recorded every 1 ps along with the stress contribution of each atom. The average Voronoi volume, $V$, was calculated for each atom, and the Von Mises effective stress value [60] was computed based on a stress tensor that is computed from the contribution of the individual atoms using the Voronoi/atom and stress/atom algorithms of LAMMPS. It is not straightforward to quantitatively compare this highly inhomogeneous stress field to experimentally measured global stress values, but it does correlate well with the observed defect structure. Fig. 5(d) and (e) nicely illustrate how grain boundaries and defects contribute to the residual stress of the grown film. A more quantitative analysis of the correlation is left out for future work, but it can be noted that most of the contribution of the stress components to the Von Mises effective stress comes from the normal stress in the X- and Y-plane, since the surface was not allowed to relax the pressure during thermal treatment (similar to real-world conditions). This aspect is in agreement with general theoretical assumptions [1].

# 4. Discussion

## 4.1. Effects of the MEAD algorithm parameters

The MEAD algorithm has two main parameters that need to be optimized, i.e., the grid spacing $a$ of the phantom atoms in the regular periodic grid (before the random displacements are applied) and the energy cutoff $\lambda$ for selecting the lowest energy phantom atoms. These parameters strongly influence the computational efficiency. Smaller values of $a$ provide a better probing of the volume around MEPS, however, it increases the phantom atom densities and thus the potential energy computing time, especially for multi-body interatomic potentials. A larger value of $\lambda$ will select a larger number of phantom atoms, some of which will be located at a larger distance from their nearby MEP, increasing the computation time of the energy minimization step. Therefore, we first study the influence of $\lambda$ and $a$ (for the ADP interatomic potential) to establish their optimized values.

The effect of parameters $\lambda$ and $a$ were benchmarked for Al deposition on the (001) Al substrate surface (4.82 × 4.82 × 2.42 nm³), shown in Fig. 7. As expected, a larger number of atoms is deposited in 300 deposition iterations for a higher value of $\lambda$ (e.g., 0.2) compared to a lower value (e.g., 0.05), at fixed $a = 0.7$ Å, see Fig. 7(a). It can be seen that $\lambda = 0.15$ is an optimum value as it allows depositing a larger number of Al atoms compared to $\lambda$ values of 0.05 and 0.10 (see Fig. 7(a)) without generating additional defects. Although the number of deposited atoms is even larger for $\lambda = 0.20$, the crystal structure analysis by the a-CNA method reveals that the fraction of defect atoms is also higher (Fig. 7(b)), while the surface roughness profile also starts deviating significantly for $\lambda > 0.15$, see Fig. 7(c). In contrast, $\lambda \le 0.15$ results in layer-by-layer deposition of a homoepitaxial crystalline structure, as one would expect. As a general rule, for other material systems, $\lambda$ should be less than or equal to a threshold value for which
deposition morphology does not generate artificial defects (here found to be 0.15) and, at the same time, results in efficient computation.

Regarding grid spacing *a*, a lower value will create a dense grid of phantom atoms that can locate the MEPs more precisely. However, this comes at the expense of a higher computational cost. To benchmark the *a* parameter, we used different values at a fixed $\lambda = 0.15$, shown in Fig. 7(d–f). For a higher grid spacing of *a* = 0.1 nm, after 300 deposition steps, the number of deposited atoms is significantly lower, and the fraction of unknown atoms (defects with non-crystalline structure) is larger as compared to the smaller grid spacing of *a* = 0.05 nm. This indicates that the deposition with a larger grid spacing of 0.1 nm is not able to find all of the MEPs efficiently and accurately. We also indicated the number of loops taken to deposit the 12 000 atoms, the maximum reached for *a* = 0.1 nm, by a horizontal dashed line in Fig. 7(d). When using 16 processors, the time taken to reach this value was 23, 14, and 24 min for *a* value of 0.05 nm, 0.07 nm, 0.1 nm, respectively. Therefore, *a* = 0.07 nm indicates an optimum magnitude as it simulates the same deposition morphology as *a* = 0.05 nm but with only 61% of the computation time, whereas it is even faster than *a* = 0.1 nm as this simulation misses too many MEPs, which not only makes it less accurate but also slower than *a* = 0.07 nm. As a general rule, for other material systems, we advise the grid spacing to be smaller than the distance between the neighboring MEPs on the PES, which should establish the identification of all MEPs on the depositing surface. At the same time, *a* should be sufficiently large to speed up the computation, as discussed above, without generating artificial defects.

## 4.2. Effect of the interatomic potential

It was observed in Fig. 1(a,b) that the PES morphology for the EAM potential is similar to that of the ADP potential. Since the potentials are very similar in form, this is not a surprising result. The Morse potential also shows a similar deposition morphology, with a linear correlation between energy and height of the surface points. However, when tested with the *a* and $\lambda$ parameters defined in the previous section, it generates some BCC phase during deposition, see Fig. 1(g). This happens because Morse potential's [41] energy difference between the BCC and FCC phases of Al is unphysically small (0.02 eV/atom), whereas *ab-initio* calculations show this difference to be 0.11 eV/atom [61]. In contrast, the other interatomic potentials, used in this work, are much closer to the *ab-initio* value with 0.06 eV/atom for EAM [8], 0.07 eV/atom for ADP [40], and 0.1 eV/atom for Tersoff [42]. Since the transition energy is much lower for the Morse potential, it is relatively easier to nucleate and grow the BCC phase, as can be seen in Fig. 1(g).

Finally, the Tersoff potential shows a significant amount of amorphous structure, see Fig. 1(h), which appears because of the minimum energy plateaus mentioned previously, which also results in a non-linear relation between height and potential energy on the surface. Note that the difference between maxima and minima potential energy (see Fig. 1(e–f)) for ADP, EAM, Morse, and Tersoff potentials is 0.96 eV, 0.6 eV, 0.26 eV, 0.21 eV, respectively. The MEAD algorithm can be run with any interatomic potential, however, naturally, the quality of the resulting deposited film is dependent on the model used to describe interatomic interactions. Therefore, the generation of the BCC structures for the Morse potential and the amorphous phase for the Tersoff potential is not an artifact of the MEAD algorithm or the MEAD parameters but exposes the inherent limitations of the use of interatomic potentials.

## 4.3. Inclusion of the temperature: MEAD coupled with tfMC

The previous examples used energy minimization to reorganize the atoms being deposited and find their final minimum energy configuration, while a post-process thermal equilibration was used to introduce thermal effects and provide thermal energy for atomic reorganizations, as was demonstrated in Fig. 5(b,c). Alternatively, to take into account

<figure><img src="image_11.png" /><figcaption>Fig. 8. Comparison of the MEAD deposited morphology obtained using (a) FIRE energy minimization and (b) tfMC equilibration at 300 K. The atomic color coding is the same as that in Fig. 5.</figcaption></figure>

the effect of temperature on the deposition process, the MEAD method can be directly coupled with time-stamped force-bias Monte Carlo (tfMC) equilibration, instead of FIRE energy minimization, after every deposition step. The tfMC algorithm allows recrystallization at a much longer timescale compared to post-process MD integration [35]. This, in turn, helps to simulate longer timescale atomic reorganization that is hardly accessible to traditional MD simulations.

A comparison of the MEAD deposition system using FIRE energy minimization and tfMC equilibration at 300 K is shown in Fig. 8. tfMC simulations were performed with maximum displacement parameter $\Delta = 1.5 \times 10^{-2}$ nm, which is just above the recommended value of $1 \times 10^{-2}$ nm from Mees et al. [34]. A quantitative comparison of the structural evolution during deposition is given in Fig. 9. It is interesting to see that tfMC equilibration results in an increase in a crystalline structure way earlier in the deposition compared to FIRE minimization and an overall notable decrease of the amorphous structures. The number of nano-crystalline grains and density of defects were, however, comparable. Then again, the structure obtained using tfMC is stable during the following thermal NVT simulations at 300 K and higher temperature, in contrast to the one obtained from the FIRE-generated structure.

In both cases, after thermal treatment, Al FCC crystals show a nearly clean interface with the Si wafer, with an amorphous interphase layer of only a few atoms thick, before growing predominantly along the [001] growth direction (compare Fig. 5(c) with Fig. 8(b)). This interphase layer is caused by the above-mentioned surface reconstruction of the equilibrated Si substrate containing dimers with bridge-like structures and valleys where Al atoms tended to get stuck (see Fig. 4(b)). This complicated surface strongly affects the formation of initial crystal nuclei and leads to an initial amorphous deposition morphology. These observations may seem to contradict the general consensus in experimental literature, which lacks sub-nanometer resolution, that Al-on-Si growth at room temperature or above occurs through immediate polycrystalline growth, however, they are in good agreement with detailed X-ray scattering studies that did detect that, for Al(111) and (001) growth on Si(001), the Al-Si interface is hard to stabilize resulting in an ultra-thin amorphous interphase layer [62,63]. As expected, for metallic substrates that do not undergo surface reconstruction, the MEAD simulations produce homoepitaxial growth. Particularly, for Al-on-Al growth, homoepitaxial FCC growth is observed in Fig. 4, while Fe-on-Fe and Ti-on-Ti growth results in, respectively, homoepitaxial BCC and HCP growth (not shown). In the case of Al-on-Al, the use of tfMC did not result in a noticeable change in the obtained structure compared to FIRE energy minimization, as the latter already produced the expected layer-by-layer homoepitaxial growth.
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Left: Comparison of the structure obtained during deposition of Al atoms on a Si(100) surface through the use of FIRE (black) and tfMC at 300 K (gray). The atoms are classified either in FCC crystal, HCP (at SFs or TBs), or unknown/amorphous structures (UNK). Values are normalized by the total number of atoms in the system at each step, including the initial $2 \times 10^4$ atoms of silicon in a diamond configuration (UNK). Compared to FIRE-based structures, tfMC-based structures are more ordered, containing much less unknown/amorphous atoms and defects. Right: Comparison of the structures during subsequent 300 K NVT thermalization. tfMC-based structures are stable in time, while the structure obtained through FIRE minimization keeps on reorganizing to some extent, especially at GBs.</figcaption></figure>

## 4.4. Implementation details and computational costs

The MEAD algorithm is coupled with other atomistic simulation elements, specifically the interatomic potentials and the energy minimization or thermal equilibration algorithm, making it a versatile simulator of atomistic thin-film deposition. For these atomistic simulation elements, in-built computation functions inside LAMMPS are used, making the method easily accessible to the community. The scripts are available upon request and will be made publicly available on repositories.

The methodology can be easily applied to simulate a large variety of material systems using readily available interatomic potentials. We have provided insights and guidelines for the selection of appropriate interatomic potentials for MEAD simulations in Section 4.2, where correct identification of the MEPs is demonstrated irrespective of the magnitude of the energy difference between the maxima and minima on the PES of a given interatomic potential. The MEAD method also provides freedom in the relaxation algorithm, either by energy minimization, thermal equilibration, or another scheme, while warranting the key elemental atomic mechanism for deposition of high-crystalline, low-defect film growth, i.e., that adatoms are able to reach the MEPs. This is achieved by directly depositing the atoms close to the MEPs, while the relaxation algorithm employed in MEAD only provides minor adjustments (of less than 1 Å for most atoms, as shown in Fig. 6(b)). For the choice of this relaxation step, we have compared conjugate gradient (CG), FIRE, and room-temperature tfMC algorithms, finding very similar deposition morphologies, but differences in computational efficiency, iterations, and operating temperature. This is non-trivial. For instance, in an interesting paper by Alvarez et al. [64], the authors managed to reproduce the very different deposition morphologies of, e.g., Au versus W, by introducing the so-called 'most coordinated relaxation' and 'random relaxation' schemes, which follow a different set of rules for atom movement; however, as stated by the authors, these algorithmic schemes are without fundamental backing. In MEAD, the underlying physics of atomic deposition is well justified.

Regarding computational costs, for the 213 nm thick Al–Si system of Fig. 5(a), 7.78 million atoms were deposited in 528 CPU hours using 1 node with 16 processors. With multiple parallel nodes or the use of GPUs instead of CPUs, and/or a faster interatomic potential (e.g., embedded atom method), the MEAD algorithm would enable the deposition of ~50 million in accessible computational times. Note that

our initial MEAD simulations with the tfMC algorithm (Fig. 8(b)) show that tfMC-based MEAD is even faster than FIRE-based MEAD. For the tfMC-based MEAD simulations, 53 nm deposition height was simulated in 43.5 h on a single node with 64 processors. In contrast, for the same resources, the fastest MD sputtering simulations deposit only up to ~1×10⁴ atoms in 1 h and at extremely high growth rates (>10⁸ nm/ps), potentially leading to artificially high residual stresses and affecting the deposited film morphology, while off-lattice kMC simulations deposit only a few atoms per hour. On-lattice kMC simulations can deposit a huge number of atoms but under the severe limitations of pre-defined lattice sites, à priori known reaction rates, and absent stresses.

## 4.5. Possibilities and limitations of MEAD

We have shown that the MEAD method captures the relevant atomic processes for different material systems, with realistic simulations of the microstructure evolution and defect mechanisms during thin-film growth, yielding high film crystallinity and low defect density, and that the simulated atomic systems can be analyzed further to study many interesting effects and mechanisms. For instance, the typically strong influence of the substrate is of clear interest; indeed, as seen in Figs. 4(a) and 5, the dominant growth direction for Al changes from <100> growth on an Al(001) substrate to <111> growth on a Si(001) substrate, in agreement with experiments [65,66]. There is also a strong incentive to better understand the effect of heat treatments, as it is one of the preferred ways of the industry to optimize the microstructure and residual stress, which can be studied by equilibration at high temperatures, as shown in Fig. 5. A largely open field of research is the build-up of intrinsic residual stresses during thin-film deposition, which is important for many applications and can be compared to experimental stress-thickness profiles [1].

Moreover, several extensions to the MEAD algorithm, some already mentioned above, can easily be implemented to expand the range of atomic and microstructure mechanisms that can be studied, including:

* **stochasticity of atomic deposition:** Instead of the so-far used hard energy cut-off parameter ($\lambda$) to determine which phantom atoms are deposited in the system, Maxwell-Boltzmann statistics can be easily implemented in the selection procedure to capture the statistical distribution of atoms over the various energy states at thermal equilibrium at a certain temperature. This will help to randomize the morphology of the deposited film corresponding to a certain temperature.
* finite adatom surface diffusion: During the growth process, the adatom surface diffusion length is the key factor determining the surface roughness evolution and resulting defect density and porosity, yet this important parameter can only be experimentally measured indirectly from the surface roughness evolution [67]. By making the above-mentioned selection criterion for the MEPs dependent on the adatom surface diffusion length (i.e., adatoms can only deposit within a certain radius from their random landing site), its influence on the surface roughness evolution can directly be studied. This is a relatively simple pre-processing step before the energy selection procedure (step (2)).

* directional/oblique deposition: So far, all MEPs, also the ones that are hidden behind a roughness peak, have been included as potential atom deposition sites, from which the actual atom deposition sites are then selected by the energy selection procedure. In directional/oblique deposition, the direction of the incoming deposition flux is controlled or manipulated to achieve desired film characteristics. This can be implemented in MEAD, also as a pre-processing step before the energy selection procedure, by only including the MEPs, which are in direct line of sight from an imagined inclined flux of incoming atoms to the list of potential atom deposition sites. Note that, in reality, the extent to which film growth will occur in the shadow regions behind the roughness peaks will depend strongly on the adatom surface diffusion length, therefore both pre-processing steps would wisely be included simultaneously.

* long-range bulk diffusion: Similar to the incorporation of energy minimization techniques (e.g., conjugate-gradient or FIRE) and thermal equilibration (tfMC) in the MEAD algorithm, the algorithm of Grand Canonical Monte Carlo (GCMC), which is readily available in LAMMPS, can be incorporated to simulate atom surface diffusion and diffusion of interstitial atoms in the bulk [68], while semi-GCMC can be incorporated to simulate the diffusion of substitutional atoms in the bulk [69,70]. This would naturally be implemented as step (4) of the MEAD method, opening a range of bulk atomic mechanisms in high-crystallinity, low-defect, 'thick' films, including solid solution mechanisms such as dislocation pinning and the formation of nano-precipitates in alloys [71].

With these demonstrated and future possibilities, the MEAD methodology shows clear potential to unlock the way to observe, analyze, extract, and calculate a wide range of useful material properties, characteristics, and atomic mechanisms during and after film deposition. This list of interesting phenomena includes (but is not limited to):

* the initial stages of nucleation and growth and the type of growth mechanism (e.g., island growth versus layer-by-layer growth),

* the overall resulting microstructure (amorphous/poly-crystalline/ single-crystal/multi-phase) and their transition temperature and mechanisms, as shown for Si-on-Si deposition in Fig. 4,

* the evolution of defects such as vacancies, voids, dislocations, SFs, GBs, TBs, and phase boundaries, etc.

* the density profile over the thickness, which is important for, e.g., optoelectronic and mechanical properties,

* the evolution of the surface roughness with film thickness, which can easily be compared to simple Atomic Force Microscopy experiments,

* the build-up of intrinsic residual stresses during film growth and the resulting residual stress versus thickness profile, see Fig. 5, which is particularly important for MEMS and NEMS,

* the 'recovery' and 'recrystallization' temperature of the deposited thin film, relevant for, e.g., multi-layer microdevices,

* the development of a thin film 'texture', i.e., the grain orientation distribution, important for, e.g., the elasto-plastic film properties,

* the evolution of specific grain shapes, resulting in, e.g., columnar growth, important for, e.g., the wear resistance of coatings,

* the dependence of a wide range of material properties on the deposition parameters, such as the temperature,

* and phase equilibria, e.g., coarsening of multi-phase microstructures and defect evolution phenomena in alloys, including dislocation pinning by solutes and nano-precipitates.

To study these phenomena and to provide guidelines to experimentalists, the MEAD simulation parameters can be tuned to match the final properties of the simulated thin films to experimental observations. These tuning parameters include, for instance:

* the temperature during deposition via the integrated tfMC protocol, to match the experimental deposition temperature

* the other tfMC parameters, such as the number of iterations, to control the degree of atomic relaxation that changes most strongly with the deposition rate,

* the temperature and other parameters of a post-process thermal equilibration routine, to match the evolution of the microstructure, defect densities, and residual stress profile to those measured in annealed specimens,

* and the energy cutoff parameter ($\lambda$), which controls the deposition morphology, from amorphous to crystalline (see Section 4.1), as well as the surface roughness evolution and defect density and porosity [67].

More parameters should become available for fine-tuning in future versions of the MEAD algorithm, though their testing is beyond the scope of this work. These parameters include, for instance:

* the surface temperature used in a Maxwell-Boltzmann statistics energy selection criterion, which enables studying the thin film defect density and its temperature dependence, which can be tuned to match the experimental degree of local heating by, e.g., a plasma-induced ion bombardment, separate of the bulk relaxation temperature (used in tfMC).

* the adatom surface diffusion length, which can be tuned to replicate the experimental surface roughness evolution, defect density and porosity,

* the angle and spread of a directional deposition extension, which can be used to study the effects and potential benefits of oblique deposition, e.g., by controlling the microstructure grain morphology towards maximum wear resistance,

* the temperature and number of iterations in a GCMC or semi-GCMC step, which controls the extent of bulk diffusion mechanisms, such as defect evolution and dislocation interaction phenomena, and phase equilibration and nanoprecipitation in alloys.

Of course, the MEAD method also has clear limitations. For now, it has only been applied to single-element type depositions, although deposition of multi-element systems, such as the relevant Al-Cu binary alloy system, are currently successfully being explored as examples of two-phase materials. As the name of the method implies, *Minimum Energy Atomic Deposition* only supports the growth from dissociating molecules where the solid reaction product at the surface consists of one or more atoms that each deposit at their low energy locations during near-equilibrium growth, while by-products are removed as gases from the reaction chamber. This means that the MEAD method cannot simulate film growth by direct deposition of molecules or molecule fragments into a film, or other cases where the molecules do not undergo initial dissociation at the surface into separate adatoms. This excludes a wide range of (exotic) deposition processes, however, it does include many of the industrially relevant applications of high-crystallinity, low-density films for microdevices, deposited with one of the common thin film deposition methods (CVD, PVD, PECVD, MBE, PLD, etc.)

Another limitation of MEAD is that the effect of the kinetic energy of the incoming atoms is not directly considered in the algorithm. If the surface undergoes kinematic bombardment by high-energy ions or
atoms, resulting in consecutive collisions with the previously deposited
bulk atoms, then this effect cannot be simulated by MEAD, and this
simulation belongs to the realm of MD sputter deposition simulation. If,
however, the arriving atoms stay at the surface but bring kinetic energy
to the top surface, then this surface energy increases the adatom surface
diffusion length, the effect of which can be simulated by MEAD indi-
rectly by changing the adatom surface diffusion length. By design, the
MEAD algorithm focuses on the final deposition state, not the interme-
diate atomic mechanisms. Therefore, if the interest lies in investigating
the transition of the atoms between adsorption sites, then the MEAD
algorithm would have to be coupled with the existing methods such as
MD, off-lattice kMC, and/or GCMC. Moreover, the MEAD method only
simulates the resulting film, and not the gas fluxes of volatile species
or resputtered atoms emitted from the surface. As a consequence,
parameters such as the deposition efficiency cannot be calculated. Nev-
ertheless, the effect of elusive deposition quantities, such as the sticking
coefficient of the depositing species or the resputtering rate of the
deposited film, can still be studied with MEAD, though only indirectly,
as can be understood as follows. A high sticking coefficient reduces the
adatom surface diffusion length as highly sticking species get trapped in
local energy minima, whereas sputtering preferentially removes surface
atoms at the higher energy states, thereby effectively increasing the
adatom surface diffusion length. This opens an indirect route to study
the effect of the average sticking coefficient or the resputtering rate on,
e.g., the surface roughness evolution or void density, by varying the
adatom surface diffusion length of the MEAD simulations (discussed
above). The effect of other deposition quantities may be studied in a
similar indirect fashion. Finally, note that, as MEAD avoids the slow
surface diffusion and reaction barrier computations, adding tempera-
ture effects through isothermal equilibration, tfMC, or both requires
careful balancing between accuracy and computational efficiency. It
should thus be checked that the obtained deposition morphology is
sufficiently insensitive to the MEAD deposition parameters.

5. Conclusion

In conclusion, this paper presents an extremely simple and compu-
tationally efficient atomistic simulation method to simulate thin-film
deposition. The method aims to overcome the limitations of existing
methods to enable in-depth studies of atomic growth mechanisms, the
evolution of crystal defects, and residual stress build-up during thin-
film deposition—an important and active research area. The method
allows depositions of tens of millions of atoms at near-equilibrium
deposition rates in reasonable computational time while allowing for
the growth of >100 nm-thick films with high crystallinity and low
defect concentration. This has been validated on three material systems,
demonstrating that the method is much faster than MD simulations.
However, the proposed method does not aim to replace the existing
algorithms such as MD or kMC, but instead, improve upon the existing
methods in certain aspects by integrating them into the MEAD method.
One can combine the MEAD algorithm with the existing techniques
such as energy minimization with FIRE, thermal equilibration with
tfMC, and MD equilibration, thereby also incorporating the capabili-
ties of these techniques. Unlike on-lattice kMC, the method does not
need pre-determined reaction rates and deposition sites—both of which
naturally emerge from the simulation methodology. Our results are
also coherent with the experimental literature on thin film deposition.
Hence, this method can be broadly used to simulate and analyze
the evolution of defect structures and the build-up of residual stress
during thin film growth and subsequent heat treatments. The use of
in-built functions in LAMMPS makes the method easily accessible to
the community.

CRediT authorship contribution statement

**Shivraj Karewar:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Germain Clavier:** Writing – review & editing, Visualization, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Marc G.D. Geers:** Writing – review & editing, Supervision, Resources, Investigation. **Olaf van der Sluis:** Writing – review & editing, Supervision, Project administration, Investigation, Funding acquisition. **Johan P.M. Hoefnagels:** Writing – review & editing, Validation, Supervision, Methodology, Investigation, Formal analysis, Conceptualization.

Declaration of competing interest

The authors declare that they have no known competing finan-
cial interests or personal relationships that could have appeared to
influence the work reported in this paper.

Acknowledgments

This work was carried out on the Dutch national e-infrastructure
with the support of SURF Cooperative. The authors would like to ac-
knowledge the Moore4Medical project, which received funding within
the Electronic Components and Systems for European Leadership Joint
Undertaking (ECSEL JU) in collaboration with the European Union's
H2020 Framework Programme (H2020/2014-2020) and National Au-
thorities, under grant agreement H2020-ECSEL-2019-IA-876190. The
authors acknowledge the fruitful discussions with Prof. Markus Hütter
at TU/e.

**Appendix A. Supplementary data**

Supplementary material related to this article can be found online
at https://doi.org/10.1016/j.surfcoat.2024.131462.

Data availability

Data will be made available on request.

References

[1] G. Abadias, E. Chason, J. Keckes, M. Sebastiani, G.B. Thompson, E. Barthel, G.L. Doll, C.E. Murray, C.H. Stoessel, L. Martinu, J. Vac. Sci. Tech. A 36 (2018) 020801.

[2] M. Tilli, T. Motooka, V.-M. Airaksinen, S. Franssila, M. Paulasto-Kröckel, V. Lindroos (Eds.), Micro and Nano Technologies, in: Handbook of Silicon Based MEMS Materials and Technologies, William Andrew Publishing, Boston, 2015, pp. 3–85.

[3] Z. Qi, H. Wang, Research (2020) https://spj.science.org/doi/pdf/10.34133/2020/2969510.

[4] D. Yang, first ed., Handbook of Photovoltaic Silicon, Springer, 2019.

[5] T.F. Kuech (Ed.), second ed., Handbook of Crystal Growth, North-Holland, Boston, 2015.

[6] J. Wang, Z. Zheng, J. Chan, J. Yeow, Microsyst. Nanoeng. 6 (2020) 73.

[7] H. Zhang, Z.N. Xia, Nucl. Instrum. Methods Phys. Res. B 160 (2000) 372.

[8] X.W. Zhou, H. Wadley, R. Johnson, D. Larson, N. Tabat, A. Cerezo, A. Petford-Long, G. Smith, P. Clifton, R. Martens, T. Kelly, Acta Mater. 49 (2001) 4005.

[9] X. Guo, P. Brault, Surf. Sci. 488 (2001) 133.

[10] C. Chu, T. Chen, Surf. Coat. Technol. 201 (2006) 1796.

[11] M. Taguchi, S. Hamaguchi, Thin Solid Films 515 (2007) 4879, the Third International Symposium on Dry Process (DPS 2005).

[12] K. Inai, Y. Kikuhara, K. Ohya, Surf. Coat. Technol. 202 (2008) 5374.

[13] J. Houska, S. Mraz, J.M. Schneider, J. Appl. Phys. 112 (2012) 073527.

[14] G. Kokkoris, P. Brault, A. Thomann, A. Caillard, D. Samelor, A.G. Boudouvis, C. Vahlas, Thin Solid Films 536 (2013) 115.

[15] L. Xie, P. Brault, A. Thomann, J. Bauchire, Appl. Surf. Sci. 285 (2013) 810.

[16] J. Houska, Surf. Coat. Technol. 254 (2014) 131.

[17] J. Houska, Surf. Coat. Technol. 304 (2016) 23.
[18] J. Houska, P. Machanova, M. Zitek, P. Zeman, J. Alloys Compd. 828 (2020) 154433.

[19] P. Brault, Front. Phys. 6 (2018) http://dx.doi.org/10.3389/fphy.2018.00059.

[20] I.V. Chepkasov, V.S. Baidyshev, E.V. Sukhanova, M.A. Visotin, P. Süle, Z.I. Popov, Appl. Surf. Sci. 527 (2020) 146736.

[21] F.V. Grigoriev, V.B. Sulimov, A.V. Tikhonravov, Comput. Mater. Sci. 188 (2021) 110202.

[22] P. Brault, A. Thomann, M. Cavarroc, Eur. Phys. J. D 77 (2023) 19.

[23] M. Ghaemi, A. Lopez-Cazalilla, K. Sarakinos, G. Rosaz, C. Carlos, S. Leith, S. Calatroni, M. Himmerlich, F. Djurabekova, Surf. Coat. Technol. 476 (2024) 130199.

[24] M. Castro, R. Cuerno, M. Nicoli, L. Vazquez, J.G. Buijnsters, New J. Phys. 14 (2012) 103039, For those deposition cases where there is no adatom surface diffusion (i.e., at sufficiently low temperatures) regular MD simulations only need to simulate the adatom to surface trajectory and the thermal spike cooling which speeds up the simulation orders of magnitude, however, without surface diffusion so-called 'cauliflower growth' occurs resulting in high roughness and high porosity films.

[25] C.C. Battaile, D.J. Srolovitz, J.E. Butler, J. Appl. Phys. 82 (1997) 6293.

[26] G.H. Gilmer, H. Huang, T. Diaz de la Rubia, J. Dalla Torre, F. Baumann, Thin Solid Films 365 (2000) 189.

[27] C.C. Battaile, D.J. Srolovitz, Annu. Rev. Mater. Res. 32 (2002) 297.

[28] F. Nita, C. Mastail, G. Abadias, Phys. Rev. B 93 (2016) 064107.

[29] E. Chason, A.F. Bower, J. Appl. Phys. 125 (2019) 115304.

[30] N. Cheimarios, D. To, G. Kokkoris, G. Memos, A.G. Boudouvis, Front. Phys. 9 (2021) 631918.

[31] P. Wang, W. He, G. Mauer, R. Mücke, R. Vaßen, Surf. Coat. Technol. 335 (2018) 188.

[32] M. Taguchi, S. Hamaguchi, J. Appl. Phys. 100 (2006) 123305.

[33] E. Neyts, A. Bogaerts, Theor. Chem. Acc. 132 (2013) 1320.

[34] M.J. Mees, G. Pourtois, E.C. Neyts, B.J. Thijsse, A. Stesmans, Phys. Rev. B 85 (2012) 134301.

[35] K.M. Bal, E.C. Neyts, J. Chem. Phys. 141 (2014) 204104.

[36] R. Namakian, B.R. Novak, X. Zhang, W.J. Meng, D. Moldovan, Appl. Surf. Sci. 570 (2021) 151013.

[37] L.B. Freund, S. Suresh, Thin Film Materials: Stress, Defect Formation and Surface Evolution, Cambridge University Press, 2004.

[38] K. Seino, A. Oshiyama, Appl. Surf. Sci. 561 (2021) 149927.

[39] M. Krone, J. Stone, T. Ertl, K. Schulten, in: M. Meyer, T. Weinkaufs (Eds.), EuroVis - Short Papers, The Eurographics Association, 2012.

[40] S. Starikov, I. Gordeev, Y. Lysogorskiy, L. Kolotova, S. Makarov, Comput. Mater. Sci. 184 (2020) 109891.

[41] L.A. Girifalco, V.G. Weizer, Phys. Rev. 114 (1959) 687.

[42] G. Plummer, G.J. Tucker, Phys. Rev. B 100 (2019) 214114.

[43] A. Stukowski, Jom 66 (2014) 399.

[44] J. Guénolé, W.G. Nöhring, A. Vaid, F. Houllé, Z. Xie, A. Prakash, E. Bitzek, Comput. Mater. Sci. 175 (2020) 109584.

[45] S. Plimpton, J. Comput. Phys. 117 (1995) 1.

[46] A.P. Thompson, H.M. Aktulga, R. Berger, D.S. Bolintineanu, W.M. Brown, P.S. Crozier, P.J. in 't Veld, A. Kohlmeyer, S.G. Moore, T.D. Nguyen, R. Shan, M.J. Stevens, J. Tranchida, C. Trott, S.J. Plimpton, Comp. Phys. Comm. 271 (2022) 108171.

[47] A. Stukowski, MSMSE 18 (2010) 015012.

[48] A. Stukowski, MSMSE 20 (2012) 045021.

[49] E. Maras, O. Trushin, A. Stukowski, T. Ala-Nissila, H. Jónsson, Comput. Phys. Comm. 205 (2016) 13.

[50] G. Springholz, A.Y. Ueta, N. Frank, G. Bauer, Appl. Phys. Lett. 69 (1996) 2822.

[51] R.A. Wolkow, Phys. Rev. Lett. 74 (1995) 4448.

[52] T. Yamasaki, T. Uda, K. Terakura, Phys. Rev. Lett. 76 (1996) 2949.

[53] T. Hsieh, K. Jung, D. Kwong, S. Lee, J. Electrochem. Soc. 138 (1991) 1188.

[54] A. Ichimiya, T. Hashizume, K. Ishiyama, K. Motai, T. Sakurai, Ultramicroscopy 42 (1992) 910.

[55] A.H. Mahan, J. Carapella, B.P. Nelson, R.S. Crandall, I. Balberg, J. Appl. Phys. 69 (1991) 6728.

[56] S. Yu, E. Gulari, J. Kanicki, Appl. Phys. Lett. 68 (1996) 2681.

[57] A.-P. Prskalo, S. Schmauder, C. Ziebert, J. Ye, S. Ulrich, Comput. Mater. Sci. 50 (2011) 1320.

[58] P. Reinig, F. Fenske, B. Selle, W. Bohne, J. Röhrich, I. Sieber, W. Fuhs, Appl. Surf. Sci. 227 (2004) 114.

[59] E. Kasper, Appl. Phys. A 28 (1982) 129.

[60] R. Hill, The Mathematical Theory of Plasticity, vol. 11, Oxford University Press, 1998.

[61] Y. Mishin, D. Farkas, M.J. Mehl, D.A. Papaconstantopoulos, Phys. Rev. B 59 (1999) 3393.

[62] F.J. Lamelas, M.-T. Tang, K. Evans-Lutterodt, P.H. Fuoss, W.L. Brown, Phys. Rev. B 46 (1992) 15570.

[63] S. Nishikawa, K. Tani, T. Yamaji, J. Mater. Res. 7 (1992) 345.

[64] R. Alvarez, A. Garcia-Valenzuela, G. Regodon, F.J. Ferrer, V. Rico, J.M. Garcia-Martin, A.R. Gonzalez-Elipe, A. Palmero, Nanotechnology 35 (2023) 095705.

[65] C. Bisch, E. Boellaard, G.C.A.M. Janssen, P.F.A. Alkemade, S. Radelaar, Thin Solid Films 336 (1998) 84.

[66] N. Jackson, Vacuum 132 (2016) 47.

[67] J.P.M. Hoefnagels, A novel diagnostic approach for studying silicon thin film growth, (Phd thesis 1 (research tu/e / graduation tu/e)), Appl. Phys. and Science Education, 2005.

[68] S.S. Moeini-Ardakani, S.M. Taheri-Mousavi, J. Li, Modelling Simul. Mater. Sci. Eng. 29 (2021) 055018.

[69] Y. Mishin, Acta Mater. 52 (2004) 1451.

[70] K. Ueno, Y. Shibuta, IOP Conf. Ser.: Mater. Sci. Eng. 861 (2020) 012064.

[71] B. Sadigh, P. Erhart, A. Stukowski, A. Caro, E. Martinez, L. Zepeda-Ruiz, Phys. Rev. B 85 (2012) 184203.