Influence of ion-induced secondary electron emission parameters in PICMC plasma simulations with experimental validation in DC cylindrical diode and magnetron discharges

The experimental system (figure 1) designed to validate the simulation results consists of a 1.1 m vertical vacuum chamber of 200 mm inner diameter. Pumping is ensured by a turbomolecular pump (60 l/s pumping speed for N₂) backed up by a primary pump. The system base pressure is in the 10⁻⁷–10⁻⁶ mbar range in the main chamber without bake-out. Argon is injected in the main chamber through a variable leak valve, and two capacitive diaphragm gauges are located at the top and bottom of the chamber in order to precisely control the argon process pressure in the chamber during the plasma discharge. A coil made of 2.8 mm diameter insulated copper cable and wound around the main chamber provides a solenoidal axial magnetic field described in section 2.4. This solenoid has an inner diameter of 220 mm, and a total height of 990 mm.

Zoom In Zoom Out Reset image size

A niobium rod cathode (18 mm diameter, 303 mm length) is coaxially mounted inside a stainless steel anode (125 mm average diameter, 500 mm length) as shown in figure 1. The anode is made of two half shells creating an octagon, such that diagnostics (quartz balance, sample for thickness analysis not used in this study) can be easily mounted on the planar faces while keeping the azimuthal symmetry close to that of a cylinder. The cathode–anode distance of 51.5 mm is chosen equal to the one between substrate and cathode in the HIE-ISOLDE coating assembly [22]. The anode is grounded through fixation threads attached to the main chamber top flange, and a negative bias voltage is applied on the cathode by an MDX 500 advanced energy power supply unit. This unit can deliver up to 1000 V–0.5 A with an accuracy of 0.2% of the full rated output, and is used to apply and read discharge voltage and current values. It is operated in power regulation mode.

A Langmuir probe is mounted on a linear motion vacuum feedthrough of 405 mm stroke length. This allows measurements of axially resolved plasma density, electron temperature and plasma potential to compare with simulation results. The probe consists of a 1 mm diameter tungsten wire (r_probe = 0.5 mm) inserted in a 1.6 mm inner diameter, 3 mm outer diameter alumina tube attached to a vertical rod outside of the anode and connected to a vacuum-compatible coaxial cable insulated with kapton, as seen in figures 1 and 2. The effective tungsten tip length L_probe immersed in the plasma is of 2 mm, and the middle of the tip is located at 27 mm from the cathode axis with an uncertainty of ±1.5 mm due to the geometrical design. A slit in the anode allows vertical motion of the probe.

Zoom In Zoom Out Reset image size

To obtain values of plasma parameters comparable to the simulated ones, the probe is operated in sweeping mode by applying a triangular voltage V_a varying between −40 V and +3 V at a sweeping frequency f = 13.1 Hz. The current drawn by the probe is recorded by a Lecroy HDO6104 oscilloscope. Current density/voltage (J/V_a) characteristics are analysed according to the procedure suggested in [23, appendix B.4]. First, the plasma potential V₀ [V] is estimated as the voltage corresponding to the maximum of the J/V_a characteristics first derivative. Then, the electron temperature T_e [eV] is estimated as ${T}_{\mathrm{e}}=\frac{{J}_{0}}{{k}_{0}}$ ⁴ where J₀ [A m⁻²] and k₀ [A m⁻² V⁻¹] are, respectively, the current density and the slope of the I/V characteristics at the plasma potential. The plasma density n₀ [m⁻³] is estimated as ${n}_{0}=\frac{{e}^{1/2}{J}_{0}}{q\sqrt{\frac{q{T}_{\mathrm{e}}}{{m}_{\mathrm{e}}}}}$ where e is the base of the natural logarithm, and V_f [V] is the floating potential. These four values are taken as initial guesses in a four-parameter fit of the J/V_a characteristic from V_f to V₀, with the approximated expression of the current density collected by a cylindrical probe J_approx(V_a) [A m⁻²] given by Guittienne et al [23] in equation (B12):

