Comparison of 1D and 2D particle-in-cell simulations for DC magnetron sputtering discharges

Direct current magnetron sputtering (DCMS) has been developed as a useful technique and successfully applied by various industries for the deposition of metallic and dielectric thin films.¹ When the discharge is driven by high power density and short duty cycle, it is referred to as high power impulse magnetron sputtering (HiPIMS) discharge.^2,3 To deposit high-quality thin films, a thorough understanding of discharge physics and the basic plasma processes is of essential significance. As an important tool, zero-dimensional global model has been adopted in investigating the DCMS discharges,⁴ specifically when the discharge is operated at high power density (i.e., HiPIMS discharges^5–7) and with complex plasma chemical reactions (reactive HiPIMS discharges^8,9) However, the plasma global model cannot provide spatially resolved information. Among other spatially resolved models,^10–15 the analytical model¹⁰ is the least computationally expensive, but its accuracy and applicability are the most limited. Fluid models^11–13,15 are computationally efficient; however, their accuracy is questionable and the applicability is limited under low pressure and strong magnetic field, which are typical conditions for magnetron sputtering discharges.

The disadvantages of the above-mentioned models make particle-in-cell/Monte Carlo collision (PIC/MCC) simulation^16–23 a better choice. Despite the cost of long computation time, PIC/MCC simulations can provide a complete and self-consistent picture of magnetron sputtering discharges. It is commonly believed that due to the three-dimensional magnetic field distribution and its large curvature and strong gradients, the geometry of a planar magnetron cannot be reduced to one dimension (1D), but at least two-dimensional (2D) coordinates in radial and axial directions should be considered, ignoring the azimuthal $E×B$ direction.²⁴ For this reason, most PIC simulations on DCMS discharges are performed in two or three dimensions.^25–35 There are some one-dimensional models for DCMS discharges being developed, including analytical,³⁶ fluid,³⁷ and PIC simulations.^38,39 However, their applicability and accuracy are unclear compared to 2D simulations. If 1D simulations can capture the primary characteristics of magnetron discharges, the high computational load can be greatly reduced when investigating the DCMS discharge qualitatively.

To investigate the applicability and accuracy of 1D PIC/MCC simulations on DCMS discharges, this work performs one-dimensional (1d3v) and two-dimensional (2d3v) PIC/MCC simulations for a typical balanced DCMS discharge. The 1D and 2D Cartesian simulations are based on our ASTRA PIC/MCC code (see the supplementary material in Ref. 40 for the description and the code benchmark with Turner et al.⁴¹). Figure 1(a) illustrates a schematic of the modeled planar magnetron sputtering setup. The 2D simulation region is enclosed by red lines. The plasma is generated between a metal target and chamber walls separated by 3 cm in the $y−direction$ with a reflection plane at x = 0 cm and confined by the right chamber wall at x = 6 cm. For the 1D simulation, the line of x = 2.9 cm is selected as the simulation region. This line is the most representative for balanced DCMS discharges, where the magnetic field lines are parallel to the target surface, and the plasma region is approximately symmetrical with this line as the mid-plane. The metal target, with a thickness of 0.6 cm and a width of 5.5 cm, is connected to a current source which supplies a constant current density of $0.2 mA/cm2$ in the 2D simulation. Note that this current density is the nominal current density, averaged over the entire area of the target surface. Since the plasma is highly nonuniform and concentrated above the “racetrack” area (where the magnetic field lines are parallel to the target surface and the plasma density is the highest), the “effective target area,” which collects the ion current, is generally about one order of magnitude smaller.^42,43 To consider this more physical relevant, “effective current density,” a constant current density of $2 mA/cm2$ is applied in the 1D simulation (ten times compared with that used in the 2D simulation). The correctness of this assumption is verified later. The time-dependent electric potential on the metal target is automatically adjusted by the charge accumulation on the target provided by the current source and the plasma through a self-consistent circuit model.⁴⁴ Figure 1(b) presents the magnetic field distribution in the 2D simulation region. The balanced magnetic field is created by two magnets located under the metal target. The maximum parallel magnetic field B_x is about 240 G. For the 1D simulation, the magnetic field component B_x along x = 2.9 cm is adopted, while the corresponding B_y component is relatively small and ignored. The boundary conditions for the 2D PIC simulation are as follows:

