Monte Carlo and Kinetic Monte Carlo Models for Deposition Processes: A Review of Recent Works

Monte Carlo (MC) and kinetic Monte Carlo (kMC) models are widely used for studying the physicochemical surface phenomena encountered in most deposition processes. This spans from physical and chemical vapor deposition to atomic layer and electrochemical deposition. MC and kMC, in comparison to popular molecular methods, such as Molecular Mechanics/Dynamics, have the ability to address much larger time and spatial scales. They also offer a far more detailed approach of the surface processes than continuum-type models, such as the reaction-diffusion models. This work presents a review of the modern applications of MC/kMC models employed in deposition processes.


INTRODUCTION
Monte Carlo (MC) and kinetic Monte Carlo (kMC) are widely used methods in many fields of science and engineering: From materials science and polymers properties [1], astrophysics and black holes mergers [2] to computational geometry and volume approximation [3]. Their popularity in materials science stems from their inherit ability to simulate the molecular level of materials seamlessly. In MC/kMC, the particles (molecules, atoms, beads) move stochastically according to specific rules (events/processes), transferring the system randomly over the phase-space and approximating the mean values of various properties. In contrast to other molecular methods such as molecular dynamics (MD), the system in MC cannot easily be trapped in local energy minima and even if it is trapped, it can be "kicked out" to other states by incorporating sophisticated events. Furthermore, kMC filters out vibrational movements, allowing it to run over much larger spatial and time scales than MD. Especially in film growth, snapshots of MC/kMC simulation can be directly related and compared to scanning tunneling microscopy images.
MC/kMC have notable applications in the study of film deposition processes-probably the most important in the fabrication of semiconductor devices. Because of their importance, deposition processes have extensively studied via MC/kMC models either standalone or combined with other models in the context of multiscale modeling [4]. Nevertheless, there is no current review of the works that use MC or kMC concerning film deposition. This work is a critical review of the works on deposition processes utilizing MC/kMC models published within the last 5 years. Following the introduction, the deposition processes are shortly presented and the works that use MC/kMC models are reviewed. The paper ends with summary and outlook.

THE MC AND KMC METHODS FOR DEPOSITION PROCESSES
In MC and kMC algorithms, sequential events are performed stochastically. MC solves the steady state Master Equation (ME) and kMC the transient one. The transient ME reads, where p j(i) is the probability of the system to be found in state j(i) at time t. T ij and T ji denote the transition rate or transition probability from state i to j and vice versa.
Each event occurs at a certain probability/rate to form a Markov chain [5]. To generate the Markov chain the desired probability distribution p i(j) must obey the detailed balanced condition, where E i(j) is the energy of the system in state i(j). Metropolis et al. [6] proposed that, so that the system will unconditionally move from state i to j if Practically, a random number ξ is chosen between (0, 1] and if ξ < exp − E k b T , the system moves to the state j, otherwise the move is rejected. In this way, different states of the system are generated and the thermodynamic average of a quantity q i reads kMC method solves Equation (1). The most popular algorithm proposed by Bortz et al. [7] is termed as the N-fold method. In the N-fold method random transitions from i to j are performed unconditionally based on the transition rates, so that more likely transitions are selected more often. Every transition event i is assigned a rate which reads, where ν i is a frequency prefactor, E i is the energy barrier and T is the temperature. Practically, the simulation starts by defining all rates (rate catalog) r i of all possible processes that describe the physical problem. The total rate, R = i r i , is first computed and then a process n is randomly chosen according to, where ξ 1 is randomly chosen in (0, 1) and a single event is performed. The time advances as t = t + t with t being, where ξ 2 is an additional random number chosen in (0, 1). R is recalculated based on the new system state. The algorithm stops when the desired time interval is reached. kMC rates must obey also the detailed balance condition (Equation 2), even if the system is not in equilibrium, to ensure the dynamic evolution will correspond to a physical system [8]. kMC rates depend on both the particle and the lattice type that participate in the process and can be calculated via Transition State and Harmonic Transition State Theories (TST -HTST) [8], Density Functional Theory (DFT) and ab initio methods (e.g., [9][10][11][12]).
Focusing on the deposition processes, the notion of lattice is of great importance. The lattice represents the deposition surface and is composed of sites upon which all events occur, simplifying the construction of the rate catalog. Depending on how the lattice is represented, the atomistic information can either be presented in full detail (e.g., [13]) or in a coarse-grained way where microscopic neighboring sites are coalesced into coarse cells (e.g., [14,15]). Off-lattice kMC [16] methods have also been proposed in atomistic representation where the rate catalog is computed "on the fly" in every step.
In MC/kMC the number and types of processes that occur upon a lattice site differ depending on the physical/chemical phenomena of interest. For deposition processes, adsorption, desorption, surface diffusion, and surface reaction events are used. This set of events is customized to capture the physical and chemical mechanisms unique to each deposition process (Figure 1). The ultimate goal of an MC/kMC model is to describe the interactions of particles with the surface, define the growth rate and predict the profile of the growing film on the surface.

