Modeling and plasma characteristics of high-power direct current discharge

As a surface deposition technique, magnetron sputtering has been widely used in many fields [1] but increasing industrial demands [2] require equipment that is capable of producing films with higher densities and adhesion strength. Compared to neutral atoms, ions have better controllability [3] as the energy and angle can be regulated by an electromagnetic field and consequently, the film structure and properties can be optimized more easily [4–6]. Therefore, it is essential to improve the plasma ionization efficiency in the instruments, for example, by means of high-power impulse magnetron sputtering (HiPIMS) [7, 8]. By introducing high-power pulses with a duty ratio of 0.5%–1.5% to the cathode, HiPIMS can generate high-intensity discharges and highly ionized plasmas. However, the instantaneous high intensity discharge tends to be unstable [3] and arcing occurs frequently, especially in reactive atmosphere. Although the problems can be mitigated by controlling the power and feedback, it is difficult to eliminate them totally. Meanwhile, owing to ion back-attraction [8, 9] and short pulse duration, the deposition rate of HiPIMS is much less than that of DC magnetron sputtering [8]. To improve the discharge stability and deposition rate, pre-ionizations by direct current (DC) [10], radio frequency [11], and medium frequency [12] are used to obtained more stable discharges, or alternate ionization and deposition are realized to obtain higher deposition rate by high and low pulsing techniques [13], modulated pulsed techniques [10], and deep oscillation pulsed techniques [14, 15]. However, the ionization rate is sacrificed and the improvement in deposition rates is limited [16]. Weakening the magnetic field strength on the target surface to reduce restrained electrons and ions can also improve the deposition rates [17–19] but it results in increased target voltage at the same time.

Compared to pulsed discharges, DC discharges are more stable and high-intensity DC magnetron sputtering is expected to increase the deposition rate. However, continuous high-intensity discharges raise the target temperature rapidly and can cause damage to the cathode [8]. In this work, the ionization and deposition rates of C-HPMS are studied by simulation and compared to HiPIMS to explore the practicality of C-HPMS. The results show that the DC discharge power density can be much lower than that of pulsed discharge for similar plasma density and ionization rate, but the DC deposition rate is nearly an order of magnitude larger than that of HiPIMS. Our results provide the theoretical basis for the realization of highly ionized and continuous high-power DC discharges.

In order to accurately calculate the evolution process of each particle in the discharge process and ionization rates of the sputtered materials, a C-HPMS global model is established based on the HiPIMS global model [20]. The simulation domain is the ionization region which is a circular semi-cylinder with the same diameter as the width of the etching ring on the target [21], as shown in figure A1.

In the Ar/Al discharge, the basic particles include cold electrons (e), hot electrons (e^H), ground state atoms (Ar and Al), energetic atoms (Ar^H), metastable atoms (Ar^m), and ions (Ar⁺, Al⁺ and Al²⁺). The particle equilibrium equations are established based on the particle behavior. By adjusting the fitting parameters (Hall parameter ωτ: the ratio of Hall current to discharge current, which is usually between 8 and 30 [22] in a pulsed discharge), the calculated current I_calc (A) is the same as the actual current, indicating that the model reflects the real physical situation [21]. The influence of spatial inhomogeneity of the magnetic field is indirectly included in the Hall parameters. In the HiPIMS global model, fast overflow of hot electrons and space temperature increase are generally ignored due to the short time of the discharge pulse [20, 21, 23]. However, in the C-HPMS, the fast overflow of the hot electrons must be taken into account. In addition, the increased temperature prevents the inward diffusion of the background gas particles from the outside, which cannot be ignored, too.

In the HiPIMS global model, the hot electron equilibrium equation is described as equation (1) [21]:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{n}_{{{\rm{e}}}^{{\rm{H}}}}}{{\rm{d}}t}=\displaystyle \frac{1}{e{U}_{{\rm{D}}}/2}\left(\displaystyle \frac{{U}_{{\rm{SH}}}\cdot {I}_{{\rm{se}}}}{{V}_{{\rm{IR}}}}-{Q}^{{\rm{H}}}\right),\end{eqnarray} \tag{ 1 }$$

where ${n}_{{{\rm{e}}}^{{\rm{H}}}}$ (m⁻³) means the density of hot electron in discharge region. ${V}_{{\rm{IR}}}$ (m³) represent the total volume of discharge region. U_D (V) is the pulse voltage, U_SH = U_D−U_IR (V) is the sheath voltage and U_IR (V) is the presheath voltage. Q^H (W m⁻³) represents the energy loss in inelastic collisions between hot electrons and atoms [24]. The generation of hot electron is the secondary electron current mainly produced by sputtering of Ar⁺ and Al²⁺ and consumption is cooling of ionization collisions. The ionization collision frequency is related to the electron temperature and the relationship is evaluated experimentally [25].

In C-HPMS discharge (duration is 100%), the difference between the generation and the ionization of hot electrons will be amplified and accumulated with the discharge time. Consequently, the consumption of hot electrons is assumed to include two processes: overflow out of the ionization region and ionization collisions [21]. The latter is considered in the HiPIMS global model and so overflow of hot electrons should be added to the DC model. Accordingly, the hot electron equilibrium equation is expressed as equation (2):

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{n}_{{{\rm{e}}}^{{\rm{H}}}}}{{\rm{d}}t}=\displaystyle \frac{1}{{\rm{e}}{U}_{{\rm{D}}}/2}\left(\displaystyle \frac{{U}_{{\rm{SH}}}\cdot {I}_{{\rm{se}}}}{{V}_{{\rm{IR}}}}-{Q}^{{\rm{H}}}\right)-{{\rm{\Gamma }}}_{{{\rm{e}}}^{{\rm{H}}}}\displaystyle \frac{{S}_{{\rm{IR}}}-{S}_{{\rm{T}}}}{{V}_{{\rm{IR}}}}.\end{eqnarray} \tag{ 2 }$$

The overflow flux is:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{{\rm{e}}}^{{\rm{H}}}}\approx \displaystyle \frac{{D}_{{{\rm{e}}}^{{\rm{H}}}}}{R}0.4{n}_{{{\rm{e}}}^{{\rm{H}}}}\left(1+\displaystyle \frac{{U}_{{\rm{IR}}}}{{T}_{{{\rm{e}}}^{{\rm{H}}}}}\right),\end{eqnarray} \tag{ 3 }$$

where R (m) is the half width of race track and ${T}_{{{\rm{e}}}^{{\rm{H}}}}$ (V) is the temperature of hot electrons. The energy distribution function of the secondary electrons is assumed to be Maxwellian with temperature ${T}_{{{\rm{e}}}^{{\rm{H}}}}$ [26]. ${D}_{{{\rm{e}}}^{{\rm{H}}}}$ (m² s⁻¹) is the Bohm diffusion coefficient of electrons, as shown in equation (4)

$$\begin{eqnarray}&&{D}_{{{\rm{e}}}^{{\rm{H}}}}=\displaystyle \frac{1}{{\omega }_{{\rm{ge}}}{\tau }_{{\rm{c}}}\cdot \sqrt{{T}_{{{\rm{e}}}^{{\rm{H}}}}/{T}_{{\rm{e}}}}}\displaystyle \frac{{T}_{{{\rm{e}}}^{{\rm{H}}}}}{B}.\end{eqnarray} \tag{ 4 }$$