$$\begin{equation}{J}_{\mathrm{a}\mathrm{p}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{x}}\left({V}_{\mathrm{a}}\right)=q{u}_{\mathrm{t}\mathrm{h}\mathrm{e}}{n}_{0}{e}^{-\frac{1}{2}-\frac{{V}_{0}}{{T}_{\mathrm{e}}}}\left[{e}^{\frac{{V}_{\mathrm{a}}}{{T}_{\mathrm{e}}}}-{e}^{\frac{{V}_{\mathrm{f}}}{{T}_{\mathrm{e}}}}\right].\end{equation} \tag{ 1 }$$

This fit provides final values for T_e, V₀, n₀ and V_f, with experimental error bars taking into account probe area uncertainty and error estimates on the determination of initial V₀, n₀ and E_e (converted from T_e by assuming a Maxwellian electron energy distribution) before the final four-parameter fit. At last, we estimate that the methodology used here to analyse (J/V_a) characteristics is associated with an error margin of ∼20% in the pressure range of interest.

Plasma simulations are performed on a high performance computing cluster installed at CERN with a PICMC parallel code developed at the Fraunhofer IST [24–26]. The simulation geometry includes cathode, anode, insulating ceramic and surrounding vacuum chamber. Surface meshing and physical surfaces definition are obtained with the open source software Gmsh [27], which is also used for visualising simulation results. Anode and chamber are defined as grounded elements, insulating ceramic as floating, and the negative bias voltage on the cathode is adjusted throughout the simulation run by a feedback loop enforcing the power setpoint.

A summary of the numerical and physical parameters is presented in table 1. To compare simulation results with experiments, a power setpoint P_dis of 20 W is applied for all cases. Higher powers would result in higher currents and plasma densities, and unfulfillable numerical constraints in terms of time step and mesh size due to computational costs. The Cartesian simulation volume mesh is refined inside the anode where the plasma is denser, and made coarser outside of it. To ensure proper collision statistics, the scaling factor representing the number of argon simulation particles with respect to real particles is adjusted for each pressure to have about ten argon numerical particles in each 1 mm³ cell inside the anode. Volume reactions used in the simulation model involving argon atoms, ions, and electrons were described in [26], while surface reactions assume that electrons and argon ions impinging on surfaces are collected there, with the emission of one neutral argon atom and a number of electrons according to the value of γ_IISEE for each collected ion. Electrons impinging on a physical surface are collected with a probability of 1 without further re-emission. Secondary electron emission under neutral, metastable or photon bombardment is also neglected.

The choice of these numerical parameters could lead to some inaccuracies in the simulation results presented in this manuscript. First, the mesh size of 1 mm is not resolving the Debye length for the typical densities and energies simulated here. This could induce so-called numerical heating [28], meaning that insufficient spatial resolution of the electric field could result in artificially too high electron kinetic energies. However, a finer mesh size would make simulation runs even longer and practically impossible to carry. To partly mitigate this issue, the stability of each simulation in terms of cathode voltage convergence, sufficient cathode sheath resolution and absence of abnormally high electron energies in the plasma volume are ensured. Then, absence of secondary electron emission due to sources other than ion bombardment and uncertainties inherent to the collision cross-sections model could also alter the presented results. The impact of these choices on the simulation inaccuracies and in particular on the assessed IISEE parameters is quite challenging to quantify and is beyond the scope of the present study. Discussions on this matter can be found in [10, 12, 13] in the case of radio-frequency discharges simulated with one-dimensional in space codes.

