Numerical modeling of enhanced reactive oxygen plasma in pulsed laser deposition of metal oxide thin films

3.1 Modeling equations

In terms of electron densities, the continuity equation is defined as follows:

where t is the time, while n e , Γ e , and R e represent, respectively, the electron density, the flux of electrons, and the term source dependent on different chemical reactions included in the deposition process [22].

Drift diffusion equations are used to compute the electron density and average energy of electrons, expressed as follows:

Here, D e and μ e represent the electron diffusion and the electron mobility coefficients, respectively, and E is the electric field.

The energy equation for electrons is given by

Here, ε and Γ ε represent the mean energy of electrons and the flux of electron energy, respectively, and R ε denotes the source term of electron energy because of inelastic collisions. The energy diffusivity, energy mobility, and electron diffusivity are given by

where T e represents the temperature of electrons.

The mass fractions of non-electron species are calculated using the following equation:

where ρ , M n , and u denote the gas density, the average molar mass, and the average velocity of the gas, respectively. J k , w k , and V k are the diffusion flux vector, the mass fraction, and the speed vector, respectively, for non-electron species k . Z k denotes the charge number, D k is the diffusion coefficient, μ k represents the mobility coefficient, and R k is the source term.

Equations for electromagnetic fields are formulated as follows [23]:

In these equations, i , ε 0 , μ 0 , and ρ ext represent, respectively, the imaginary number, the permittivity of free space, the vacuum permeability, and the space charge density expressed on the basis of a chemistry model and given by:

B is the magnetic field vectors and ω = 2 π f denotes the angular frequency, where f represents the electromagnetic field frequency. K is the dielectric tensor written as

where σ p is the plasma conductivity given by

and ω pe is the electron plasma frequency expressed as

Here, n e , m e , ϑ , and e represent, respectively, the period averaged electron density, the electron mass, the total collision frequency, and the electron charge.

Combining Eqs. (9) and (10), the coil current density J ext ( s ) can be calculated as

where s represents the position along the induction coil, and the capacitive current flowing through the coil I c can be derived from the volume integral of the charge on the coil

For neutral background gas, Navier–Stokes equations are as follows:

Mass conservation:

Momentum conservation:

where u , T , ρ , μ , p , I , B , and J represent, respectively, the velocity, the temperature, the gas density, the viscosity, the pressure, the density tensor, the magnetic field, and the induced current density in the plasma.

Species on the surface are determined from their deposition height [24]:

where h , R , M , and ρ b , represent, respectively, the total growth height, the surface rate expression, the molecular weight, and the density of the bulk species.

3.2 Plasma boundary conditions

The electron flow into the reactor wall is determined by:

where v e , th , K B , and n represent, respectively, the thermal electron velocity, the Boltzmann constant, and the normal unit vector. γ p is the coefficient of the secondary emission.

Electrons transfer energy to walls through:

As a result of surface reactions, heavy species lose ions to the wall since the electric field is directed toward the wall [25].

Here, ( E ∙ n ) and ( n ∙ j k ) are the normal components of the wall of the electric field and the heavy species current density, respectively. e k and M ω represent the elementary charge of the species k and the mass of the ionic species, respectively. In addition, the reactor walls are grounded.

3.3 Chemical model

Identification of species and reactions in plasma chemistry is critical. The species and reaction used for inductively coupled oxygen plasma enhanced atomic layer deposition are taken from the study of Christophorou and Olthoff [26]. The chemical model includes electrons and 11 different oxygen species, such as O₂, O, O⁺, O⁻, O₂ ⁺, O₂ ⁻, O^*(¹D), O^*(1S), O₂ ^*(v), O₂ ^*(1Δ), and O₂ ^*(1Σ). In this case, the O^*(1S) and O^*(1D) states are the atomic oxygen exciting states. The vibrationally excited state of O₂ is represented by O₂ ^*(v), which refers to the O245 state. The excited states of molecular oxygen are given by (O₂a1d) and (O₂b1s), which refer to O₂ ^*(1Δ) and O₂ ^*(1Σ), respectively. O₃ is not included in the model since it has a minor effect on these plasma types [27]. A total of 62 reactions are included in the model [24]. The electron impact reaction rates are given in the study of Kropotkin and Voloshon [28,29].

The reactor walls neutralize and reflect all the O⁺, O⁻, and O₂ ⁺ ions in the plasma. In addition, O^*(1D), O^*(1S), O₂ ^*(v), O₂ ^*(1Δ), and O₂ ^*(1Σ) are returned to the ground state after being de-excited and reflected back into the plasma after colliding with the wall. A reaction rate of 0.2 is assumed for atomic oxygen to recombine at the wall to produce O₂ [30].