B(T) represents the transverse magnetic field of the target surface. The overflow flux and diffusion coefficient of hot electrons can be derived from the electron movement [27]. However, owing to the energy difference between cold and hot electrons [25], the Hall parameter (${\omega }_{{\rm{ge}}}{\tau }_{{\rm{c}}}$) of hot electrons should be multiplied by a parameter: $\sqrt{{T}_{{{\rm{e}}}^{{\rm{H}}}}/{T}_{{\rm{e}}}},$ where ${T}_{{\rm{e}}}$ (V) and ${T}_{{{\rm{e}}}^{{\rm{H}}}}$ (V) represent the electron temperature of cold and hot electrons, respectively. ${T}_{{\rm{e}}}$ (V) can be calculated by the energy balance equation. Since hot electrons are converted from secondary electrons which are accelerated by the sheath voltage, the hot electron temperature ${T}_{{{\rm{e}}}^{{\rm{H}}}}$ (V) is approximately 1/3 of the sheath voltage ${U}_{{\rm{SH}}}$ (V) [21] due to energy loss of collisions [28].

The cathode is assumed to have sufficient cooling capacity and the target does not melt or evaporate. Consequently, the effect of heat is mainly on the average velocity of the two particles Ar and Ar^m (the average temperature or velocity of the other particles are considered in the pulsed model). The particle fluxes (${{\rm{\Gamma }}}_{{\rm{Ar}},{\rm{diff}}}$ and ${{\rm{\Gamma }}}_{{{\rm{Ar}}}^{{\rm{m}}},{\rm{diff}}}$) of Ar and Ar^m and average thermal velocity (${u}_{{\rm{Ar}},300{\rm{K}}}$) are:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{Ar}},{\rm{diff}}}=({N}_{{\rm{Ar}}0}-{N}_{{\rm{Ar}}}){u}_{{\rm{Ar}},300{\rm{K}}}/4,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{{\rm{Ar}}}^{{\rm{m}}},{\rm{diff}}}={N}_{{{\rm{Ar}}}^{{\rm{m}}}}{u}_{{\rm{Ar}},300{\rm{K}}}/4,\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&{{u}_{{\rm{Ar}},300{\rm{K}}}=\left({T}_{{\rm{Ar}},300{\rm{K}}}\cdot \displaystyle \frac{8{k}_{B}}{\pi {M}_{{\rm{Ar}}}}\right)}^{1/2},\end{eqnarray} \tag{ 7 }$$

where ${N}_{{\rm{Ar}}}$ (m⁻³) and ${N}_{{\rm{Ar}},0}$ (m⁻³) are the densities of Ar atoms in and outside the ionization region, respectively. Considering the ionization region temperature, the average velocity of neutral Ar moving outward from the ionization region is:

$$\begin{eqnarray}&&{{u}_{{\rm{Ar}},{\rm{IR}}}=\left({T}_{{\rm{Ar}},{\rm{IR}}}\cdot \displaystyle \frac{8{k}_{B}}{\pi {M}_{{\rm{Ar}}}}\right)}^{1/2},\end{eqnarray} \tag{ 8 }$$

where ${T}_{{\rm{Ar}},{\rm{IR}}}$ (K) is the temperature of the ionization region. The diffusion fluxes of Ar atoms become:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{Ar}},{\rm{diff}}}=({N}_{{\rm{Ar}}0}{u}_{{\rm{Ar}},300{\rm{K}}}-{N}_{{\rm{Ar}}}{u}_{{\rm{Ar}},{\rm{IR}}})/4,\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{{\rm{Ar}}}^{{\rm{m}}},{\rm{diff}}}={N}_{{{\rm{Ar}}}^{{\rm{m}}}}{u}_{{\rm{Ar}},{\rm{IR}}}/4,\end{eqnarray} \tag{ 10 }$$

where the temperature of Ar entering the ionization region is 300 K and the temperature of the ionization region is calculated by the three-dimensional heat transfer field model for different power.

To verify the effects of hot electron accumulation and temperature rise on the discharge process, the Ar/Al discharge is simulated by four different models: (1) pulsed global model, (2) temperature modified model, (3) hot electron modified model, and (4) both hot electron and temperature modified models. Generally, the DC discharge becomes stable at about 300 and so 2000 μs is used to ensure that the discharge is stable so that the model can be fitted with the calculated current. The current and voltage used for fitting are obtained by experiments in which the Al target with dimensions of 121 mm × 294 mm × 12 mm is used in the DC discharge, discharge power is 5 kW, 10 kW, 15 kW, and 20 kW, respectively, and working pressure is 0.8 Pa.

The current evolution obtained by simulation and experiments versus discharge time is shown in figure 1. The current calculated by the original pulsed global model cannot converge in 2000 μs and continues to increase exponentially. At this time, the density of hot electrons is on the order of 10¹⁸ m⁻³, which is more than 10 times that of cold electrons, as shown in figure 1(a). After temperature modification, the evolution of the discharge current changes from being exponential to logarithmic within 2000 μs and the discharge process is close to equilibrium as shown in figure 1(b), indicating that temperature modification is indeed helpful to balance the discharge balance. However, it is not the key because extending the simulation time and adjusting the Hall parameter still cannot match the calculated and actual currents. Figure 1(c) shows that the discharge current reaches a stable state at about 1000 μs after some fluctuations and the numerical value is consistent with the actual current indicating the critical effect of hot electrons in the balance. After modifying the temperature and hot electrons, the discharge current converges rapidly before 300 μs and the fitted result is optimal as shown in figure 1(d). The convergence time is further shortened by increasing the power, revealing that the modified global model has a better convergence.

Zoom In Zoom Out Reset image size

The temperature increase affects plasma diffusion, ion sputtering, and collisions. Both the hot electron and temperature modified models are implemented to analyze the 20 kW DC discharge of the Al target. It is assumed that the temperature in the ionization region is 300–2000 K with an interval of 50 K. The evolution of the particle density and ionization with temperature is derived and presented in figure 2. For neutral atoms, the total Ar atom density decreases with increasing temperature because the Ar supplement from the outside to the ionization region decreases and Ar^m overflow increases. Meanwhile, the density of metal particles (Al and Al⁺) decreases because the high temperature increases the average free path of the particles to weaken sputtering and ionization. However, the ionization rate of metal increases by more than 10% as the temperature goes up from 300 to 2000 K. Owing to the high temperature in C-HPMS, the gas is rarefied before the target, resulting in that the collision occurs mostly between electrons and metal atom whose first ionization energy is much lower than that of the gas neutrals. Thus, the energy loss decreases, leading to an increase of the average electron temperature T_e (V).

Zoom In Zoom Out Reset image size

According to the voltage-current expression of magnetron sputtering by Westwood [29]:

$$\begin{eqnarray}&&J=\beta ^{\prime} {(V-{V}_{0})}^{2},\end{eqnarray} \tag{ 11 }$$

where J $({\rm{A}}\,{{\rm{m}}}^{-2})$ is the current density, β' is the gas constant of the discharge system, V (V) is the voltage, and ${V}_{0}$ (V) is the ignition voltage, the voltage–current relationship of the abnormal magnetron glow system is quadratic. Therefore, the high-power discharge V–I curve can be obtained by the low-power discharge characteristics [30] obtained by experiments involving Cu and Al targets with a transverse magnetic field intensity on the target surface about 40 mT. The discharge race track width is 5 cm and the length of the race track is 51 cm. Thus the race track area is about 255 cm². The input power is 1 kW, 5 kW, 10 kW, 15 kW, 20 kW, 25 kW, and 30 kW, respectively (provided by AE Ascent 60 K), and the working pressure is 0.8 Pa as shown in table 1. The high-power discharge V–I curves shown in figure 3 are obtained from the fitting of the table 1 data and the high-power voltage/current data shown in table 2 are obtained from the fitting of the curves in figure 3 at an interval of 5 kW.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Using the parameters in table 2, the high-power DC discharges of Al and Cu can be simulated with both the hot electron and temperature modified global models. Figure 4 shows the calculation results including the evolution of the Hall current and Hall parameter versus discharge power [22, 31]. The Hall current generated by the cyclotron motion of electrons on the target surface represents the discharge intensity of the planar cathode. It increases almost linearly with discharge power on both the Al target and Cu target. The Cu target shows a larger slope than the Al target, indicating that the discharge intensity increases faster. The Hall parameter of both targets fluctuates when the power is small and then increases gradually and stabilizes. The stabilized Hall parameter of the Al target is about 11.5 and that of Cu is over 20 when the power density is 400 W cm⁻², demonstrating that the Cu discharge is more intense. The sputtering yields of Al and Cu and particle energies calculated by SRIM [32] are shown in table 3. Cu has a larger sputtering yield than Al. The first ionization energy of Cu is 7.4 eV suggesting that Cu can be easily ionized and so the Cu discharge is strong. Although the first ionization energy of Al is smaller (5.8 eV), the sputtering yield of Al is less than that of Cu. Under the high power conditions, the metal discharge is dominant in the plasma [33] and the intensity of Al is less than that of Cu.