DEPOSITION PROCESSES
During any deposition process the material is deposited upon a surface either by physical or chemical processes. In physical processes, the material is injected as a gas and sticks on the deposition surface with a probability (or rate). Physical Vapor Deposition (PVD), which encompasses sputtering and evaporation, belongs to this category of deposition processes. In chemical processes, the material is grown on the deposition surface through surface reactions. Examples of chemical deposition processes are Chemical Vapor Deposition (CVD), Atomic Layer Deposition (ALD) and electrochemical deposition. 1 | Schematic representation at the molecular level of the basic principles in deposition processes. In PVD, particles are adsorbed in the substrate. In CVD, particles (here an arbitrary molecule is shown) react on the substrate surface to grow the film. In ALD, multiple precursors (here A and B) are injected into the reactor though pulses in a cyclic way of pulses-purges and grow the film through self-limiting surface reactions. In electrodeposition, the voltage (V) is applied causing the particles from the cathode to move to the anode through an electrolyte where the film is grown through adsorption or surface reactions. In MC/kMC methods, an activation energy is assigned to each event. In MC, if the event leads the system to a smaller energy, it is accepted unconditionally. If not, it is selected with a probability according to Equation (2). In kMC, an event is selected from a predefined rate catalog containing all possible rates and performed unconditionally. Then the system evolves in time according to Equation (4). The four basic events -adsorption, reaction, desorption and diffusion-whose combination can describe a deposition process, with their activation energies, E a , E r , E d , and E dif are shown.

Physical Vapor Deposition
PVD is performed in vacuum chambers where the coating material is evaporated-via a method such as electron or laser beam, arc discharge or sputtering-and is guided to the substrate. Upon hitting the substrate, the vapor condenses to form a film. It is used extensively for the production of metallic nanorods. The applications of the nanorods, e.g., integrated circuit, metallic glue, lithium-ion batteries and fuel cells to name a few, demand conformal shell layer with minimized defects. Compared to CVD and ALD, which offer highly conformal and uniform films on nanostructures, PVD entails a lower cost, higher availability of materials, greater growth rates, and no need for high substrate temperatures. PVD is thus the preferred method for large-scale fabrication. During PVD, the basic growth mechanisms are adsorption and surface diffusion, making MC/kMC models very popular to study the formation of nanorods.
In a series of works, Cansizoglu et al. [17,18] and later Yurukcu et al. [19], studied the growth of nanorods via MC. They applied their MC model in PVD of Ag on In 2 S 3 nanorods and were able to assess the favorable conditions that lead to a conformal shell coating (Figure 2A). They concluded that a wider angular distribution of incident atomic flux can provide conformal coatings around nanowires. Wang et al. [20] coupled Computational Fluid Dynamics (CFD) simulations with an MC model to predict the columnar growth of Cu nanorods for different operating pressures of a plasma reactor. Du and Huang [21,22] used kMC to verify their theory for the transition from thin film to nanorod Cu growth. They found a critical coverage condition at which multiple-layer surface growth form for the onset of nanorod growth. Yang et al. [23] created a kMC model for metallic nanorods growth considering 3D Ehrlich-Schwoebel (ES) diffusion effects. They concluded that larger 3D ES barrier leads to multiple steps formation and eventually nanorods formation.
Alongside nanorods works, Pflug et al. [24] studied the growth of ZnO:Al during reactive sputtering by expanding direct simulation MC to include surface reactions. Chernogor et al. [25,26] performed MC computations to study the properties of TiCrN-Mo 2 N-Ni films during arc discharge PVD. Chen et al. [27] proposed a kMC model for predicting the grain boundary during MoS 2 growth. Finally, Evrard et al. [28] proposed a multiscale framework based on an MC model to predict film's thickness as a function of the mass flux reaching the deposition surface during magnetron sputtering.

