Simulation of the optical coating deposition

2.1 Input data for the simulation

In the classical simulation of the thin-film production, the following data are used as input data: parameters characterizing the beam of deposited particle parameters of the residual gas in a vacuum chamber, substrate temperature, etc. Unfortunately, such important input data as energy and velocity distributions of deposited particles and flux density are hardly available now for the most part of deposition processes.

There are two ways to handle this problem. The first one is based on the development of multi-scale models [3] that are able to simulate interactions of ion beam with target and resulting energy and velocity distribution of atoms arriving at the substrate. The second way is the investigation of the dependence of the simulation results on various parameters of the deposition process. In this case, the characteristic values of the deposition parameters should be taken from the intervals that are typical for the simulated deposition techniques. For instance, for the low-energy deposition techniques like thermal evaporation, the characteristic values of deposited particle energy should be taken around 0.1 eV, while for the high-energy techniques like ion beam sputtering (IBS), the energy of the sputtered atoms may achieve several tens of eV.

2.2 Classical molecular dynamics and Monte-Carlo methods

In frame of the full-atomistic approaches every atom of the simulation cluster is considered as a force center interacting with other atoms in accordance with the chosen expressions for the potential energy of interatomic interaction. In the case of MD simulation processes atoms of the growing film move accordingly to the Newton’s laws. In the case of Monte-Carlo (MC) method atoms move in agreement with the Metropolis scheme [4]. MD and MC algorithms are able to perform numerical experiments with millions of atoms and film thicknesses of more than hundred nanometers [5] at the modern high performance supercomputers. Time intervals of MD and MC modeling are limited by tens of nanoseconds and for this reason these modeling methods are not suitable for the simulation of such long-time processes as the surface diffusion. However, in the very recent years the combined MD and MC atomistic technique providing an essential acceleration of simulation was proposed [6]. In perspective this technique can be applied for optical coatings growth simulation.

Full-atomistic MD and MC approaches for deposition simulation are organized as a step-by-step procedure that, in general, consists of the following steps:

1. Substrate preparation. Both crystalline [7], [8], [9] and glassy or amorphous [10], [11], [12], [13], [14], [15], [16] structures are used as substrates. The annealing procedure can be used to prepare the glassy substrate from the crystalline structure [14].

2. Deposited atoms are injected into random positions above the substrate surface (Figure 1). Initial velocities are directed to the substrate according to the angular distribution of the deposited atom velocities. Velocity values are specified accordingly to the energy distribution of the deposited particles, that is, either discrete [8], [10], [11], [12], [13], [14], [15] or continuous [3].

3. MD (MC) simulation starts from the described configuration. As a rule, the NVT (constant number of particles, volume, and temperature) ensemble with the periodic boundary conditions is used at this simulation step. Also, the NVE (constant number of particles, volume, and energy) ensemble is used for thermal equilibration before NVT simulation [3], [8].

4. Steps 2 and 3 are repeated until the growing film thickness achieves a specified value.

Figure 1:

Simulation box structure.

This scheme was applied to the simulation of the silicon dioxide film deposition [7], [8], [10], [11], [12], [13], [14] as well as titanium dioxide [9], and other material [17] deposition.

The choice of the force field describing the interatomic interactions is one of the key issues in the classical atomistic simulation. Usually, the two-body interatomic potentials [18], [19], [20], [21], [22] are used for the thin-film growth simulation. In this case, the potential energy of the system is represented as the sum of terms depending on the distances between the pairs of interacting force centers. These terms are the electrostatic and the Van-der-Waals contributions. The first is the Coulomb interaction of point charges centered at the atoms. For the Van-der-Waals contribution, the Morse potential [23] was used for silica in Refs. [18], [19], [24], [25] and for TiO₂ [26]. The Lennard-Jones potential was applied for the silicon dioxide thin-film simulation [10–14] as well as for Si-(ZrCu) films [27]. In Ref. [17], the pairwise ionic potential was applied to the study of the Mg–Al–O thin-film deposition.

More sophisticated many-body potentials [28], [29], [30], [31] were also used for the deposition process simulation [16]. The polarizable potentials were proposed to study the deposition of Ti_x O_y clusters with an impact energy of up to 40 eV onto a rutile substrate [32], [33].

The force field parameters are defined either from the quantum chemistry calculations [34], [35], [36] or from the fitting procedure comparing the calculated and experimental condense-phase properties. In Ref. [37], the ab initio-based force field was adapted for the condense-phase simulation and was applied for the simulation of the Ni film formation on the silicon dioxide substrate.