Magnetron simulations require a magnetic field map as an input to the plasma simulation code. This map can either be computed by the boundary element method module of the PICMC code in the case of permanent magnet assemblies, or imported as a 3D text file from an external source. In the present case, we choose the latter option and compute the magnetic field generated by the solenoid using the magnetostatic module of the Opera commercial simulation software [29]. Figure 3 shows the good agreement between the simulated axial component of the magnetic field and the corresponding quantity measured on the solenoid axis with a gaussmeter for different coil currents. The simulated axial component of the magnetic field corresponding to a simulation current density of j = 0.581 A mm⁻² is found to match the experimental value generated by a current I_coil = 10 A, as shown by the blue curves in figure 3. A linear scaling of the simulation current density and of the experimental I_coil at, respectively, j = 2.324 A mm⁻² and I_coil = 40 A shows a consistent agreement between simulation and measurement (red curves of figure 3). Therefore, the simulated field generated with j = 4.648 A mm⁻² (whose axial component is plotted with a green straight line in figure 3) is selected for the magnetron simulations which will be presented in section 3.2 as an accurate 3D representation of the experimental field generated with a current I_coil = 80 A (not measured). The axial positioning of the anode at 200–700 mm below the main chamber top flange is chosen such that the anode lies in the uniform field region (B ∼ 160 Gauss, see figure 3), thus minimising edge effects on the plasma due to non-uniformity of the magnetic field.

Zoom In Zoom Out Reset image size

In this section, a suitable value of γ_IISEE is numerically assessed in the diode configuration by comparing global discharge currents and voltages between simulations and experiments as a function of the argon gas pressure.

Three simulations are performed in the DC diode configuration for different γ_IISEE with p_Ar = 0.3 mbar. As can be seen in table 2, increasing γ_IISEE leads to a discharge voltage decrease (respectively a current increase) since it entails more electrons emitted from the cathode towards the plasma region. Moreover, a simulation yield of 0.13 fits the experimental current and voltage values for this given pressure within the experimental uncertainties.

Then, we perform a parameter scan in process pressure p_Ar while keeping γ_IISEE = 0.13. In diode coating processes such as the one of the HIE-ISOLDE cavity [22], a pressure increase leads to a large number of scattering collisions between sputtered niobium and argon atoms and hence high niobium redeposition on the cathode and very low deposition rates [30]. Even though we do not consider sputtering processes in the present study, we choose to vary p_Ar between 0.2 and 0.4 mbar to be relevant with typical diode process pressures.

It is shown in table 3 that the simulation I/V trend and absolute values match the experimental data. Validity of the IISEEY value is further confirmed by a loss of plasma sustainability for both simulation and experiment at p_Ar = 0.05 mbar, with voltages saturating at 1000 V and weak currents of, respectively, 1.2 mA and 2 mA being delivered, both failing to reach the 20 W power setpoint.

Physical parameters can be extracted from each simulation case, such as electron and ion densities, electric potential and species energies. A map of the electron density distribution is displayed in figure 4(a) for p_Ar = 0.3 mbar, along with a photograph (visible light spectrum) of the experimental diode plasma taken through a viewport located below the vacuum chamber (figure 4(b)). Figure 5 presents radial and axial plots of the electron density and energy, and radial plot of the electric potential for the three pressure cases. Radial plots represent azimuthally averaged values in the horizontal mid-plane of the cathode, in the discharge gap (9–60.5 mm). Axial plots are similarly obtained from azimuthally averaged values over a 27 mm ± 1 mm radius cylinder around cathode axis, with its axis origin at mid height of the cathode. Displayed electron energies correspond to mean electron energies within each cell.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Increasing the pressure leads to a radial confinement of the plasma as shown by the radial density peak position shift (figure 5(a)), along with a cathode sheath width reduction (figure 5(e)). Indeed, a higher pressure means that electrons accelerated through the cathode sheath experience more collisions with neutral argon particles, hence losing their energy closer to the cathode (figure 5(b)). As seen in figure 5(c), the axial electron density profiles are almost flat around the cathode and quickly drop at the cathode edges (±151.5 mm). This confinement can be explained by two reasons. First, in the absence of a magnetic field, the plasma is electrically driven. Due to the aspect ratio between cathode length (303 mm) and anode average diameter (120 mm), the electric field is almost constant and radially directed, which ensures the axial plasma uniformity around the cathode. Second, the relatively high process pressures limit the axial drift of charged species because of collisions, thus leading to the abrupt electron density axial drops at the cathode edges.