Zoom In Zoom Out Reset image size

The densities of the typical particles in the plasmas of the Al and Cu discharges as a function of power are calculated and presented in figure 5. The electron density increases with discharge power and the plasma densities of Al and Cu are on the order of 10¹⁹ m⁻³ for power densities of 183 and 117 W cm⁻². The density of Al atoms shows a small increase when the power density is low and maintains at 10¹⁸ m⁻³ when the power density exceeds 50 W cm⁻², whereas the densities of Al⁺ and Al²⁺ increase continuously with power density. This phenomenon is mainly caused by the small sputtering yield and small ionization energy, inducing the increased sputtering atoms produced by higher power can be further ionized completely. Different from Al, both the densities of Cu and Cu⁺ continue to increase rapidly with increase power but that of Cu²⁺ changes relatively slowly up to the maximum before declining. As the power density goes up, more Cu atoms are sputtered and most of them are ionized because of the small first ionization energy of Cu. However, the large amount of first ionization arises from enhanced collisions between neutral atoms and electrons in the discharge causing the electrons to cool [23, 34]. The ionization energy of Cu⁺ to Cu²⁺ is higher and therefore, it is almost impossible for the electron energy to reach the condition of divalent ionization resulting in a smaller Cu²⁺ density. The density curve of Cu⁺ in figure 5(b) is almost coincidental with the electron density curve when the power density is large, suggesting that Cu⁺ ions are the dominant cations in the Cu discharge and sputtering is mainly self-sustained.

Zoom In Zoom Out Reset image size

Figure 6 shows the ionization rates of metal sputtering in the ionization region calculated by the particle densities. Here, α represents the metal ionization rate and β stands for the probability of ions returning to the target. The metal ionization rate increases logarithmically with discharge power density. The ionization rate of Al exceeds 90% when the power density is 400 W cm⁻² but that of Cu is only about 70%. For comparison, HiPIMS with a pulse width of 300 μs and peak power density of 564 W cm⁻² are simulated by the pulsed global model as shown in figures 6(c) and (d). After stabilization of the HiPIMS discharge, the ionization rates of Al and Cu are 91.9% and 74.6%, respectively. Considering the instability of HiPIMS, DC and HiPIMS are assumed to have the same ionization levels when the difference is less than 10%, as shown by the purple region in figures 6(a) and (b). When the ionization rates of Al and Cu in the high-power DC discharges reach 90% of that of the HiPIMS discharge, the power densities are only 183 W cm⁻² and 176 W cm⁻², respectively, indicating that the power density required by C-HPMS is much less than that of HiPIMS for similar metal ionization rates. Due to the temperature increase of C-HPMS, the gas is rarefied and the proportion of metal atoms increases, resulting in a larger ionization collision probability of metal atoms. Besides, according to figures 6(a) and (b), the higher the power density, the slower the growth rate of the ionization rate. Thus, the ionization of C-HPMS can catch up that of HiPIMS with the increase of power density. With increasing discharge power densities, the probability of ion returning for Al⁺ decreases slowly but that of Cu⁺ increases initially, reaches a maximum, and then stabilizes. The increase in the probability of ion returning is smaller than that of the ionization rate, demonstrating that more ions diffuse outward from the ionization region. Therefore, it can be concluded that the effects of discharge power on deposition rate are positive.

Zoom In Zoom Out Reset image size

The total diffusion fluxes of the metal ions at the edge of the ionization region are used to determine the deposition rates. The ratio of the deposition rates for different DC power to the HiPIMS (peak power is 564 W cm⁻²) discharge is obtained by simulation and shown in figure 7 This ratio increases rapidly with discharge power densities. For the DC discharge at 183 W cm⁻² (Al) and 176 W cm⁻² (Cu), the metal ionization rate of DC is close to that of HiPIMS. The total diffusion fluxes of Al and Cu are 26 and 30 times those of HiPIMS, respectively, and increase rapidly as the power density is increased. The results indicate that not only the DC discharge can provide highly ionized plasmas like HiPIMS with a smaller power density, but also the deposition rates of DC discharge are one order of magnitude higher than those of HiPIMS because of the longer duration.

Zoom In Zoom Out Reset image size

A Cu target is used to carry out the C-HPMS by a homemade cathode under a powerful cooling, the glow pictures and I–V parameters are presented in figure 8. The glow becomes brighter significantly with increasing of the power. The maximum power density can reach 235 W cm⁻² supplied by a high-power AE device, under which the glow is bright green indicating a stable and strong metal discharge. All the data reveals the feasibility to develop a long duration high-power glow discharge to obtain the high ionization and rapid deposition simultaneously.

Based on the pulsed global model for HiPIMS, a C-HPMS global model that incorporates the effects of hot electron accumulation and temperature rise is established. The model adds overflow of hot electrons and overflow of gas atoms at the high temperature. According to the voltage–current characteristics obtained from the low-power V–I data, the plasma discharges are simulated using Al and Cu as examples. The results show that similar plasma densities and ionization rates can be obtained by DC discharges but the power density of 183 W cm⁻² (Al) and 176 W cm⁻² (Cu) are smaller than that of HiPIMS (peak power is 564 W cm⁻²). Metals such as Al with smaller sputtering yields and ionization energy show weaker discharge intensity but higher ionization rates than those with larger sputtering yields and ionization energies such as Cu. The C-HPMS increases the effective sputtering time, resulting in the deposition rates of metal ions increased by more than 10 times. The results described in this paper provide the theoretical basis for the realization of high ionization and rapid deposition of coatings and insights into the development of deposition equipment for C-HPMS.

This work was financially supported by Guangdong Science and Technology Research Grants (No. 2015B090927003), Shenzhen Science and Technology Research Grants (JCYJ20150828093127698 and JCYJ20170306165240649), Hong Kong Innovation and Technology Fund (ITF) (No. ITS/452/17FP), Postdoctoral Innovative Talent Support Program (No. BX20190001), as well as Guangdong—Hong Kong Technology Cooperation Funding Scheme (TCFS) (GHP/085/18SZ).

The schematic of the particle equilibrium reactions in the ionization region is shown in figure A1. The simulation area is approximately a semicircle. The subscripts are defined as follows: diff (diffusion), ex (excitation), Dex (de-excitation), iz (ionization), P (Penning ionization), chexc (charge exchange), coll (collision), recom (recombination), sput (sputtering) and modification of the hot electrons and temperature is highlighted in red in the figure.

Zoom In Zoom Out Reset image size

Table A1.  Parameters in the global model.

Symbol Meaning Unit

${\omega }_{{\rm{ge}}}{\tau }_{{\rm{c}}}$	Hall parameter (fitting parameter)
B	Average magnetic field at the ionization region boundary	mT
ID	Experiment discharge current	A
UD	Experiment discharge voltage	V
${I}_{{\rm{s}}{\rm{e}}}$	Secondary electron current,	A
USH	Sheath voltage.	V
Icalc	Simulated discharge current	A
P	Pressure	Pa
T	Gas temperature	K
RIR	The radius of the ionization region	M
ST	The contact area of the ionization region with the target	m2
SIR	Ionization region superficial area	m2
VIR	Ionization region volume	m3
k	Reaction coefficient	m3 s−1
γ	The secondary electron emission coefficient
Y	Sputtering yield
Te	Electron temperature	eV
n	Number density	m−3
α	Ionization rate
β	Ion return probability
β'	Magnetic control system gas constant