The universal ReaxFF force field based on the parameterization of the density functional theory (DFT) quantum chemistry calculations was developed in the recent year [38], [39], [40]. The ReaxFF has a functional form allowing one to describe the chemical bond formation and dissociation in terms of bond-order parameter depending on the interatomic distances. As the bond formation is an important component of the film growth process, the ReaxFF looks promising for the simulation of the deposition processes.

The following parameters are usually used in classical MD simulations of the silica deposition process [10], [11], [12], [13], [14], [15]: the time step of the MD modeling is 0.5 fs, duration of one deposition cycle (step 3 in the algorithm described in this section) is 4–6 ps for NVT simulation. The duration of the NVE simulation is equal 1 ps [8]. The Berendsen thermostat [41] is used to keep the simulation box temperature constant. The electrostatic component of the interatomic energy is calculated using the particle mesh Ewald [42] method.

2.3 Kinetic Monte Carlo

In the above section, we joined discussions of the MD and MC methods because both methods have similar thermodynamic background and describe processes at the same level of sophistication. The kinetic Monte Carlo (kMC) stays separately and requires specifying the transition rates between the initial and possible final states of the system. These rates should be obtained from the results of the simulation using other methods such as MD and MC, mentioned in Section 2.2.

The kMC approach is specially oriented to the simulation of rare events. This approach is not accurate enough for the description of fast processes such as collisions of high-energy deposited atom with the substrate and the film, but it is suitable for the simulation of long-time processes like the diffusion of deposited atoms on the growing film surface. The kMC algorithms are able to simulate the processes having a duration of hundreds and thousands of seconds in the clusters consisting of millions of atoms [43], [44].

The kMC was applied for the simulation of Cu thin film growth [45], [46]. The diffusion of atoms on the growing surface was taken into account in the frame of the applied model. The Morse potential was used for the calculation of interatomic potential energy with parameters that are defined in Refs. [24], [47]. The film density dependence on the substrate temperature was investigated in Ref. [45]. The multi-scale approach including the kMC simulation level was developed in Ref. [3] and was applied for the simulation of TiO₂ thin-film deposition. The density and roughness dependencies on the number of deposited atoms were investigated. In these simulation experiments, the maximum number of deposited atoms achieved two million.

2.4 Quantum methods

The QC approaches are used for the simulation of optical and electronic thin-film properties. As the QC modeling requires a lot of computational recourses, the simulation cluster dimensions are essentially smaller than those in the case of the MD simulations (Figure 2). The atomic positions in the QC cell should be optimized before the optical and electronic parameter calculations.

Figure 2:

Comparison of the QC and classical MD simulation levels.

In Ref. [3], the quantum chemistry approach was applied for the optical property calculation of the resulting classical MD TIO₂ films. The stoichiometric TIO₂ supercell structure consisting of 96 atoms was prepared from the deposited film. All simulations were performed using the Vienna ab initio simulation package [48]. The electronic correlation and the exchange interaction were considered at the density functional theory level [49] in the frames of the PBE generalized gradient approximation [50] and the hybrid functional HSE06 [51].

The self-consistent charge density-functional-based tight-binding (DFTB) method was used in Ref. [52] for the simulation of the highest occupied molecular orbital–lowest unoccupied molecular orbital (HOMO–LUMO) gap and the electronic density of the states of the TIO₂ clusters.

2.5 Other methods

In Ref. [53], the simplified approach referred to as the discrete approach was proposed. This approach describes the growing film as a set of sites that can be either empty or occupied by atoms [53]. The discrete approach is suitable for simulations of crystalline substrate and films. Several models of the deposition process – random, ballistic, and solid on solid – were considered [53]. In the frames of these models, every injected atom was deposited at a random position of a regular network representing the film surface. The models differ by equations describing the height dependence of a network site on the deposition time. In the simplest case of random deposition without diffusion, the model has an analytical solution predicting the unlimited surface roughness growth with the increase in the number of deposited atoms.

In Ref. [54], the simulation of the SiO₂ thin-film deposition was performed based on the assumption that silicon oxidation on the growing film surface can be described in terms of Brownian motion (three-dimensional random walk) of oxygen atoms. The probability of oxygen atoms moving from one local minimum to another was calculated as exp(−E/k_BT), where E is the activation energy. Moving was performed if the probability was greater than a random value in the interval [0; 1]. This scheme is close to the Metropolis algorithm that is widely used in MC simulations.