Langmuir probe measurements are performed at the three pressures along the vertical direction, according to the probe configuration described in section 2.2. The probe surface used to convert the probe current into current density is the geometrical area of the probe tip A_p = $2\pi {r}_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{e}}{L}_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{e}}+\pi {r}_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{e}}^{2}$, where r_probe = 0.5 mm and L_probe = 2 mm. Similarly to the simulation results of figure 5, experimental plasma parameters taken at 27 mm radius from the cathode axis are almost pressure independent in the studied pressure range. Therefore, we only present in figure 6 the comparison of simulated and experimental plasma density, electron energy and plasma potential at p_Ar = 0.3 mbar. We use a synthetic diagnostic to compare simulations and experiments, where simulated values are extracted from a cylinder of radius 27 mm ± 2.5 mm such that they match the probe spatial resolution. As such, averaged simulated results are presented with error bars corresponding to ±3 standard deviations over the 5 mm radius averaging. Analysis of the J/V_a probe characteristics yields the electron temperature T_e [eV], which is converted into electron energy E_e [eV] for comparison with the simulated profile of figure 6(b). This conversion assumes a Maxwellian energy distribution for the electrons, i.e. E_e [eV] = $\frac{3}{2}{T}_{\mathrm{e}}$ [eV].

Zoom In Zoom Out Reset image size

Figure 6 shows that experimental electron density (electron energy and plasma potential, respectively) is about twice (respectively about half) of its simulated equivalent in the central axial region (z ∈ [−150; 150] mm). The discrepancy between simulation and experiment could be partly explained by plasma perturbations induced when the probe voltage V_a is positively swept above the floating potential V_f. Indeed, electron peak currents collected by the probe are of ∼6 mA, compared to a global discharge current of ∼51 mA (see table 3 for p_Ar = 0.3 mbar). This could mean that for positive probe voltages V_a, electrons are collected from a plasma region beyond the probe vicinity. For a given positive voltage V_a above V_f, the electron current is therefore overestimated. This results in an overestimate of the plasma density n_e along with an underestimate of the plasma potential V₀ and electron temperature T_e. This could also explain why the plasma density decay observed in the simulation profile for z < −150 mm and z > 150 mm is not well reproduced in the measurements, since the probe collects electrons from a large plasma region independently of its position. Real plasma parameters without probe perturbation should be closer to the simulated ones. It is also worth mentioning that the probe perturbation does not impact the values given in table 3, since those are taken when the probe is electrically floating and axially positioned at the anode extremity.

Electron trajectories are modified in the presence of a magnetic field. In particular and contrarily to the diode configuration, secondary electrons emitted from the cathode surface can be recollected depending on their initial energy distribution, which therefore plays a role in the global I/V discharge parameters. To validate the use of a Gaussian energy distribution of secondary electrons emitted from the cathode, magnetron simulations are run at different pressures with γ_IISEE = 0.13 according to the results of the previous section 3.1. The magnetic field profile used as input for the simulations was described in section 2.4, corresponding to an experimental axial field of 160 Gauss peak value. Four pressures are chosen, as given in table 1, as a trade-off between minimum pressure needed to sustain the plasma discharge and process-limited higher pressure (same criteria as for the diode study). The other physical and numerical parameters are unchanged.