The particle density is described by considering the generation and loss of each particle. Based on the model described in our previous work [30], the model is modified by hot electron accumulation and vacuum chamber temperature increase and is calculated according to process shown in figure A2. Table A1 presents the parameters in the global model, in which n (m⁻³) is the density, m (kg) is the mass, and Γ (m⁻² s⁻¹) is the flux of the corresponding species. k is the reaction coefficient that can be calculated from table A2.

Table A2.  Rate coefficients of the Ar/Cu and Ar/Al plasmas.

Reaction Rate coefficient (m³ s⁻¹) Threshold (eV) References

${\rm{e}}+{\rm{A}}{\rm{r}}\to {\rm{A}}{{\rm{r}}}^{+}+2{\rm{e}}$	${k}_{{\rm{i}}{\rm{z}}}=2.3\times {10}^{-14}{T}_{{\rm{e}}}^{0.59}\exp (-17.44/{T}_{{\rm{e}}})$	15.76	[20]
${{\rm{e}}}^{{\rm{H}}}+{\rm{A}}{\rm{r}}\to {\rm{A}}{{\rm{r}}}^{+}+2{\rm{e}}$	${k}_{{\rm{i}}{\rm{z}}}^{{\rm{H}}}=8\times {10}^{-14}{T}_{{\rm{e}}}^{0.16}\exp (-27.53/{T}_{{\rm{e}}})$	15.76	[2, 23]
${\rm{e}}+{\rm{A}}{\rm{r}}\to {\rm{A}}{{\rm{r}}}^{{\rm{m}}}+{\rm{e}}$	${k}_{{\rm{e}}x}=2.5\times {10}^{-15}{T}_{{\rm{e}}}^{0.74}\exp (-11.56/{T}_{{\rm{e}}})$	11.56	[20]
${{\rm{e}}}^{{\rm{H}}}+{\rm{A}}{\rm{r}}\to {\rm{A}}{{\rm{r}}}^{{\rm{m}}}+{\rm{e}}$	${k}_{{\rm{e}}x}^{{\rm{H}}}=3.84\times {10}^{-14}{T}_{{\rm{e}}}^{-0.68}\exp (-22.32/{T}_{{\rm{e}}})$	11.56	[2, 23]
${\rm{e}}+{\rm{A}}{\rm{r}}\to {\rm{A}}{\rm{r}}+{\rm{e}}$	$\begin{array}{l}{k}_{{\rm{e}}{\rm{l}}}=2.336\times {10}^{-14}{T}_{{\rm{e}}}^{1.609}\\ \times \exp (0.0618{({\rm{ln}}{T}_{{\rm{e}}})}^{2}-0.1171{({\rm{ln}}{T}_{{\rm{e}}})}^{3})\end{array}$		[27]
${\rm{e}}+{\rm{A}}{{\rm{r}}}^{{\rm{m}}}\to {\rm{A}}{\rm{r}}+{\rm{e}}$	${k}_{{\rm{d}}{\rm{e}}{\rm{x}}}=4.3\times {10}^{-16}{T}_{{\rm{e}}}^{0.74}$	−11.56	[20]
${\rm{e}}+{\rm{A}}{{\rm{r}}}^{{\rm{m}}}\to {\rm{A}}{{\rm{r}}}^{+}+{\rm{e}}$	${k}_{{\rm{m}}{\rm{i}}{\rm{z}}}=6.8\times {10}^{-15}{T}_{{\rm{e}}}^{0.67}\exp (-4.2/{T}_{{\rm{e}}})$	4.2	[20]
${{\rm{e}}}^{{\rm{H}}}+{\rm{A}}{{\rm{r}}}^{{\rm{m}}}\to {\rm{A}}{{\rm{r}}}^{+}+{\rm{e}}$	${k}_{{\rm{m}}{\rm{i}}{\rm{z}}}^{{\rm{H}}}=5.7\times {10}^{-13}{T}_{{\rm{e}}}^{-0.33}\exp (-6.82/{T}_{{\rm{e}}})$	4.2	[35]
${\rm{e}}+{\rm{C}}{\rm{u}}\to {\rm{C}}{{\rm{u}}}^{+}+2{\rm{e}}$	$\begin{array}{l}{k}_{{\rm{M}}{\rm{i}}{\rm{z}}}={k}_{{\rm{M}}{\rm{i}}{\rm{z}}}^{{\rm{H}}}=3.898\times {10}^{-14}\\ \times {T}_{{\rm{e}}}^{0.484}\exp (-7.1344/{T}_{{\rm{e}}})\end{array}$	7.68	[36]
${\rm{e}}+{\rm{C}}{{\rm{u}}}^{+}\to {\rm{C}}{{\rm{u}}}^{2+}+2{\rm{e}}$	${k}_{{\rm{M}}{\rm{iz}}2}={\rm{2}}{{\rm{.763}}}^{-{\rm{19}}}{{T}_{{\rm{e}}}}^{{\rm{2}}}+1{{\rm{.382}}}^{-{\rm{16}}}{T}_{{\rm{e}}}-{\rm{3}}{{\rm{.3}}}^{-{\rm{15}}}$(Te > 36 eV)	19.57	[36]
${\rm{e}}+{\rm{C}}{\rm{u}}\to {\rm{C}}{\rm{u}}+{\rm{e}}\begin{array}{l}\left({\rm{elastic}}\right)\end{array}$	${k}_{{\rm{M}}{\rm{e}}{\rm{l}}}=1.03\times {10}^{-11}{T}_{{\rm{e}}}^{-0.607}\exp (-8.98/{T}_{{\rm{e}}})$		[37]
${\rm{A}}{{\rm{r}}}^{{\rm{m}}}+{\rm{C}}{\rm{u}}\to {\rm{A}}{\rm{r}}+{\rm{C}}{{\rm{u}}}^{+}+{\rm{e}}$	${k}_{p}=2.36\times {10}^{-16}$	−4.46	[38]
${\rm{A}}{{\rm{r}}}^{+}+{\rm{C}}{\rm{u}}\to {\rm{A}}{\rm{r}}+{\rm{C}}{{\rm{u}}}^{+}$	${k}_{{\rm{c}}{\rm{h}}{\rm{e}}{\rm{x}}{\rm{c}}}=2\times {10}^{-16}$		[39]
${\rm{e}}+{\rm{A}}{\rm{l}}\to {\rm{A}}{{\rm{l}}}^{+}+2{\rm{e}}$	${k}_{{\rm{M}}{\rm{iz}}2}={\rm{1.3467}}\times {{\rm{10}}}^{-{\rm{13}}}{{T}_{{\rm{e}}}}^{0.3576}\exp \left(-{\rm{6.7829}}/{T}_{{\rm{e}}}\right)$	5.57	[21]
${\rm{e}}+{\rm{A}}{{\rm{l}}}^{+}\to {\rm{A}}{{\rm{l}}}^{2+}+2{\rm{e}}$	${k}_{{\rm{M}}{\rm{iz}}2}={\rm{2.34}}\times {{\rm{10}}}^{-{\rm{14}}}{{T}_{{\rm{e}}}}^{0.59}\exp {\left(-{\rm{17.44}}/T\right)}_{{\rm{e}}}$	18.17	[21]
${\rm{e}}+{\rm{A}}{\rm{l}}\to {\rm{A}}{\rm{l}}+{\rm{e}}$	${k}_{{\rm{M}}{\rm{e}}{\rm{l}}}={10}^{-0.0104{({\rm{ln}}{T}_{{\rm{e}}})}^{2}+0.1134{\rm{ln}}{T}_{{\rm{e}}}-11.7}$		[20]
${\rm{A}}{{\rm{r}}}^{{\rm{m}}}+{\rm{A}}{\rm{l}}\to {\rm{A}}{\rm{r}}+{\rm{A}}{{\rm{l}}}^{+}+{\rm{e}}$	${k}_{p}=5.9\times {10}^{-16}$	−5.6	[20]
${\rm{A}}{{\rm{r}}}^{+}+{\rm{A}}{\rm{l}}\to {\rm{A}}{\rm{r}}+{\rm{A}}{{\rm{l}}}^{+}$	${k}_{{\rm{c}}{\rm{h}}{\rm{e}}{\rm{x}}{\rm{c}}}=1\times {10}^{-15}$		[20]