Chemical Vapor Deposition
CVD is the most widely used deposition process in industry. Applications range from coatings and integrated circuits to powders and nanomaterials. During CVD, a compound (precursor) containing the material to be deposited, is transferred to the substrate for deposition where the film is grown via surface reactions. The deposition of two-dimensional (2D) materials is one of the main applications of CVD in the literature of the past 5 years. 2D materials, such as graphene and transition metal dichalcogenides (TMDs e.g., MoS 2 , WSe 2 , WS 2 ), have attracted considerable interest due to their potential in the semiconductor industry, particularly applied to electronic and spintronic devices. Such devices demand high quality orientationally-ordered films with low nucleation density and low number of defects, properties supported by CVD.
Most published works focus on graphene growth. During CVD, graphene can grow in different shapes, such as hexagons, multiple-lobed islands, snowflake-like structures or hexagons with dendritic edges. The shape and size of these graphene domains are affected by many process and physical parameters such as deposition flux, temperature and state/orientation of the substrate. To predict the growth shape, new diffusion processes, whose activation energies are calculated via DFT, are designed. Gaillard et al. [9] showed that zigzag edges of graphene flakes arise from low fluxes and higher deposition rates and lower temperature conditions produce less compact islands formations. Their predicted shapes are in qualitative agreement with experimental reports for hexagons, multiplelobed and snowflakes with fractal shapes. Additives, such as N 2 , does not affect the shape of graphene flakes [29]. By introducing etching process on their rate catalog, Chen et al. [30,31] predicted dodecagon with a hexagonal hole, double hexagonal rings, honeycomb-like networks and nanoribbons. In a recent work [32], they studied graphene growth on Cu (1 1 1) as a function of deposition flux and temperature. They concluded that at high temperatures, an increase in the deposition flux leads to graphene transforming from a circular shape to an intermediate compact hexagonal shape, before finally transitioning to a fractal with sixfold symmetry. Göltl et al. [33] by combining experiments and kMC computations proposed a qualitative description for the growth of high aspect ratios nanoribbons on a Ge (0 0 1).
Because of the 2D nature of graphene, most works only consider the early stages (monolayer) of growth. The mechanism of graphene growth depends on the specific substrate. Graphene growth on transition metals such as Ni, Co, and Fe, is a precipitate of carbon atoms dissolved in the substrate, making carbon solubility a key indicator for the growth rate. Graphene growth on a Cu substrate is due to the diffusion of carbon atoms on the Cu surface with the latter having low solubility and catalyst activity. Enstone et al. [34] proposed an MC model for graphene growth on a Cu substrate to study the effects of the substrate's roughness on graphene nucleation. They showed that the key factors which determine the island size and formation are roughness amplitude and mobility parameters ( Figure 2B). Conversely, on active catalytic surfaces, such as Ir, Rh, and Ru, the interaction of carbon-substrate atoms is strong. Since the graphene and substrate lattices are of differing structures, the graphene spreads like a carpet, forming moiré patterns of super unit cells of hundreds of atoms. Jiang and Hou [35], used a multiscale "standing-on-the-front" (SOF) kMC [36] approach for graphene growth on an Ir surface. In SOF-kMC the growth kinetics are determined by the attachment and detachment of different carbon clusters at a growth front. Their computations showed that the rate-dominating event for concave growth-front segments is carbon monomer attachment-and five-atom carbon clusters attachment for other segments-leading to time-resolved growth behavior that was consistent with scanning tunneling microscopy experiments. TMDs, such as MoS 2 , WSe 2 , and WS 2 , have gain the attention of the semiconductor industry due to their relatively inert surfaces and thickness-dependent electrical and optical properties [37]. The factors that control the size and shape of the evolving monolayer during CVD of TMDs is still an open issue. The first works [10,38] attempted to connect the reactor operating conditions to the geometric characteristics of TMDs. Both works applied kMC models that take into account the explicit formation of chemical bonds. Govind Rajan et al. [38] developed a kMC model based on the terrace-ledge-kink formalism. They were able to accurately predict the geometric shape, size, and aspect ratio of triangular and hexagonal MoS 2 and WS 2 structures and link them to the CVD reactor operating conditions. Nie et al. [10] performed kMC coupled with DFT computations to study the deposition of a WSe 2 monolayer on graphene. They constructed a phase diagram of regions of triangular compact, fractal and dendritic structures of WSe 2 based on various operating conditions. Yue et al. [37] detailed the growth of WSe 2 grains and Wu et al. [39] included growth anisotropy during the growth of WS 2 on quartz by incorporating local substrate effects in their kMC model. Chen et al. [11] developed a combined DFT-kMC model to elucidate the mechanism for ultrafast growth of regular triangular WSe 2 monolayer with compact edges on Au (111). They concluded that this ultrafast growth was due to fast kink nucleation and ultrafast kink propagation along the edge. Li et al. [12] studied MoS 2 growth and concluded that Mo concentration gradient is the key factor for the morphological evolution of MoS 2 from dendritic shape to compact triangular geometry (Figure 2D).