Table 4 shows the experimental and simulated discharge currents and voltages for the four different pressures. The simulation trend matches the increase of current with pressure and corresponding decrease of voltage. Discrepancies of ∼25–35  V (∼10% of the discharge voltage) can be noticed between measurement and simulation for the same pressure. This could be partly attributed to longer simulation stabilization times in the magnetron configuration with respect to the diode case. Indeed, voltage convergence exhibits an exponential decay behaviour in time such that simulation values given at 30 μs are slightly overestimating voltages and thus underestimating currents. This slow convergence of the magnetron cases could be explained by the drift of the charged species outside of the cathode–anode region. This can be seen in figure 7(a) showing electron density losses at the anode ends, which lead to a larger effective plasma volume and longer computation times as well as longer physical time needed to reach a steady-state. This phenomenon of electron end-losses was described in [31] for such a configuration, and is usually avoided in coating systems such as planar or post-magnetrons. As a trade-off between reasonable computational time and required physical accuracy, the achieved agreement between simulated and experimental values meets the precision required for our practical purpose.

Zoom In Zoom Out Reset image size

Figure 8 presents the electron density and energy radial and axial profiles, and electric potential radial profile for the four pressure cases. Pressure-dependent variations of plotted quantities are less pronounced than for the diode cases (figure 5). The decay of plasma characteristics (n_e, T_e) is smoother along the z-axis because of electron end-losses bringing the plasma outside of the axial cathode region. Electron energies are about twice higher than in the diode configuration due to lower working pressures. The apparent saturation of the electron density axial profile when the pressure is increased (figure 8(c)) is attributed to a mitigation of the magnetic field influence when more collisions between charged particles and process gas atoms occur, thus limiting the transport of charged particles and tending towards the diode case.

Zoom In Zoom Out Reset image size

Langmuir probe measurements are taken at the four pressures in the vertical axis direction, similarly to the ones described in section 3.1. For the magnetron discharge, the probe surface used to convert the probe current into current density is taken as the projection of the probe tip geometrical surface onto the plane perpendicular to the magnetic field [32]: A_p = 2r_probe L_probe. Indeed, with B = 160 Gauss and E_e ∼ 2 eV according to simulated values of figure 8(d), the electron gyroradius r_L ∼ 0.3 mm is smaller than the probe typical size. Experimental plasma parameters taken at 27 mm radius from the cathode axis are almost identical for p_Ar ∈ [0.05; 0.07; 0.09] mbar. As such, we only compare experimental measurements with simulated parameters using the synthetic diagnostic described in section 3.1 for the two extreme pressures p_Ar = 0.03 mbar and 0.09 mbar, as shown in figure 9. Due to the influence of the magnetic field, a non-uniformity of the azimuthal plasma densities, as seen in figure 7(a), is leading to local fluctuations of the azimuthal electric potential profiles. Because of the very different time scales between simulations and Langmuir probe measurements, comparison of azimuthal electric potential profiles cannot be made, as opposed to what was showed in the diode case in figure 6(c).

Zoom In Zoom Out Reset image size

Figure 9(a) shows that simulated and experimental electron density profiles agree in shape except for z < −100 mm. The experimental increase in density for these axial positions is consistently measured for all pressures, but remains unexplained. The simulated densities in the central axial region are 6.5 and 3.5 times higher than the experimental ones, respectively, for p_Ar = 0.03 mbar and p_Ar = 0.09 mbar, but the trend with pressure matches: a larger pressure results in a larger electron density in both simulation and experiment. Figure 9(b) also reveals an agreement of relative electron energy profiles between simulation and experiment with a matching trend as a function of pressure. However, simulated electron energies are 2.6 to 1.8 times lower than their experimental counterparts, respectively, for p_Ar = 0.03 mbar and p_Ar = 0.09 mbar. The currents drawn by the probe are ten to twenty times smaller than the ones in diode and the plasma is therefore not disturbed in the vicinity of the probe. This could explain the reasonable agreement of measurements with simulations in magnetron despite the fact that the theory used for the analysis of Langmuir probe characteristics J/V_a does not account for the presence of the magnetic field.