A.1. Particle balance for Ar and metal species

M stands for metal in the following equations. For the working gas Ar, the rate equation is:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{d}}{n}_{{\rm{A}}{\rm{r}}}}{{\rm{d}}t}=-({k}_{{\rm{i}}{\rm{z}}}{n}_{{\rm{e}}}+{k}_{{\rm{i}}{{\rm{z}}}^{{\rm{H}}}}{n}_{{{\rm{e}}}^{{\rm{H}}}}){n}_{{\rm{A}}{\rm{r}}}\\ \,\,-\,({k}_{{\rm{e}}{\rm{x}}}{n}_{{\rm{e}}}+{k}_{{\rm{e}}{{\rm{x}}}^{{\rm{H}}}}{n}_{{{\rm{e}}}^{{\rm{H}}}}){n}_{{\rm{A}}{\rm{r}}}\\ \,\,-\,{{\rm{\Gamma }}}_{M,{\rm{c}}{\rm{o}}{\rm{l}}{\rm{l}}}\displaystyle \frac{{{m}}_{M}}{{{m}}_{{\rm{A}}{\rm{r}}}}\displaystyle \frac{{n}_{{\rm{A}}{\rm{r}}}}{{n}_{{\rm{A}}{\rm{r}}}+{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}}\displaystyle \frac{{S}_{{\rm{I}}{\rm{R}}}-{S}_{T}}{{V}_{{\rm{I}}{\rm{R}}}}\\ \,\,+\,{k}_{{\rm{c}}{\rm{h}}{\rm{e}}{\rm{x}}{\rm{c}}}{n}_{{\rm{A}}{{\rm{r}}}^{+}}{n}_{M}+{k}_{p}{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}{n}_{M}\\ \,\,+\,({k}_{{\rm{d}}{\rm{e}}{\rm{x}}}{n}_{{\rm{e}}}+{k}_{{\rm{d}}{\rm{e}}{{\rm{x}}}^{{\rm{H}}}}{n}_{{{\rm{e}}}^{{\rm{H}}}}){n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}+{{\rm{\Gamma }}}_{{\rm{A}}{\rm{r}},{\rm{d}}{\rm{i}}{\rm{f}}{\rm{f}}}\displaystyle \frac{{S}_{{\rm{I}}{\rm{R}}}-{S}_{T}}{{V}_{{\rm{I}}{\rm{R}}}}.\end{array}\end{eqnarray} \tag{ A1 }$$

Considering the effect of heat on particle overflow, the diffusion flux of Ar is ${{\rm{\Gamma }}}_{{\rm{A}}{\rm{r}},{\rm{d}}{\rm{i}}{\rm{f}}{\rm{f}}}=\tfrac{{n}_{{\rm{A}}{\rm{r}},0}{u}_{{T}_{{\rm{A}}{\rm{r}}},300{\rm{K}}}\,-{n}_{{\rm{A}}{\rm{r}}}{u}_{{T}_{{\rm{A}}{\rm{r}}},{\rm{IR}}}}{4}.$ ${n}_{{\rm{A}}{\rm{r}},0}$ is the density of Ar outside the ionization region calculated by the Clapper equation. ${n}_{{\rm{A}}{\rm{r}}}$ is the Ar density calculated by the model. The average velocities of Ar particles in different temperatures regions are ${u}_{{\rm{A}}{\rm{r}},300{\rm{K}}}=\sqrt{\tfrac{8{k}_{B}{T}_{300{\rm{k}}}}{\pi {m}_{{\rm{A}}{\rm{r}}}}}$ and ${u}_{{\rm{A}}{\rm{r}},{\rm{IR}}}=\sqrt{\tfrac{8{k}_{B}{T}_{{\rm{IR}}}}{\pi {m}_{{\rm{A}}{\rm{r}}}}},$ respectively. ${{\rm{\Gamma }}}_{M,{\rm{c}}{\rm{o}}{\rm{l}}{\rm{l}}}$ is the collision part of the metal and metal ion sputtering flux calculated by ${{\rm{\Gamma }}}_{M,{\rm{c}}{\rm{o}}{\rm{l}}{\rm{l}}}=\left({{\rm{\Gamma }}}_{M,{\rm{s}}{\rm{p}}{\rm{u}}{\rm{t}}}+{{\rm{\Gamma }}}_{{M}^{+},{\rm{s}}{\rm{p}}{\rm{u}}{\rm{t}}}+{{\rm{\Gamma }}}_{{M}^{2+},{\rm{s}}{\rm{p}}{\rm{u}}{\rm{t}}}\right){F}_{{\rm{c}}{\rm{o}}{\rm{l}}{\rm{l}}},$ wherein ${{\rm{\Gamma }}}_{M,{\rm{s}}{\rm{p}}{\rm{u}}{\rm{t}}}=\tfrac{{n}_{M}{u}_{{T}_{M}}}{4},$ ${{\rm{\Gamma }}}_{{M}^{+},{\rm{s}}{\rm{p}}{\rm{u}}{\rm{t}}}=\tfrac{{n}_{{M}^{+}}{u}_{{T}_{M}}}{4}$ and ${{\rm{\Gamma }}}_{{M}^{2+},{\rm{s}}{\rm{p}}{\rm{u}}{\rm{t}}}=\tfrac{{n}_{{M}^{2+}}{{\rm{u}}}_{{T}_{M}}}{4}$ are the overflow flux of metal and metal ions. u_TM is the root-mean-square velocity of metal particle. ${F}_{{\rm{c}}{\rm{o}}{\rm{l}}{\rm{l}}}=1\,-\exp (-{{\rm{R}}}_{{\rm{I}}{\rm{R}}}{\sigma }_{{\rm{A}}{\rm{r}}-M}({n}_{{\rm{A}}{\rm{r}}}+{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}+{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{H}}}}))$ is the collision probability with Ar in the ionization region. ${\sigma }_{{\rm{A}}{\rm{r}}-M}\,=\pi {\left({a}_{{\rm{A}}{\rm{r}}}+{a}_{M}\right)}^{2}$ is the estimated cross-section of Ar-metal collisions based on the hard sphere model, a_Ar = 1.88 Å, a_Cu = 1.4 Å and a_Al = 2.51 Å are the atomic radii of Ar, Cu and Al.

For Ar⁺, the rate equation is described as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{d}}{n}_{{\rm{A}}{{\rm{r}}}^{+}}}{{\rm{d}}t}=-{k}_{{\rm{c}}{\rm{h}}{\rm{e}}{\rm{x}}{\rm{c}}}{n}_{{\rm{A}}{{\rm{r}}}^{+}}{n}_{M}\\ \,-\,{{\rm{\Gamma }}}_{{\rm{A}}{{\rm{r}}}^{+}}\displaystyle \frac{{S}_{T}+({S}_{{\rm{I}}{\rm{R}}}-{S}_{T})(1-{\beta }_{{\rm{A}}{{\rm{r}}}^{+}})}{{V}_{{\rm{I}}{\rm{R}}}}\\ \,+\,({k}_{iz}{n}_{{\rm{e}}}+{k}_{iz}^{{\rm{H}}}{n}_{{{\rm{e}}}^{{\rm{H}}}})({n}_{{\rm{A}}{\rm{r}}}+{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{H}}}})\\ \,+\,({k}_{{\rm{m}}iz}{n}_{{\rm{e}}}+{k}_{{\rm{m}}iz}^{{\rm{H}}}{n}_{{{\rm{e}}}^{{\rm{H}}}}){n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}\end{array}\end{eqnarray} \tag{ A2 }$$