Atomic Layer Deposition
ALD is a subset of CVD that is based on sequential and self-limiting surface reactions. In contrast to CVD, multiple precursors are cyclically pulsed and purged to grow the desired material. ALD offers exceptional conformality on high aspect ratio micro-structures, high thickness and film composition control [53] which makes it prevalent in modern complementary metal-oxide-semiconductor (CMOS) and dynamic randomaccess memory (DRAM) fabrication processes. Concerning deposition inside micro-structures, MC/kMC methods have been used for over three decades in the context of CVD [4]. Nevertheless, CVD could not provide good step coverage in high aspect ratio structures, for which ALD came as a solution.
Schwille et al. [54,55] proposed a new MC method to study the deposition inside micro-structures. Their direct simulation MC-derived method models only the walls of the geometrical structure, rather than partitioning the entire domain into grid cells. This assumes that deposited particles can only collide with the molecules of the carrier gas and can only interact with the walls of the structure [55]. The simplification drastically decreased the computational demands. Poodt et al. [56] proposed an MC model to study the step coverage inside pores by considering reactor pressure effects through a gas-phase collision model whose frequency is determined by the mean free path of the precursor molecules. They concluded that ALD is in the diffusion-limited regime for various experimental conditions, even when low reactor pressures are used. Cremers et al. [57] used a 3D MC model to compare step coverage in arrays of holes and pillars and concluded that the latter are suitable for fabrication processes where large surface areas are desired, such as in sensors, solar cells, fuel cells, and batteries. Muneshwar et al. [58] used a scalable kMC model to study the effects of parasitic reactions (no-ALD) in the step coverage inside high aspect ratio features ( Figure 2E). They concluded that it is essential that parasitic side surface reactions must be restricted, either by operating a lower temperatures or adjusting dosages.
Furthermore, a series of MC/kMC studies regarding ALD of 2D materials [59,60] for conductive bridge random access memory cells, quantum dots [61] and nanoparticle reactivity [62] were also published. In a series of works, Christofides et al. [63][64][65][66][67][68] conducted multiscale simulations of plasma assisted ALD processes, combining CFD and kMC models to characterize and control the process. Taken a step further, they enhanced their computations with an ANN to characterize the microscopic domain film growth dynamics.