3.4 Geometry of the discharge model

Figure 2 (right) illustrates the simulated geometry. Plasma is confined in a vacuum chamber mounted by a metal electrode and a five-turn planar copper coil antenna. Using a quartz cylinder (1 cm thick and 3.2 cm in radius), electrical power is inductively coupled to the plasma volume. A grounded electrode with a radius of 4.12 cm and made of stainless steel is used. The quartz cylinder and the lower electrode are separated by a space of 4 cm. Oxygen gas is injected into the discharge volume at 400 W applied power with 13.56 MHz sinusoidal frequency and pressures between 3 and 100 Pa. Figure 2 (left) shows the meshing of the reference cell with 14,460 triangular elements.

Figure 2

Cell coupled GEC reference structure (right) and structure meshing of GEC-ICP (left).

3.5 Computational process

COMSOL Multiphysics is used to implement a 2D axisymmetric simulation model of ICP generated in a GEC reference cell reactor [31]. Models with axial symmetry perform better than their isotropic counterparts. Nonetheless, due to its inability to scale, the axisymmetric-based approach is impractical when dealing with large data sets. As a result, subsampling or providing computationally practical likelihood estimates is necessary [32]. It is typically found that inductively coupled discharges generated at low pressure provide high charge densities (>10¹⁶ m⁻³) [33,34]. Furthermore, surface anisotropy is achieved with low pressure ion bombardment, which makes high density plasma sources appealing. The negative ion temperature is equal to 0.3 eV when calculating negative ion mobility and diffusivity since the ions are held within the plasma core, where they acquire energy from the electromagnetic field at high frequencies [35,36]. Three steps are involved in calculating plasma properties. As a first step, the electromagnetic field inside the discharge reactor must be determined using Maxwell’s equations. In the second step, based on these fields, Boltzmann solvers are used to determine electron energy distribution functions and electron impact reaction rates. In the third step, on the basis of the reaction rates, different species of oxygen plasma gas and distribution of electron density are determined. Finally, electromagnetic fields are re-calculated based on the electrostatic field derived from Poisson’s equation, and the calculation loop is closed. An iterative process is carried out until a model is convergent. The model simulation requires a 64-bit machine with at least 4 GB of memory because of the large number of reactions and species.

4.1 Validation of the numerical model

Figure 3 shows the variation of electron density with the applied power at various pressures 0.5, 7, and 14 Pa (3, 50, and 100 mTorr) in inductively coupled oxygen plasma. The results are compared with the experimental electron density obtained from the literature [38] using a Langmuir probe with uncertainty measurements between 5 and 10%. There is a linear relationship between electron density and power, while a weak decrease of electron density is observed with pressure. The agreement is good; however, the chemical model did not adjust experimental chemical coefficient rates and there are no inelastic electron collisions, which can explain the simple disagreement between simulation and experiment.

Figure 3

Variation of numerical and experimental electron densities with the applied power in RF-inductively coupled oxygen plasma.

4.2 Reactive oxygen plasma

Figure 4 presents the 2D spatial distribution in electron density, electron temperature, collisional power loss, electric potential, and velocity field in the discharge volume at 400 W power, 10 Pa (75 mTorr) pressure, and oxygen gas.

Figure 4

2D distribution of (a) electron density, (b) electron temperature, (c) collisional power loss, (d) electric potential, and (e) velocity field.

The maximum electron density is about 1.21 × 10¹⁸ m⁻³ and is located at the center of the reactor under the RF coil (Figure 4a). Azimuthal electric field shielding occurs and is high in this case because of the density of electrons. The temperature of electrons is displayed in Figure 4b. In the reactor center, the density of electrons and temperature are elevated, and the dissipated power (Figure 4c) is increased at 0.77 × 10⁶ W/m³. Most power deposition occurs at an arc length of 35 mm below the third and fourth coils of the reactor. Figure 4d displays the distribution of electric potential, and Figure 4e presents the velocity field, where gas cannot penetrate deep into the plasma core at these flow rates. Through the pump, the gas flows down toward the wafer, then around it, and out to the pump. Gas flow has a relatively small effect on ICPs which operate at very low pressures.

The surface distribution of reactive oxygen density at 400 W RF input power and 10 Pa (75 mTorr) is shown in Figure 5. The distributions of charged particles (electrons, O⁺, O⁻, and O₂ ⁺) indicate that most of the power is dissipated just below the coil, in front of the quartz cylinder. The predominant positive oxygen ion is O⁺, with a maximum density of 1.17 × 10¹⁸ m⁻³. O₂ ⁺ has a density that is nearly two orders of magnitude lower. Due to the presence of a considerable number of negative ions O⁻, the plasma is electronegative.