${{\rm{\Gamma }}}_{{\rm{A}}{{\rm{r}}}^{+}}=0.4{n}_{{\rm{A}}{{\rm{r}}}^{+}}{u}_{{\rm{b}}{\rm{o}}{\rm{h}}{\rm{m}},{\rm{A}}{{\rm{r}}}^{+}}$ is the flux of Ar⁺ leaving the ionization region. The Bohm speed of Ar⁺ is defined as ${u}_{{\rm{b}}{\rm{o}}{\rm{h}}{\rm{m}},{\rm{A}}{{\rm{r}}}^{+}}=\sqrt{\tfrac{e{T}_{{\rm{e}}}}{{{m}}_{{\rm{A}}{\rm{r}}}}}.$

For metastable atom Ar^m, the rate equation can be written as

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}\displaystyle \frac{{\rm{d}}{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}}{{\rm{d}}t}=-({k}_{{\rm{m}}{\rm{i}}{\rm{z}}}{n}_{{\rm{e}}}+{k}_{{\rm{m}}{\rm{i}}{\rm{z}}}^{{\rm{H}}}{n}_{{{\rm{e}}}^{{\rm{H}}}}){n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}-{k}_{p}{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}{n}_{{ \mathcal M }}\\ \,\,\,\,\,-\,({k}_{{\rm{d}}{\rm{e}}{\rm{x}}}{n}_{{\rm{e}}}+{k}_{{\rm{d}}{\rm{e}}{\rm{x}}}^{{\rm{H}}}{n}_{{{\rm{e}}}^{{\rm{H}}}}){n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}\\ \,\,\,\,\,-\,{{\rm{\Gamma }}}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}},{\rm{d}}{\rm{i}}{\rm{f}}{\rm{f}}}\displaystyle \frac{{S}_{{\rm{I}}{\rm{R}}}-{S}_{T}}{{V}_{{\rm{I}}{\rm{R}}}}-{{\rm{\Gamma }}}_{M,{\rm{c}}{\rm{o}}{\rm{l}}{\rm{l}}}\displaystyle \frac{{m}_{M}}{{m}_{{\rm{A}}{\rm{r}}}}\\ \,\times \,\displaystyle \frac{{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}}{{n}_{{\rm{A}}{\rm{r}}}+{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}}\displaystyle \frac{{S}_{{\rm{I}}{\rm{R}}}-{S}_{T}}{{V}_{{\rm{I}}{\rm{R}}}}+({k}_{{\rm{e}}{\rm{x}}}{n}_{{\rm{e}}}+{k}_{{\rm{e}}{\rm{x}}}^{{\rm{H}}}{n}_{{{\rm{e}}}^{{\rm{H}}}})({n}_{{\rm{A}}{\rm{r}}}+{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{H}}}}).\end{array}\end{array}\,\end{eqnarray} \tag{ A3 }$$

Considering the effect of heat on particle overflow, the diffusion flux of Ar^m is ${{\rm{\Gamma }}}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}},{\rm{d}}{\rm{i}}{\rm{f}}{\rm{f}}}=\tfrac{{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}{u}_{{T}_{{\rm{A}}{\rm{r}}},{\rm{IR}}}}{4}.$ ${n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}$ is the Ar^m density of the ionization region.

For Ar^H, the rate equation is

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{d}}{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{H}}}}}{{\rm{d}}t}=-({k}_{{\rm{i}}{\rm{z}}}{n}_{{\rm{e}}}+{k}_{{\rm{i}}{\rm{z}}}^{{\rm{H}}}{n}_{{{\rm{e}}}^{{\rm{H}}}}){n}_{{\rm{A}}{{\rm{r}}}^{{\rm{H}}}}\\ \,\,\,\,-\,({k}_{{\rm{e}}{\rm{x}}}{n}_{{\rm{e}}}+{k}_{{\rm{e}}{\rm{x}}}^{{\rm{H}}}{n}_{{{\rm{e}}}^{{\rm{H}}}}){n}_{{\rm{A}}{{\rm{r}}}^{{\rm{H}}}}\\ \,\,\,\,-\,{{\rm{\Gamma }}}_{{\rm{A}}{{\rm{r}}}^{{\rm{H}}},{\rm{d}}{\rm{i}}{\rm{f}}{\rm{f}}}\displaystyle \frac{{S}_{{\rm{I}}{\rm{R}}}-{S}_{T}}{{V}_{{\rm{I}}{\rm{R}}}}+{{\rm{\Gamma }}}_{{\rm{A}}{{\rm{r}}}^{+}}\displaystyle \frac{{S}_{T}}{{V}_{{\rm{I}}{\rm{R}}}}.\end{array}\end{eqnarray} \tag{ A4 }$$

The rate equation for metal atoms is

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{d}}{n}_{M}}{{\rm{d}}t}=-({k}_{{M}{\rm{i}}{\rm{z}}}{n}_{{\rm{e}}}+{k}_{{M}{\rm{i}}{\rm{z}}}^{{\rm{H}}}{n}_{{{\rm{e}}}^{{\rm{H}}}}){n}_{M}-{k}_{p}{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}{n}_{M}\\ \,\,\,\,\,\,\,-\,{k}_{{\rm{c}}{\rm{h}}{\rm{e}}{\rm{x}}{\rm{c}}}{n}_{{\rm{A}}{{\rm{r}}}^{+}}{n}_{M}-{{\rm{\Gamma }}}_{M,{\rm{d}}{\rm{i}}{\rm{f}}{\rm{f}}}\displaystyle \frac{{S}_{{\rm{I}}{\rm{R}}}-{S}_{T}}{{V}_{{\rm{I}}{\rm{R}}}}\\ \,\,\,+\,({{\rm{\Gamma }}}_{{\rm{A}}{{\rm{r}}}^{+}}{Y}_{{\rm{A}}{{\rm{r}}}^{+}}+{{\rm{\Gamma }}}_{{M}^{+}}{Y}_{{M}^{+}}+{{\rm{\Gamma }}}_{{M}^{2+}}{Y}_{{M}^{2+}})\displaystyle \frac{{S}_{T}}{{V}_{{\rm{I}}{\rm{R}}}},\end{array}\end{eqnarray} \tag{ A5 }$$

where ${{\rm{\Gamma }}}_{M,{\rm{d}}{\rm{i}}{\rm{f}}{\rm{f}}}$ represents the diffusion flux of metal atoms that have not collided. It can be expressed by ${{\rm{\Gamma }}}_{M,{\rm{d}}{\rm{i}}{\rm{f}}{\rm{f}}}={{\rm{\Gamma }}}_{M,0}(1-{F}_{{\rm{c}}{\rm{o}}{\rm{l}}{\rm{l}}}).$ ${{\rm{\Gamma }}}_{{M}^{+}}=0.4{n}_{{M}^{+}}{u}_{{\rm{b}}{\rm{o}}{\rm{h}}{\rm{m}},{M}^{+}}$ and ${{\rm{\Gamma }}}_{{M}^{2+}}\,=0.4{n}_{{M}^{2+}}{u}_{{\rm{b}}{\rm{o}}{\rm{h}}{\rm{m}},{M}^{2+}}$ are the fluxes of metal ions leaving the ionization region. ${u}_{{\rm{b}}{\rm{o}}{\rm{h}}{\rm{m}},{M}^{2+}}$ is $\sqrt{2}$ times as much as ${u}_{{\rm{b}}{\rm{o}}{\rm{h}}{\rm{m}},{M}^{+}}.$ ${Y}_{{\rm{A}}{{\rm{r}}}^{+}},$ ${Y}_{{M}^{+}}$ and ${Y}_{{M}^{2+}}$ are the sputtering yields of Ar⁺ and metal ions with target and calculate it in advance through the TRIM software. By fitting, the functional relationship between the discharge voltage and sputtering yield $Y=a{U}^{b}$ is obtained. For Cu ions, the fitting parameters are $a=0.0691$ and $b=0.556.$ For Ar⁺, the fitting parameters are $a=0.1421$ and $b=0.468.$ For Al ions, the fitting parameters are $a=0.0849$ and $b=0.408.$