The dependence of the nanostructured TiO₂ film thickness on the deposition conditions was studied with the help of the artificial neural network (ANN) and regression method [26]. The ANN was used to determine the relationships between the deposition parameters and the deposited film properties. The author of Ref. [26] concluded that the application of the ANN in combination with MD simulations allows one to reduce the number of experiments.

2.6 Interface between different theory levels

One of the main tasks of the thin-film growth modeling is providing a proper interface between the different simulation levels. The problem is that the most accurate QC level providing the calculations of the key thin-film parameters (refractive index, electronic density gap, and so on) can be applied only to small clusters consisting of no more than several hundreds of atoms. The concept of the virtual coater [55] was proposed to solve the problem. The idea is to combine simulation techniques having different length scales to a multiple scale model. The corresponding approach is outlined in Ref. [55]. The vacuum chamber geometry, ion source, and target parameters were the input data for the first simulation level. The output of this level (deposited atom energy and velocity distribution, density of their flow, charges of the deposited atoms, and so on) were then used as the input data for the classical simulation approaches: MD, MC, and kMC. The structural parameters (density, point defect concentration, porosity, and structural ring distribution) can be used to prepare the simulation cells for the QC level. The first results of the application of this scheme in studying titanium oxide films is presented in Ref. [3].

4.1 Density

The results of the silicon dioxide density calculations using the classical MD simulations with different force fields are summarized in Table 1.

It is seen from Table 1 that the most part of the applied force fields reproduce the experimental value of the silicon dioxide glassy density (this value is equal to 2.2 g/cm³). It is likely that the high density value of 2.58 g/cm³ reported in Ref. [58] for the Beest Kramer and van Santen (BKS) force field [21] is connected with the BKS fitting scheme that is oriented to the calculations of the crystalline cell parameters.

Simulations of the thin-film deposition processes in the frame of the deposited films of silica (DESIL) force field [10–14] show that the deposited film density is somewhat higher than the density of glassy silica. These simulations also show that the film density is higher when the energy of the deposited silicon atoms is higher.

Thermal annealing results in the essential reduction of the densities of the deposited SIO₂ films (see the bottom row of Table 1). This result is in a good agreement with the experimental data [60].

The density profiles of the deposited silicon dioxide films are shown in Figure 3. Simulations were performed in the frame of the DESIL force field for the three different values of silicon atom energy. Silica glass substrate was prepared from α-quartz by heating it from 300 K to the temperature above the quartz melting point followed by quenching in the MD simulation with GROMACS [56]. Every deposition cycle was simulated in the NVT (constant number of atoms, volume, and temperature) ensemble. The deposition temperature was equal to 300 K, and the initial velocities of the oxygen atoms corresponded to 0.1 eV. Essential fluctuations of the film density were observed in the case of low-Si atom energy. The increase in deposition energy results in a smoother density profile.

Figure 3:

Density profiles of the deposited SiO₂ films.

The density profiles of the silicon dioxide thin films for the temperature values in the interval from 3400 K to 5200 K are reported in Ref. [61]. The BKS force field [21] was used for the simulation. The maximum cluster size in Ref. [61] achieved 2.7 nm. The density fluctuations with amplitudes up to 0.2 g/cm³ were observed. The thickness of the transition layer between the film and vacuum achieved 0.5–0.7 nm.

The dependence of the TIO₂ film density on the number of deposited atoms was investigated in Ref. [3] using the kMC approach. The MA pairwise potential [31], including the Coulomb and Buckingham terms, was used. It was found that the density reduces with the growth of the deposited atom number N (the N values of up to ~10⁶ were considered). It was also found that the density reduces from 3.4 to 2.7 g/cm³ with the growth of incidence angle from 0° to 50° in the case of the TI atom energy of 0.1 eV. The density increased approximately for 0.6 g/cm³ when the TI atom energy was increased from 0.1 eV to 10 eV [3].

4.2 Roughness

In Ref. [62], the surface roughness was calculated as a root mean square deviation of the vertical coordinates of the surface atoms. The surface roughness of the deposited TIO₂ films was calculated in Ref. [3] using the kMC approach. In all the cases, the roughness was increased with the growth of the number of the deposited atoms, and it achieved 1.5 nm for N=2 · 10⁶. The roughness was also increased with the growth of the angle of incidence. It was reduced from 1.3 nm to 0.2 nm with an increase in the deposited Ti atom energy from 0.1 eV to 10 eV (angle of incidence equal to zero).