Figure 5

Surface distribution of reactive oxygen plasma (electron, O₂, O, O₂ ⁺, O⁺, O₂ ⁻, O⁻, O^*, O₂ ^*(1∆), O₂ ^*(v), and O₂ ^*(1∑)) in 10 Pa oxygen gas at 400 W in a GEC-ICP reactor.

The neutral oxygen atom (O), excited atomic oxygen species (O^*(1D) and O^*(1S)), and excited molecular oxygen species (O₂ ^*(1Δ), O₂ ^*(1Σ), and O₂ ^*(v)) identified in the model as O2a1d, O2b1s, and O245, respectively, have peak densities that are two or three orders of magnitude higher than those of charged particles [39]. In addition, the neutral species distribution across the inter-electrode gap is more homogeneous. These reactive neutral species dominate the impact of charged particles in PE-PLD, emphasizing higher importance in thin film deposition. The singlet delta oxygen O₂ ^*(1Δ) identified as O2a1d and atomic oxygen (O) are chemically highly reactive. So, in PLDs with an O₂ gas background, these species are critical in thin film deposition.

Plasma plumes as well as background gases react to create reactive O and O₂ ^*(1Δ) species. As a result, reactive oxygen species are directly correlated with the ablation process limiting the control of reactive oxygen species properties [40].

O⁻ ions have a higher density than O₂ ⁻ ions, since O₂ is produced by the dissociative attachment reaction (e + O₂ → O⁻ + O) which is greater than the three-body attachment reaction (e + 2O₂ → O₂ + O₂ ⁻).

Figure 6 illustrates the axial distribution of oxygen density at 400 W-RF and varied O₂ pressures of 10 Pa (75 mTorr) and 50 Pa (375 mTorr). In Figure 6a, a doughnut shape is formed by the charged species, which is even more localized in front of the quartz cylinder in the 10 Pa case. In Figure 6b, the reactor-averaged density of all charged species is higher at 50 Pa than at 10 Pa. Reactive neutral oxygen atom species have lower peaks, and reactor average densities are reduced by a factor of 10 orders for O. The maximum oxygen species densities are localized under the quartz window.

Figure 6

Axial density distribution in reactive oxygen gas (electron, O₂, O, O₂ ⁺, O⁺, O₂ ⁻, O⁻, O^*, O₂ ^*(1∆), O₂ ^*(v), and O₂ ^*(1∑)) with 400 W input power for (a) 10 Pa and (b) 50 Pa.

O₂ ^*, O^*, and O are no longer homogeneously distributed across the reactor. They are localized in the same region as the charged species. During the PLD process, the oxygen content is optimized by exploring these differences in density with pressure.

During deposition, oxygen atoms are generated by both the target and collisions in the plasma plume and have a greater oxidation power compared to molecular oxygen [41]. Figure 7 shows the evolution of oxygen atoms O and electron density under various pressures. Figure 7a shows that the electron density varies linearly with pressure. Figure 7b shows that with increasing pressure, the oxygen atom density increases and then decreases.

Figure 7

Radial distribution of (a) electron density and (b) atomic oxygen density for different pressures.

To explain more the variation of atomic oxygen density with pressure, the pressure is varied between 3 and 100 Pa (22–750 mTorr), which is a suitable domain to activate the deposition substrate. Figure 8 displays the density distribution of oxygen atoms and metastable molecular oxygen close to the metal electrode. It indicates that the O₂ ^* density increases from 4 × 10¹³ cm⁻³ to 1 × 10¹⁵ cm⁻³ when the pressure is increased from 3 Pa to 100 Pa. In contrast, O exhibits a different behavior with pressure, varying between 3 × 10¹³ cm⁻³ and 2 × 10¹⁴ cm⁻³, with an elevated density at 10 Pa (75 mTorr) and the lowest at 100 Pa. Hence, the dominant species at 3 Pa is O, which has a density that is twice as high as O₂ ^*. By increasing the pressure to 100 Pa, the dominant reactive species is O₂. This change in reactive oxygen species impacts on thin film properties. A similar result is found in simulations by Blackwell et al. [42], wherein the crystal composition, stoichiometry, and quality of metal oxide films are affected by the O/O₂ ratio.

Figure 8

Density distribution of O₂* and O with pressure in front of the lower electrode.