1. The left boundary at x = 0 cm is symmetric, where particles are reflected and the electric field E_x is 0.

2. The top and right boundaries at y = 3.6 cm and x = 6 cm are grounded and the electric potential $φ=0$⁠. Particles reaching these boundaries are absorbed.

3. At the gap at y = 0 cm and $x=5.5−6$ cm, particles are absorbed and the potential varies linearly.

4. The secondary electron emission (SEE) and the electron reflection are only considered on the cathode target surface.

The simulations were performed in argon gas, consisting of electron–neutral collisions (elastic, excitation, ionization) and ion–neutral collisions (elastic, charge exchange) with the same cross section data used in Ref. 45. Only electrons and argon ions are traced as particles. The background gas is spatially uniform and temporally independent, with a temperature of 300 K and a pressure of 5 mTorr. Initially, the electron and ion superparticles are uniformly distributed in space with a density of $1015 m−3$ and a Maxwellian distribution of 2 eV and 0.026 eV, respectively. Note that for DC discharges, when a voltage source is adopted, the addition of a ballast resistor in series is often preferred to limit the discharge current and solve the stability problem.³¹ Here, we use a current source in the simulation and select an initial density that is close to the steady state; otherwise, the high steady-state current density may result in a failure of igniting the discharge due to an overly fast voltage rising. One can also gradually increase the current density from low to high for pre-ionization and obtain the same steady-state results. Ions are advanced once every ten time steps to speed up the calculation. Electrons are reflected at the target surface with a probability of 0.5, while an energy-dependent SEE coefficient induced by argon ions on dirty metal surfaces is adopted as^46,47

where ε is the energy of argon ions in eV. Around $4.6×105$ superparticles are followed in the 2D simulation and around $7×104$ in the 1D simulation. The time step $Δt$ and the grid spacing $Δx$ are set as $1×10−10 s$ and 0.2 mm, respectively. To eliminate the major constraints of grid spacing and time step in PIC simulation, the implicit PIC^18,48–52 has been developed, which damps high frequency modes while maintains most of the kinetic effects. Here we adopt a standard direct implicit algorithm^48,50,52 for all the simulations; therefore, the major constraints on resolving the Debye length and the electron plasma period are eliminated.^51,53 Courant criterion requires $vΔt/Δx≤1$⁠, where v is a characteristic velocity of particles. For the electrons emanated from the target surface and accelerated in the cathode sheath, a short $Δt$ is required to fulfill the Courant criterion, which could result in a long computational time. However, since these fast electrons only account for a small fraction of the entire electron population, $Δt$ is fixed in the simulation and the following procedure is implemented. The electron travel distance, $Δs$⁠, in a time step $Δt$⁠, is checked for each electron. If $Δs>Δx$⁠, the time step $Δt$ for this electron is split and the electron is advanced twice with $0.5Δt$⁠. This procedure iteratively repeats until $Δs<Δx$⁠. The typical convergence time is about tens of microseconds.

The left column of Fig. 2 presents the calculated 2D spatial density distributions of ions, electrons, and space charges. The plasma is well confined by the magnetic field and has a maximum density of about $1.2×1016 m−3$ above the racetrack. The electron density follows the argon ion density in the bulk plasma region (i.e., the quasi-neutral region), leaving an electropositive sheath region with a thickness of approximately 3 mm near the racetrack [Fig. 2(c1)]. The 2D simulation results at x = 2.9 cm are compared with their 1D counterparts in the right column of Fig. 2. The maximum densities calculated by 1D and 2D PIC simulations are very close. In the 2D simulation, the local current density at x = 2.9 cm at the target surface is about $1.85 mA/cm2$⁠, i.e., nearly one order of magnitude higher than the nominal current density of $0.2 mA/cm2$ and close to the current density of $2 mA/cm2$ used in the 1D simulation. Therefore, in 1D and 2D simulations, the similar maximum plasma density above the target racetrack comes from the similar discharge intensity, i.e., the similar local discharge current density. For the 1D case, electrons cannot escape to the target along the magnetic field lines, resulting in a higher plasma density in the bulk plasma region.