The dependence of the roughness R_h on the SIO₂ thin-film thickness and deposition conditions was investigated in the frame of the DESIL force field in Refs. [10], [11], [12], [13], [14]. To calculate R_h, the horizontal plane of the simulation box was divided into square cells having the same surface areas (see Figure 4). For every cell, a search for the atom with a maximum vertical coordinate h_i was performed.

Figure 4:

Dependencies of thin-film roughness on film thickness for various energies of deposited Si atoms E(Si) and different substrate temperatures T.

The R_h value was then calculated as a rms of the h_i values. The calculated R_h values were in a good agreement with the experimental data [63]. It is seen from Figure 4 that the decrease in the substrate temperature results in the increase in the R_h values in the case of the silicon atom energy of 0.1 eV. The growth of the deposition energy leads to a reduction in the surface roughness.

The theoretical models of the evolution of the surface characteristics including the roughness characteristics are discussed in Ref. [53]. The horizontal plane of the growing film was divided into cells. The equations describing the time dependence of the film height in the frame of every cell were derived. The model taking into account surface diffusion was applied in studying the diamond-like carbon film growth. It was found that the increase in the deposition energy leads to the decrease in the surface roughness.

4.3 Porosity and point defects

The porosity analysis of the TiO₂ films was performed in Refs. [3], [44]. The kMC approach was used to obtain a film cluster consisting of two million deposited atoms. Only the pores that were able to contain at least 15 atoms were taken into account. It was found that the closed pores were located near the growing film surface, and this was considered as the evidence of a rather dense structure of the film. The total volume of all the closed pores having diameters of at least two atomic dimensions was found to be equal to 7% of the total cluster volume. The pore concentration rapidly decreased with the pores size growth.

The porosity of the silicon dioxide films was studied in the frame of the classical MD simulation using the DESIL force field in Refs. [10], [11], [12], [13], [14]. The number of the deposited atoms in the simulation procedure achieved 2.4 · 10⁵.

Only the pores having dimensions of at least 0.1 nm were taken into account. It is seen from Figure 5 that low deposition energy results in the formation of pores with dimensions of up to several nanometers. These pores are not observed in the case of high deposition energy. The concentration of the pores for all dimensions is reduced with the growth of the deposition energy.

Figure 5:

Pores in the SIO₂ films atomistic structure. Values of the deposited Si atom energies are shown.

The simulation of the point defects in the silicon dioxide thin films and fused silica was performed in Ref. [14]. It was found that the main types of point defects in SiO₂ were one- and three-coordinated oxygen atoms O₁ and O₃. The concentrations of these defects were calculated for the two films deposited at different energies of the deposited Si atoms. The influence of annealing on the point defect concentration was also studied. The results are shown in Figure 6. The defect concentration in the ‘as-deposited’ films is higher in the case of low deposition energy. The annealing procedure results in the essential reduction of the point defect concentration for all defect types and for all deposition energies.

Figure 6:

Point defect concentration in ‘as-deposited’ and annealed SIO₂ films.

4.4 Ring statistics

Rings consisting of different numbers of atoms (n-membered rings) are formed in different compounds (SiO₂, TiO₂, B₂O₃, and so on). The relative weights of the n-membered rings is defined as f_n=N_n/N_tot, where N_n is the number of the n-membered rings and N_tot is the total number of rings.

The relative weights for the different n were calculated for the deposited SiO₂ films using the MD approach with the DESIL force field. The crystalline silicon dioxide consists of only six- and eight-membered rings, while in the fused silica and SiO₂ films, rings with different n values are formed. Values of f_n were calculated using the shortest path analysis [64].

The results of the ring statistic calculations are shown in Figure 7. The geometry of the chemical bonds is strained in the rings with n<5 as the values of the interatomic Si–O distances and the valence angles O–Si–O differ from the equilibrium ones. The most strained structure is formed in the case of the low substrate temperature and low energy of the deposited Si atoms.

Figure 7:

Ring statistic for different deposition conditions.

The obtained distributions of the f_n values are in agreement with the experimental results for the fused silica [65], [66]. The detailed comparison of the f_n distribution for the different force fields is performed in Ref. [67]. For all the force fields, the distribution maximum corresponds to the six-membered ring. The relative weight of the strained rings (N<5) is minimal in the case of the Tersoff potential for SiO₂ [30].