The plasma medium has a high density of ionized particles which affects its diagnosis. Figure 9 illustrates the radial distribution of charged and metastable species of oxygen gas with 400 W input power at 10 Pa (75 mTorr) and 50 Pa (375 mTorr). The peaks of reactive ion oxygen density increase with pressure. Also, with increasing pressure, the ion drift flux depending on the electron temperature and the diffusion flux of ions decline together, resulting in a decline of the total ion flux.

Figure 9

Radial distribution of (a) oxygen particles and (b) metastable oxygen particles for different pressures at z = 0.02 m with 400 W input power.

Figure 10 shows how charged, neutral, and metastable oxygen species density evolve with applied power. In Figure 10a, electron densities vary linearly with power. This trend is also supported by experiments [43]. According to Figure 10b, increasing the power values under the RF coil increases the density of excited oxygen atoms and oxygen molecules. Increasing power causes oxygen atoms to gain energy, leading to higher atom excitation.

Figure 10

Radial distribution of densities of (a) oxygen particles and (b) neutral oxygen particles for different powers at z = 0.02 m with 10 Pa pressure.

The electron temperature and ion density are important parameters influencing the deposition of oxide thin films. The electron temperature is proportional to bombardment energy, causing oxide films homogeneity defect. Furthermore, ion density affects the deposition rate. Hence, the oxygen ion density and electron temperature in Figure 11 are displayed as functions of power and pressure in the middle of the reactor. Figure 11a illustrates the applied power affecting the ion density and electron temperature at 10 Pa (75 mTorr). As the applied power increases between 100 and 600 W, the oxygen ion density increases and the electron temperature decreases. Increasing power causes higher ionization rates, which result in a plasma rich in electrons and ions. While a decline in the electron temperature is observed in the plasma zone because of energy loss [44,45].

Figure 11

Ion oxygen density and electron thermal distribution in films deposited by PE-PLD under different (a) applied powers and (b) pressures of the ICP.

According to Figure 11b, when the pressure varies from 10 to 40 Pa (75 to 300 mTorr), the oxygen ion density increases, while the electron temperature decreases. Furthermore, the oxygen density in the plasma region becomes higher as the pressure rises. As collisions increase at high pressure, particles per unit volume increase [46,47]. A proportional increase in collision rate with pressure results in energy loss by collisions, which lowers the temperature of electrons. There is an agreement between these results and those of Ojeda et al. [48] in terms of density and electron temperature.

A simulation of oxide film growth rate in response to applied power and gas pressure is presented. Figure 12 shows how applied power affects ZnO film deposition rates for different operating pressures. With increasing power and pressure, the deposition rate increases almost linearly [49].

Figure 12

Variation of deposition rate with applied power at different values of pressure.

4.3 Correlation analysis

This article focuses on the growth of metal oxide thin films on a polymer substrate through PE-PLD to overcome some limitations of the conventional PLD method, such as the need to use multi-element targets and elevated substrate temperatures. The standard PLD setup with a metal target will be combined with an electrically produced low-temperature oxygen background plasma. Crucially, oxygen plasma is a pulsed non-equilibrium plasma to maintain a low temperature, which prevents the substrate from being significantly heated by the plasma’s conductivity and enables the deposition of coatings onto delicate substrates.

The results show that the direct contact between the active plasma and the ablation plume results in more reactive, energetic plasma particles impinging on the substrate, instead of an afterglow plasma beam. In contrast, in the plasma beam assisted PLD system presented by Dinescu et al. [17], an additional oxygen plasma beam impinges on the substrate to form ZnO films. This system provides high-quality ZnO films only when combined with 800 K substrate heating.

In addition, since these particles are more energetic in PE-PLD than an afterglow plasma or a neutral background gas, where electrons present a high temperature proportional to the energy of bombardment of the substrate, their diffusion lengths are longer, resulting in more crystalline films [50].

A background plasma plume also provides particles for a longer period of time than a laser-produced plume. As a result, polycrystalline films can be deposited at room temperature because incoming particles are more likely to find good sites on the surface for crystalline growth. Our results are in accordance with the work of Tricot et al. [51], which concludes that using pulsed-electron beam deposition systems delivers plasma species to substrates for much longer periods of time than conventional PLD systems.

In addition, De Giacomo et al. [52] studied TiO₂ films made via PE-PLD, while Huang et al. [15] reported ZnO films deposited at ambient temperature using the RF-PE-PLD technique. These systems similar to our work use a plasma backdrop to decrease the droplets in the PLD plume rather than using the background plasma to provide oxygen for the film.