Although a higher plasma density can be achieved via the implicit approach, we select a relatively low density for the following reasons. First, when the density is approaching HiPIMS discharges, additional physical processes, such as Coulomb collisions, sputtering wind, ionization of sputtered species, SEE induced by multiply charged metal ions, etc., need to be considered. Second, the cathode sheath width is on the order of sub-millimeter at higher densities, and a very fine mesh is required to achieve good resolution in the sheath area. In this case, nonuniform or unstructured mesh in the sheath region may be a better choice, which requires further investigation and development. More complete physical processes and more advanced simulation techniques will be reported separately. Therefore, a low plasma density is selected to meet the basic physical processes and sufficient sheath resolution without losing the spirit of this research topic.

As illustrated in Fig. 3(a1), the narrow region above the racetrack corresponds to the cathode sheath, which is much wider outside the racetrack. A contour line of 0 V potential appears close to the top wall, indicating a small anode sheath (approximately 0.6 V) between the contour and the grounded wall. Comparing Fig. 3(a1) with (a2), the discharge voltage required for the 1D simulation is about −500 V, higher than −420 V of the 2D simulation. The potential of the metal target is automatically adjusted by the charge accumulation on the target provided by the current source and the plasma, as described in detail by Vahedi and DiPeso.⁴⁴ However, neglecting the difference in absolute value, both potential profiles have a sheath region of about 3 mm and a potential drop across the bulk plasma region with about 15% of the discharge voltage, consistent with previous theoretical and experimental results.^54–57 The ionization rate, as shown in Figs. 3(b1) and 3(b2), reaches a maximum of more than $2.5×1022/(m3·s)$ above the racetrack region. The maximum ionization rate appears at around 3 mm from the racetrack region, where both the electron density and the electron mean energy are relatively high. As shown in Fig. 3(a2), the electric potential of the 1D case is about 20% greater than that of the 2D case. A similar increase also appears in the electron mean energy $εe$ as shown in Fig. 3(c2). The greater $εe$ in 1D case can be primarily linked to the larger voltage drop; additional loss of electrons in 2D simulation to the target surface along the magnetic field lines may play a role as well. In the 2D case, the secondary electrons emanated from the target surface may reach the bulk plasma region and cause nonlocal ionization. However, the high $εe$ in the cathode sheath region is mainly contributed by electrons emitted from the racetrack region; they may not have an influential impact on the DCMS discharge characteristics. The $εe$ in the bulk plasma region is more important.

From Fig. 4(a), the electron mean energy $εe$ is approximately 6–9 eV in the bulk plasma region. $εe$ has the same tendency as the electron density and decreases as it moves away from the target surface, which has been observed decades ago.⁵⁸ The increase in $εe$ as approaching the anode wall has been observed in previous PIC simulations as well.^25,26,39 This increase can be attributed to the repulsion of low energy electrons by the anode sheath.