The rate equations for M⁺ and M²⁺ can be described as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{d}}{n}_{{M}^{+}}}{{\rm{d}}t}=({k}_{M{\rm{i}}{\rm{z}}}{n}_{{\rm{e}}}+{k}_{M{\rm{i}}{\rm{z}}}^{H}{n}_{{{\rm{e}}}^{H}}){n}_{M}+{k}_{p}{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}{n}_{M}\\ \,\,\,+\,{k}_{{\rm{c}}{\rm{h}}{\rm{e}}{\rm{x}}{\rm{c}}}{n}_{{\rm{A}}{{\rm{r}}}^{+}}{n}_{M}-({k}_{M{\rm{i}}{\rm{z}}2}{n}_{{\rm{e}}}+{k}_{M{\rm{i}}{\rm{z}}2}^{H}{n}_{{{\rm{e}}}^{{\rm{H}}}}){n}_{{M}^{+}}\\ \,\,\,-\,{{\rm{\Gamma }}}_{{M}^{+}}\displaystyle \frac{{S}_{T}+({S}_{{\rm{I}}{\rm{R}}}-{S}_{T})(1-{\beta }_{{M}^{+}})}{{V}_{{\rm{I}}{\rm{R}}}},\end{array}\end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{d}}{n}_{{M}^{2+}}}{{\rm{d}}t}=({k}_{M{\rm{i}}{\rm{z}}2}{n}_{{\rm{e}}}+{k}_{M{\rm{i}}{\rm{z}}2}^{{\rm{H}}}{n}_{{{\rm{e}}}^{{\rm{H}}}}){n}_{{M}^{+}}\\ \,\,\,-\,{{\rm{\Gamma }}}_{{M}^{2+}}\displaystyle \frac{{S}_{T}+\left({S}_{{\rm{I}}{\rm{R}}}-{S}_{T}\right)\left(1-{\beta }_{{M}^{2+}}\right)}{{V}_{{\rm{I}}{\rm{R}}}}.\end{array}\end{eqnarray} \tag{ A7 }$$

A.2. Particle balance for electrons

The rate equation for the electrons is

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{d}}{n}_{{\rm{e}}}}{{\rm{d}}t}=({k}_{{\rm{i}}{\rm{z}}}{n}_{{\rm{e}}}+{k}_{{\rm{i}}{\rm{z}}}^{{\rm{H}}}{n}_{{{\rm{e}}}^{{\rm{H}}}})({n}_{{\rm{A}}{\rm{r}}}+{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{H}}}})\\ \,\,\,+\,({k}_{{\rm{m}}{\rm{i}}{\rm{z}}}{n}_{{\rm{e}}}+{k}_{{\rm{m}}{\rm{i}}{\rm{z}}}^{{\rm{H}}}{n}_{{{\rm{e}}}^{{\rm{H}}}}){n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}+({k}_{M{\rm{i}}{\rm{z}}}{n}_{{\rm{e}}}+{k}_{M{\rm{i}}{\rm{z}}}^{{\rm{H}}}{n}_{{{\rm{e}}}^{{\rm{H}}}}){n}_{M}\\ \,\,\,+\,({k}_{M{\rm{i}}{\rm{z}}2}{n}_{{\rm{e}}}+{k}_{M{\rm{i}}{\rm{z}}2}^{{\rm{H}}}{n}_{{{\rm{e}}}^{{\rm{H}}}}){n}_{{M}^{+}}+\,{k}_{p}{n}_{M}{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}-{{\rm{\Gamma }}}_{{\rm{e}}}\displaystyle \frac{{S}_{{\rm{I}}{\rm{R}}}-{S}_{T}}{{V}_{{\rm{I}}{\rm{R}}}},\end{array}\end{eqnarray} \tag{ A8 }$$

where the electron flux can be expressed as ${{\rm{\Gamma }}}_{{\rm{e}}}\approx 0.4{n}_{{\rm{e}}}\tfrac{{D}_{{\rm{e}}}}{{{\rm{R}}}_{{\rm{I}}{\rm{R}}}}\left(1+\tfrac{{U}_{{\rm{I}}{\rm{R}}}}{{T}_{{\rm{e}}}}\right).$ ${D}_{{\rm{e}}}=\tfrac{1}{{\omega }_{{\rm{g}}{\rm{e}}}{\tau }_{{\rm{c}}}}\tfrac{{T}_{{\rm{e}}}}{B}$ is the Bohm diffusion coefficient of electrons and B is the average magnetic field at the boundary of the cathode sheath. A simplified Poisson equation is used to calculate the potential drop U_IR in the ionization region

$$\begin{eqnarray}&&{U}_{{\rm{I}}{\rm{R}}}=-\displaystyle \frac{e{R}_{{\rm{I}}{\rm{R}}}^{2}}{{\varepsilon }_{0}}({n}_{{\rm{A}}{{\rm{r}}}^{+}}+{n}_{{M}^{+}}+2{n}_{{M}^{2+}}-{n}_{{\rm{e}}}-{n}_{{{\rm{e}}}^{{\rm{H}}}}).\end{eqnarray} \tag{ A9 }$$

Introducing the hot electron diffusion effect, the rate equation of the hot electrons is as follows

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{n}_{{{\rm{e}}}^{{\rm{H}}}}}{{\rm{d}}t}=\displaystyle \frac{1}{e{U}_{{\rm{D}}}/2}\left(\displaystyle \frac{{U}_{S{\rm{H}}}{I}_{{\rm{s}}{\rm{e}}}}{{V}_{{\rm{I}}{\rm{R}}}}-{Q}^{{\rm{H}}}\right)-{{\rm{\Gamma }}}_{{e}^{{\rm{H}}}}\displaystyle \frac{{S}_{{\rm{I}}{\rm{R}}}-{S}_{T}}{{V}_{{\rm{I}}{\rm{R}}}}\end{eqnarray} \tag{ A10 }$$