Electrochemical Deposition (Electrodeposition)
During electrodeposition, a thin layer of one metal is coated on top of a different metal. An electrical current within a conductive substrate reduces the cations within an electrolyte, growing them as a thin film. Electrodeposition has become an increasingly popular method of fabrication for the next generation of batteries. Conventional Li-ion batteries are reaching their theoretical energy density limits, pushing researchers to investigate alternative materials. Li metal anodes are conventionally used because of their high energy density. However, Li metals show irregular deposition patterns such as dendrite formations. The irregularities result in poor cycling efficiency, capacity fade, and can cause batteries to short-circuit [69], disqualifying its commercial viability. Sitapure et al. [69] studied the multiscale formation of dendritic structures during Li metal anode electrodeposition, using MD computations for the heterogeneous solid-electrolyte interphase (SEI) and kMC for the growth of the dendrites. This formulation considered the effects of the mechanical properties of the heterogeneous SEI in dendritic formation ( Figure 2C). Vishnugopi et al. [70] constructed a kMC model based on self-diffusion processes without accounting for the existence of SEI. Rather, they considered three types of self-diffusion: terrace diffusion, diffusion away from a step and interlayer diffusion. Their computations showed that, neglecting interlayer diffusion, the growth mode of metal undergoes two transitions as the deposition rate increases. The transitions are from film-type to mossy and from mossy to dendritic. When included in the model, interlayer diffusion favors the deposition of smooth films, even at high deposition rates.
Besides Li anodes metal studies, Carim et al. [1,71] used an MC model along with experiments of photoelectrochemical growth to study the spontaneous growth of highly ordered, nanoscale lamellar morphologies of Se-Te films. In a series of works, Li et al. [72][73][74] applied a 3D kMC model to study the electrodeposition of Cu 2 ZnSnS 4 for photovoltaic applications. Zargarnezhad and Dolati [75] combined experiments and kMC computations for the electrodeposition of Ni. Crevillén-García et al. [76] proposed a kMC model based on Gaussian process (GP) emulation to predict the final shape of the film. The authors concluded that the use of GP greatly accelerates the kMC computations without substantial precision loss.

SUMMARY AND OUTLOOK
MC and kMC models that have been used during the last 5 years to study the physicochemical phenomena of deposition processes were discussed. The various categories of deposition processes (PVD, CVD, ALD and electrodeposition) were introduced and selected applications were highlighted. This work maintains a broad perspective of process types and applications, and contributes to recent reviews regarding the computational approaches for the theoretical investigation of graphene growth [77,78], multiscale modeling in CVD [4] and the implementation of surface reactions in kMC [79].
Though old in conception, MC/kMC methods continue to provide a deep insight of the physical/chemical mechanisms during deposition of novel high-end materials by introducing simple, yet effective, customized processes. MC/kMC are expected to play an important role on the study of Xenes materials [80] which-in contrast to the other 2D materials-can be produced with deposition processes alone. Algorithmically, emergent deep learning techniques provide promising methods to produce detailed and realistic rate catalogs with minimum input from the user. This, combined with parallel processing techniques, is expected to increase the robustness and applicability of kMC methods.

AUTHOR CONTRIBUTIONS
NC and GK had the initial idea for this review. NC and DT wrote the initial version of the manuscript. NC, DT, and GM reviewed the MC/kMC methods. AB contributed in shaping up the general framework of the review in the context of deposition processes. All authors contributed in the literature review and in the final version of the manuscript.