An interesting phenomenon that has been observed in DCMS discharges is the deviation of the electron energy from the Maxwellian distribution. Previous diagnostic investigations have found that the electron energy probability functions (EEPFs) are Maxwellian like near the cathode target and bi-Maxwellian like further away from the cathode target.^59–62 Figures 4(b)–4(d) illustrate the EEPFs obtained by 1D and 2D simulations at various regions indicated in Fig. 4(a). EEPF with a boosted high-energy tail is observed in all positions. Note that the Coulomb collision is not considered in our simulation. Previous PIC simulations without Coulomb collision^25,27,35 have shown a distribution substantially differed from Maxwellian with a boosted high-energy tail for magnetron discharges at a relatively low plasma density of equal to or less than 10¹⁷ m⁻³. Similar distributions have been experimentally observed in lower density magnetron discharges, such as in the substrate vicinity of DCMS discharges^59,62 or at the initial stage of HiPIMS discharges.⁶³ Our EEPF results are similar to the above-mentioned numerical and experimental results. This non-Maxwellian distribution, with a boosted high-energy tail, is generally referred to as two-temperature⁵⁹ or bi-Maxwellian like distribution.⁶² However, this is usually obtained by fitting a two-temperature model and may not necessarily mean that the electrons are composed of two fully Maxwellianized electron groups. Generally, a PIC simulation without Coulomb collision does not lead to an exact Maxwellian distribution. The similarity between simulation and diagnostic results indicates that Coulomb collisions may not play a significant role in lower density magnetron discharges, and PIC simulations without Coulomb collision are qualitatively sufficient. The strong electric field in the cathode sheath has a strong effect on electron heating, while the electric field in the bulk plasma region is greatly reduced. The maximum electron energy in the cathode sheath can reach hundreds of eV (data not shown here). The relatively high $εe$ in the cathode sheath region as shown in Fig. 3(c2) originates from this high energy tail. Actually, one can obtain the $εe$ at a certain position by integrating the EEPF $fp$ as $εe=∫0∞ε3/2fpdε∫0∞ε1/2fpdε$⁠, and the results would be consistent with the $εe$ that directly calculated by PIC simulation. Since the nonlocal effect is suppressed by the magnetic field, the 1D EEPF outside the cathode sheath changes slowly along the $y−direction$⁠, similar to the change in electric field. For 2D EEPFs, the low energy electrons cool down significantly near the wall, as shown in Fig. 4(d), This is because in the 2D case, x = 2.9 cm is not always the position with the highest potential. As shown in Fig. 3(a1), mid-energy electrons of several eV near the top wall are attracted by a local potential maximum of more than 2 V in the upper-right corner, resulting in a low-energy electron cooling as shown in Fig. 4(d).

Some previous studies have diagnosed and provided a cold electron temperature from 1 eV (Refs. 60 and 62) to 2.7 eV⁶⁴ and a hot electron temperature from 3 eV (Refs. 60 and 62) to 14 eV.⁶⁴ Because the specific distribution of electron energy is influenced by various factors, such as the discharge system structure, the magnetic field, and the discharge intensity, these results are not directly compared here. However, the general characteristics of EEPFs in the low density DCMS discharges, i.e., a boosted high-energy tail in EEPF, are captured by both 1D and 2D PIC/MCC simulations. We speculate that high and low energy electrons are generated at positions with different electric field strengths, i.e., energetic electrons mainly originate from the cathode sheath, and low-energy electrons originate from the bulk plasma region. The separation of high and low energy electron domains may originate from the drastic change in the electric field between the cathode sheath and the bulk plasma region. Similar explanation has been proposed by Sheridan et al.⁵⁹ Sahu et al.⁶⁴ provided another possible mechanism, i.e., the electron heating by Landau damping in the plasma, which has been self-consistently included in the PIC simulation (if the mechanism is effective). However, the details of the explicit mechanisms under different conditions remain to be further studied.

In summary, we demonstrate the accuracy and applicability of one-dimensional PIC/MCC simulations in investigating the balanced DCMS discharges. For a typical DCMS discharge, a PIC simulation study is performed with a current density of $2 mA/cm2$ in 1D and $0.2 mA/cm2$ in 2D to maintain a similar local discharge current density. The basic plasma parameters, such as the plasma density, the potential profile, the ionization rate, and the electron mean energy in the bulk plasma region obtained by both simulations share similar features and magnitudes above the racetrack region. The boosted high-energy tail in EEPF, which is a typical characteristic of DCMS discharges, is also observed in both simulations. The results verify the accuracy of the 1D PIC/MCC simulation and confirm its applicability in qualitatively analyzing the balanced DCMS discharges. One-dimensional PIC/MCC simulations can largely alleviate the computational load while capturing the key physical mechanisms in DCMS discharges.