${I}_{{\rm{s}}{\rm{e}}}=e{S}_{T}({\gamma }_{{\rm{A}}{{\rm{r}}}^{+}}{{\rm{\Gamma }}}_{{\rm{A}}{{\rm{r}}}^{+}}+{\gamma }_{{M}^{+}}{{\rm{\Gamma }}}_{{M}^{+}}+{\gamma }_{{M}^{2+}}{{\rm{\Gamma }}}_{{M}^{2+}})$ is the secondary electron current, U_D is the pulse voltage and U_SH = U_D−U_IR is the sheath voltage. $\gamma$ is the secondary electron emission coefficient. In this model, ${\gamma }_{{\rm{A}}{{\rm{r}}}^{+}}=0.112,$ ${\gamma }_{{M}^{+}}=0$ and ${\gamma }_{{M}^{2+}}=0.082\,16.$ Q^H represent the energy loss in inelastic collisions between hot electrons and atoms [24]

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{Q}^{{\rm{H}}}}{e}=({E}_{{\rm{A}}{\rm{r}},{\rm{c}}}+{E}_{{\rm{h}}{\rm{t}}{\rm{c}}}){k}_{{\rm{i}}{\rm{z}}}^{{\rm{H}}}{n}_{{{\rm{e}}}^{{\rm{H}}}}({n}_{{\rm{A}}{\rm{r}}}+{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{H}}}})\\ \,\,\,\,+\,({E}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}},{\rm{c}}}+{E}_{{\rm{h}}{\rm{t}}{\rm{c}}}){k}_{{m}{\rm{i}}{\rm{z}}}^{{\rm{H}}}{n}_{{{\rm{e}}}^{{\rm{H}}}}{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}\\ \,\,\,\,+\,({E}_{M,{\rm{c}}}+{E}_{{\rm{h}}{\rm{t}}{\rm{c}}}){k}_{M{\rm{i}}{\rm{z}}}^{{\rm{H}}}{n}_{{e}^{{\rm{H}}}}{n}_{M}+({E}_{{M}^{+},{\rm{c}}}+{E}_{{\rm{h}}{\rm{t}}{\rm{c}}})\\ \,\,\,\,\times \,{k}_{M{\rm{i}}{\rm{z}}2}^{{\rm{H}}}{n}_{{e}^{{\rm{H}}}}{n}_{{M}^{+}}-{E}_{{\rm{d}}{\rm{e}}{\rm{x}}}{k}_{{\rm{d}}{\rm{e}}{\rm{x}}}^{{\rm{H}}}{n}_{{{\rm{e}}}^{{\rm{H}}}}{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}\end{array}\end{eqnarray} \tag{ A11 }$$

A.3. Energy balance

T_e is obtained by solving the energy balance equation, as equation (A11)

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\left(\displaystyle \frac{3}{2}e{n}_{{\rm{e}}}{T}_{{\rm{e}}}\right)={F}_{{\rm{P}}{\rm{W}}{\rm{R}}}\displaystyle \frac{{P}_{{\rm{D}}}}{{V}_{{\rm{I}}{\rm{R}}}}-Q-(\displaystyle \frac{3}{2}e{T}_{{\rm{e}}}{{\rm{\Gamma }}}_{{\rm{e}}}\displaystyle \frac{{S}_{{\rm{I}}{\rm{R}}}-{S}_{T}}{{V}_{{\rm{I}}{\rm{R}}}}\\ \,+\,\displaystyle \sum _{i={\rm{A}}{{\rm{r}}}^{+},{M}^{+},{M}^{2+}}\displaystyle \frac{1}{2}e(q{T}_{{\rm{e}}}){{\rm{\Gamma }}}_{i}\displaystyle \frac{{S}_{T}+({S}_{{\rm{I}}{\rm{R}}}-{S}_{T})(1-{\beta }_{i})}{{V}_{{\rm{I}}{\rm{R}}}})\end{array}\end{eqnarray} \tag{ A12 }$$

P_D = U_D × I_D is the input power of the DC discharge, and I_D and U_D are obtained by experiments and volt–current characteristics. Q includes the energy losses by electron impact ionization with neutral atoms (Ar, Ar^m, Ar^H, and metal atoms) and energy losses of Penning ionization, plus the energy to de-excitation release as shown in equation (A13)

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{Q}{e}={E}_{{\rm{A}}{\rm{r}},{\rm{c}}}{k}_{{\rm{i}}{\rm{z}}}{n}_{{\rm{e}}}({n}_{{\rm{A}}{\rm{r}}}+{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{H}}}})+{E}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}},{\rm{c}}}{k}_{{\rm{m}}{\rm{i}}{\rm{z}}}{n}_{{\rm{e}}}{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}\\ \,\,+\,{E}_{M,{\rm{c}}}{k}_{M{\rm{i}}{\rm{z}}}{n}_{{\rm{e}}}{n}_{M}+{E}_{{M}^{+},{\rm{c}}}{k}_{M{\rm{i}}{\rm{z}}2}{n}_{{\rm{e}}}{n}_{{M}^{+}}\\ \,\,-\,{E}_{{\rm{d}}{\rm{e}}{\rm{x}}}{k}_{{\rm{d}}{\rm{e}}{\rm{x}}}{n}_{{\rm{e}}}{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}-{E}_{p}{k}_{p}{n}_{M}{n}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}}}.\end{array}\end{eqnarray} \tag{ A13 }$$

Zoom In Zoom Out Reset image size

The collisional energy loss for each electron–ion pair are described as equation (A14)

$$\begin{eqnarray}\begin{array}{ccl}{k}_{{\rm{i}}{\rm{z}}}{E}_{{\rm{A}}{\rm{r}},c} & = & {k}_{{\rm{i}}{\rm{z}}}{E}_{{\rm{i}}{\rm{z}}}+{k}_{{\rm{e}}{\rm{x}}}{E}_{{\rm{e}}{\rm{x}}}+{k}_{{\rm{e}}{\rm{l}}}\displaystyle \frac{3{{\rm{m}}}_{{\rm{e}}}}{{{\rm{m}}}_{{\rm{A}}{\rm{r}}}}{T}_{{\rm{e}}}\\ {k}_{{\rm{m}}{\rm{i}}{\rm{z}}}{E}_{{\rm{A}}{{\rm{r}}}^{{\rm{m}}},c} & = & {k}_{{\rm{m}}{\rm{i}}{\rm{z}}}{E}_{{\rm{i}}{\rm{z}}}+{k}_{{\rm{e}}{\rm{l}}}\displaystyle \frac{3{{\rm{m}}}_{{\rm{e}}}}{{{\rm{m}}}_{{\rm{A}}{\rm{r}}}}{T}_{{\rm{e}}}\\ {k}_{M{\rm{i}}{\rm{z}}}{E}_{M,c} & = & {k}_{M{\rm{i}}{\rm{z}}}{E}_{M{\rm{i}}{\rm{z}}}+{k}_{M{\rm{e}}{\rm{l}}}\displaystyle \frac{3{{\rm{m}}}_{{\rm{e}}}}{{{\rm{m}}}_{M}}{T}_{{\rm{e}}}\\ {k}_{M{\rm{i}}{\rm{z}}2}{E}_{{M}^{+},c} & = & {k}_{M{\rm{i}}{\rm{z}}2}{E}_{M{\rm{i}}{\rm{z}}2}+{k}_{M{\rm{e}}{\rm{l}}}\displaystyle \frac{3{{\rm{m}}}_{{\rm{e}}}}{{{\rm{m}}}_{{M}^{+}}}{T}_{{\rm{e}}}.\end{array}\end{eqnarray} \tag{ A14 }$$

The ion return probability β can be expressed as (A15)

$$\begin{eqnarray}\begin{array}{cc}{\beta }_{i}=\left\{\begin{array}{cc}0 & {U}_{{\rm{I}}{\rm{R}}}\leqslant {T}_{i}\\ 1-\tfrac{{T}_{i}}{q{U}_{{\rm{I}}{\rm{R}}}} & {U}_{{\rm{I}}{\rm{R}}}\gt {T}_{i}\end{array}\right. & i={\rm{A}}{{\rm{r}}}^{+},{M}^{+},{M}^{2+}\end{array}\end{eqnarray} \tag{ A15 }$$

and reveals the determined method of the free parameter ${\omega }_{{\rm{ge}}}{\tau }_{{\rm{c}}},$ where I_calc is the simulated discharge current used to fit the experimental current I_D

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}{I}_{{\rm{c}}{\rm{a}}{\rm{l}}{\rm{c}}}=e{{\rm{\Gamma }}}_{{\rm{A}}{{\rm{r}}}^{+}}{S}_{T}(1+{\gamma }_{{{ \mathcal M }}^{+}})+e{{\rm{\Gamma }}}_{{{M}}^{2+}}{S}_{T}(2+{\gamma }_{{ \mathcal M }{}^{2+}})\\ \,+\,e{{\rm{\Gamma }}}_{{ \mathcal M }{}^{+}}{S}_{T}(1+{\gamma }_{{ \mathcal M }{}^{+}}).\end{array}\end{array}\end{eqnarray} \tag{ A16 }$$