Data processing techniques for ion and electron-energy distribution functions

Distribution functions offer a detailed description of plasma charged particles: from obtaining some useful plasma parameters (e.g., density, temperature, and potential) to deducing heating and transport mechanisms, collision and reaction rates as well as detecting the presence of particles sub-populations.^1–4 The ion-energy distribution function (IEDF) is extensively employed to characterize the ion kinetics in plasma processing techniques and in the development of spacecraft electric propulsion technologies, where the detection of accelerated ion populations, such as ion beams, is important.^2,5,6 The electron-energy distribution function (EEDF) gives an extensive picture of the dynamics of electrons. This is essential for characterizing the plasma generation mechanisms in laboratory and commercial plasma devices, e.g., to determine reaction rates⁷ as well as in investigating fundamental aspects of plasma physics, e.g., the thermodynamics of magnetized electrons.^8,9

The ion and energy distribution functions can be obtained from the first and second derivatives of the current–voltage (I-V) characteristics measured with appropriate diagnostics. Electrostatic probes, such as retarding field energy analyzers (RFEA) and single Langmuir probes (LP), are some of the most common instruments for measuring ion and electron distribution functions, respectively. The plasma potential, and the ion beam potential and density are usually obtained through the first derivative of an RFEA I-V characteristic. From a Langmuir probe I-V characteristic, the floating potential, electron temperature, plasma density, and plasma potential can be inferred via the classical Langmuir method.^10,11 While popular, thanks to its apparent simplicity, this method is only valid for plasmas with Maxwellian electron distributions and can be inaccurate in determining the plasma potential and ion current.^12,13 Alternatively, the Druyvesteyn method determines the electron distribution function from the second derivative of the LP I-V characteristic, a process inherently more robust and information-rich than the Langmuir method, yet more challenging owing to the double differentiation.¹⁴ 

A variety of analog and numerical differentiation routines have been used to evaluate IEDFs and EEDFs. Since measurement noise amplification is inherent to the differentiation process, filtering/fitting methods may be necessary to smooth the measured I-V data without causing significant distortion and to obtain accurate first and second derivatives. The aim of this paper is to review different signal processing and differentiation methods for the acquisition of optimal ion- and electron-energy distribution functions. The study discusses analog and numerical processing techniques that are routinely used to obtain EEDFs and IEDFs. With respect to the latter, the focus is on processing IEDFs that show the presence of multiple ion populations. The different numerical processes reviewed to obtain the first and second derivatives from the I-V characteristics are Savitzky–Golay (SG) filter, B-spline piecewise polynomial fitting (BS), Gaussian filter (GF), the Blackman window filter (BW), and the AC (alternating current) superimposition. The ion-energy distribution function is also computed using the Gaussian deconvolution method. The analysis presented is intended as a baseline approach for the manipulation of I-V characteristics obtained from experimental data in a plasma discharge. Specifically, data collected in a magnetized radio frequency plasma discharge at low gas pressures are used, which contain a higher noise level than, e.g., DC plasma discharges. This work gathers the lessons learned in the acquisition and processing of electric-probe data to ease the reader's process of obtaining charged-particles energy distribution functions of suitable quality.

The work is organized as follows. Sec. II describes the main theory and concepts behind ion and energy distribution functions as well as discusses the relevant state-of-the-art work. Section III presents a review of the most commonly employed analog and numerical filtering/differentiation methods for the evaluation of IEDFs and EEDFs in experimental plasmas. Section IV describes the experimental setup, including the plasma diagnostic probes. Finally, Sec. V provides an example application of data processing analysis for ion- and electron-energy distribution functions obtained from the different methods starting from raw experimental data.

IEDF and EEDF are usually obtained through a variety of analog and numerical differentiation routines. The implementation of filtering and/or fitting methods is often required since the significant noise amplification introduced by differentiating the measured I-V characteristic could result in inaccurate energy distributions. Among the numerical methods, the Savitzky–Golay filter is one of the most commonly employed for IEDF/EEDF derivations.^22–26 Bétchu et al.²⁴ obtained the second derivative of the electron current measured in a multi-dipolar microwave plasma source using the Savitzky–Golay filter with second order polynomials. The latter was also used by Gulbrandsen and Fredriksen²⁷ to smooth the I-V curve acquired with an RFEA and obtain the ion-energy distribution function. Fernández Palop et al.²⁸ were the first to implement a Gaussian filter to smooth the I-V characteristics measured in an argon DC discharge and evaluate the electron-energy distribution function tested for different pressures and discharge currents. This method was compared with the Savitzky–Golay filter and the B-spline approximation methods for a simulated dataset and noise level. In their study, the Gaussian filter performed better than the other methods (see Figs. 5–10) in Ref. 28. Magnus and Gudmundsson²⁹ compared the use of the Savitzky–Golay filter as well as the Gaussian filter, the Blackman window, and polynomial fitting, to smooth both simulated and experimental I-V curves. The mean-squared error, the correlation coefficient, and the residuals of the different methods were compared, coupled with visual assessments of the EEPF. Their study proved that the accuracy of the EEDF is affected by the chosen smoothing method. It was observed that the Savitzky–Golay filter did perform well for noisy data, while the polynomial fitting technique was deemed inadequate to process measured I-V characteristics. The Gaussian and Blackman window filters gave good results for both simulated and experimental data; however, the effective electron temperature was overestimated due to the flattening of the EEPF maximum.

Roh et al.³⁰ studied the use of both the Savitzky–Golay and the Blackman filter to obtain EEPFs from Langmuir probes measurements taken in a low-pressure plasma. They developed a method that used two Langmuir probes with different collecting areas coupled with the two filtering techniques to accurately resolve both the low electron energies and the inelastic energy regimes.

An AC superimposition technique was first proposed in 1934 by Sloane and MacGregor³¹ as a new way to acquire a signal proportional to the first and second derivatives of Langmuir probe characteristics. Different studies applied this technique to derive EEDFs either analogically^31–34 or numerically.³⁵ The analog AC superimposition method was also used to obtain IEDFs by Hatakeyama et al.³⁶ Bang and Chung³⁷ improved the numerical AC superimposition method of Jauberteau and Jauberteau³⁵ by using AC signals whose amplitudes increase with the electron energy as well as by correcting for distortions induced from high amplitude perturbations by taking into account higher order derivative terms (see Figs. 1 and 5 in Ref. 37).

Alternatively, an analog differentiation technique described by Schoenberg³⁸ has been often used to obtain EEPFs using appropriate electronics circuitry.^3,39–42 Takahashi et al.⁴⁰ evaluated electron-energy distribution functions using analog differentiation to estimate the temperature of trapped and free electrons in a low-pressure RF plasma. The derivatives of the measured current were obtained via two resonant circuits tuned to 800 Hz. The same technique was used by Takahashi et al.¹⁸ to compute the ion-energy distribution function of an accelerated ion population.

The IEDF has also been obtained using a Gaussian deconvolution method.^2,17,20 For instance, Charles and Boswell² evaluated the ion-energy distribution function by fitting three Gaussian functions to reconstruct the measured I-V characteristic of a plasma containing an ion beam.

Obtaining an accurate IEDF or EEDF requires choosing between different techniques to compute the first and second derivatives of the I-V characteristics. An ideal method would provide an undistorted derivative with a signal-to-noise ratio (S/N) equal or lower to the raw data S/N. However, in practice, there is a trade-off between noise level reduction and data distortion. Fortunately, analog and numerical methods can typically combine both noise filtering and differentiation of the raw I-V data. Some of the most commonly used numerical routines are presented in this section, namely, the Savitzky–Golay filter, the B-spline fitting, and the Gaussian and the Blackman filters, together with an analog differentiation technique. The Gaussian deconvolution method to obtain IEDFs is also described.

The use of analog differentiator circuits predates numerical differentiation methods, and they are still regularly employed for obtaining EEDFs and occasionally IEDFs.^{4,18,38–40,42} The active analog differentiator used in this study, designed to work with a sweeping frequency of 10 Hz, follows the design described in Takahashi et al.⁴¹ and based on the work by Schoenberg.³⁸ Two cascading operational amplifiers (Renesas CA3140) are each tuned to resonate at $∼ 800$ Hz by adjusting the resistors and capacitors values. This results in a linear gain frequency increase in 20 dB per decade up to the resonance frequency, i.e., the input signals of frequencies up to $∼ 800$ Hz are time-differentiated. Beyond the resonance frequency, the circuit acts as an integrator as the gain linearly decreases until becoming less than unity at $∼ 80$ kHz. Therefore, the analog differentiator configuration also works as a bandpass filter, attenuating high frequency noises. Since real-world circuits deviate from the ideal-case, the linearity was checked to be valid up to at least $∼ 150$ Hz.

The output signals of the two cascading stages of the analog differentiator are proportional to the first and second time derivatives of the probe I(t) signal. Knowing the biasing/discriminator voltage sweep period and voltage range, the output of the differentiator first stage can be converted from being with respect to time to being with respect to voltage as $d I / d V = d t / d V × d I / d t .$ Repeating this operation for the output of the second stage gives the second derivative with respect to voltage. Finally, the absolute value of the probe current derivatives can be retrieved after correcting for the non-unity gain of the circuit at 10 Hz. Alternatively, the analog differentiator outputs can also be converted into absolute value derivatives by numerically integrating the outputs once or twice and scaling the result to the probe I-V characteristic, granted that the characteristic was recorded. Differentiator circuits have the advantage of producing repeatable outputs since no tuning of parameters for different experimental conditions is needed. This simplifies and speeds up the data processing.

A polynomial least squares regression to model the data using a single polynomial of degree N via a least squares method has been previously attempted giving poor results.²⁹ Since the I-V characteristic of interest can be visualized as part exponential and part linear, a single polynomial regression would require a high order to fit the data to a satisfactory level, e.g., N > 20, with N being the polynomial degree.²⁹ However, this leads to Runge's phenomenon and over-fitting,⁴³ which is visible in the results presented by Magnus and Gudmundsson.²⁹ Another shortcoming is the non-locality of the polynomial fit, i.e., distant data points directly impact the local fit.⁴⁴ The SG filter and the B-spline methods both circumvent these problems by fitting low degree polynomials to subsets of the data.

The Savitzky–Golay filter is frequently used for smoothing and differentiating current–voltage characteristics of Langmuir probes and RFEAs.^24,26,27,45 The filtering is achieved by doing a running fitting of a polynomial of degree N using the least squares method on a subset of the (x, y) experimental data of width M > N centered around a given point x_i (forcing M to strictly be an odd number). Evaluating the resulting polynomial at x_i gives the smoothed value Y_i, and the process is repeated for $x i + 1$⁠, etc. Owing to the nature of this process, the first and last $M − 1 2$ data points should be disregarded, or the original data-set extrapolated by $M − 1 2$ extra points at both ends.

Savitzky and Golay²² showed that for equally spaced data points, an analytical solution to the least squares polynomial smoothing exists. Each one of the N + 1 coefficients of the least squares polynomial fitting is computed as a linear combination of the y data points inside the subset of width M. It was further demonstrated that the coefficients of the linear combinations are function of (N, M) only.²² These can be tabulated for all (N, M) pairs, and applying the filter only requires convolution of the coefficients with the raw data. The tabulated coefficients can be found in the original paper by Savitzky and Golay,²² yet most numerical implementations of the method (e.g., in Python or Matlab) readily compute the coefficients when the filtering is performed. A benefit of this property is that the smoothed signal can be differentiated N times by convolution with the Nth derivative of the fitting polynomial, i.e., with the appropriate convolution coefficients.⁴⁶ This avoids having to recourse to numerical finite difference methods which tend to decrease the S/N. This can be understood by considering that the differential operator in the frequency domain is proportional to the frequency itself. For a given M, the pairs $N = 0 – 1 , 2 – 3 , 4 – 5$⁠, etc., give the same convolution coefficients for smoothing the raw data and evaluating the even derivatives, while $N = 1 – 2 , 3 – 4 , 5 – 6$ give the same coefficients for odd derivatives.⁴⁶ Therefore, for a given M, using N = 2 would give an identically smoothed second derivative as using N = 3.

The SG filter acts as a low-pass filter with a flat passband whose cutoff frequency is a function of N and M. The cutoff frequency typically increases with N and decreases with M.⁴⁶ For large M, the benefit of lower cutoff frequencies is compromised by distortions of the higher frequency content of the data. On the other hand, a large N can preserve the narrow features of a signal (such as peaks) but at the cost of increasing the cutoff frequency.⁴⁷ Obtaining an optimal filtering, therefore, requires an iterative process and a trade-off between the values of N and M.

The method employed herein is a numerical alternative proposed by Jauberteau and Jauberteau³⁵ and improved by Bang and Chung.³⁷ Once an I-V characteristic is measured, a Lagrange interpolation is performed on intervals centered on each bias value V, and the polynomial is then evaluated with $V + u ( t )$ to simulate the AC superimposed I-V trace. A fast Fourier transform isolates the harmonics $I ν$ and $I 2 ν$⁠. This is an iterative process running on each individual data point at a time. Granted that u₀ is small (⁠ $≤ 1$ V), the derivatives are well approximated by only considering the first right-hand terms of Eq. (14).

The noise filtering property of the method is a function of u₀. The larger the amplitude u₀, the stronger the noise filtering, but at the cost of potentially distorting the real signal.³⁵ To palliate to this limitation that particularly affects the computation of EEPFs, Bang and Chung³⁷ proposed using an AC signal whose amplitude u₀ linearly increases with the electron energy. They showed that the distortions to the elastic part of the EEPF were minimized while achieving a good filtering on the inelastic region where the signal-to-noise ratio is smaller.

The Gaussian deconvolution method, or Gaussian fitting, has been used to obtain the ion-energy distribution function from the I-V characteristic measured by a retarding field energy analyzer.^2,17,20 This involves an iterative process of integrating a Gaussian curve that models the ion-energy distribution function until the raw I-V data are reconstructed. The fitting can either be manual, or an automatic Gaussian fitting algorithm can be implemented that minimizes the sum of the squared errors between the measured and the reconstructed data.

An advantage of this technique is that it does not involve any filtering or differentiation steps since the IEDF is obtained in a retrograde approach.

The data used in this study have been acquired using Moa,⁵² a radio frequency plasma reactor that is an expansion of Huia.^53,54 Figure 1 shows a schematic of the experimental setup, including the magnetic field topology on axis. It consists of a 150 cm long, 9 cm inner diameter borosilicate glass tube plasma source connected to a 70 cm long, 50 cm diameter steel expansion chamber hosting the pumping system and the pressure gauges. The location of the interface between the plasma source and the expansion chamber is defined as $( r , z ) = ( 0 , 0 )$ cm. A base pressure of $∼ 10 − 7$ Torr is routinely achieved, while the argon working pressure ranges between $0.5 × 10 − 3$ and $1 × 10 − 3$ Torr. The plasma is triggered and energized in the plasma source using a 1–1/3 turns loop antenna connected through an L-type matching network to a 1 kW RF generator working at 27.12 MHz. A static magnetic field is applied either with a pair of movable Helmholtz solenoids placed concentrically with the glass tube which position along the plasma source length can be adjusted, or with a set of three fixed inclined solenoids that can produce a symmetric or deflected magnetic field.⁵² 

Two different experimental configurations are setup to collect the energy distributions. For the collection of the ion-energy distribution functions, the antenna is placed at $z = 18$ cm, and the fixed coils are used to generate an axisymmetric field that divergences in the expansion chamber. This configuration allows the generation of an accelerated ion beam during plasma expansion. The electron-energy distribution functions are measured in the plasma source with the loop antenna located at $z = 50$ cm, and the magnetic field is applied with the movable solenoids placed at $z = 80$ cm. It has to be noted that the measurements are not done simultaneously.

The retarding field energy analyzer used to measure the ion-energy distribution function consists of four mesh grids (namely, earth, repeller, discriminator, and secondary grid) and a nickel collector plate that collects the incoming ion flux. The grids are made of a nickel mesh with a transmission factor of 59%, which are attached to a copper support. The probe orifice has a diameter of 2 mm, and a 0.1 mm thick polyimide sheet is used to electrically insulate the grids and the collector plate. The design of this probe has been extensively used in similar experiments,^2,16,19,55 and a detailed description of the RFEA and its electric circuit used herein are presented in Ref. 56. In order to map the energy of the ions, the voltage of the discriminator grid is swept from 0 to 80 V to obtain the I-V characteristic. It has to be noted that from RFEA measurements, only the energy distribution function of the ions that fall through the probe sheath within an acceptance angle of $± 40 °$ is obtained.² When taking measurements of the ion current in an RF plasma with an RFEA, the collected data can be affected by different mechanisms. In particular, in an RF excited plasma, broadening of the energy distribution caused by RF modulation can be observed^2,16,57 that could lead to peak separation, giving false data on the presence of different ion populations. Additionally, possible ion-neutral collisions (elastic and charge-exchange) occurring in the probe sheath and inside the analyzer could also affect the collected ion energy. However, for the low-pressure case analyzed, collisions inside the RFEA and in the sheath can be considered negligible, i.e., $λ i ≫ s , d$ (where $λ i$ is the ion-neutral mean free path, s is the sheath width in front of the RFEA, and d is the repeller-discriminator grid distance).

As from Langmuir's work, sizing a single Langmuir probe requires the plasma volume disturbed by the probe to be much smaller than the electron mean free path $λ e$ (Ref. 12). This ensures that the probe only disturbs the plasma in ways that can be accounted for by theory. Following recommendations in Ref. 13, Langmuir's theory is satisfied if the probe tip radius $r p$⁠, the tip length $l p$⁠, the tip holder radius $r h$⁠, and the local Debye length $λ D$ all satisfy the condition $r p ln ( π l p / 4 r p ) , r h , λ D ≪ λ e$⁠. Moreover, to guarantee Druyvesteyn's isotropic plasma condition in the presence of an applied magnetic field, the probe tip should be kept perpendicular to the streamlines, and the tip radius is much smaller than the electron Larmor radius.¹² To meet these conditions, the LP used in this work is made from a tungsten tip with $r p = 25 μ$ m and $l p = 3$ mm, mounted at the extremity of a glass pipette tapering down to $r h = 0.75$ mm.

Measuring an EEPF requires biasing the LP above the local plasma potential and can typically draw electron currents on the order of a few milliamperes. Therefore, in order to close the LP electric path and to ensure that all the biasing voltage is developed across the probe sheath, the existence of a low impedance sheath at conductive walls should be verified (see Ref. 12 for more details). Serving this purpose, a grounded aluminum sleeve of adequate surface area was inserted inside the dielectric plasma source tube at the end opposite to the diffusion chamber [see Fig. 1(a)].⁵³ 

With the loop antenna capacitive coupling to the plasma, oscillations of the plasma potential at the applied radio frequency and its harmonics are known to cause severe distortions to LP I-V characteristics.^12,58 While this effect is not easily visible in the probe I-V characteristics, it becomes evident with the second derivative, and, without proper mitigation, it severely compromises EEDF studies. Because the root mean square of $V p$ at the first and second harmonics is larger than $T e / 2$ in the present apparatus, RF compensation was added to the Langmuir probe following design steps described in detail in Sudit and Chen⁵⁸ and Godyak and Demidov.¹² Four resonant chokes (TE Connectivity SC30100KT) and a reference electrode were incorporated in the probe head. The chokes self-resonant frequencies (SFR) were tuned to 27.12 and $54.24$ MHz with shunt capacitors of appropriate values. The filter is housed inside the glass pipette, placed as close as possible to the probe tip to minimize stray capacitance. Shown in Fig. 2, the attenuation spectrum of the completed probe is checked with a spectrum analyzer for eventual detuning induced by the housing and adjacent wiring. The achieved attenuation of at least $− 80$ dB at the first and second harmonics provides sufficient filter impedance to suppress the $V p$ RF oscillations. The reference electrode, connected to the probe tip through a capacitor, reduces the probe sheath impedance to ensure that the measuring tip follows the plasma potential RF oscillations. For chokes containing ferrite cores, care should be taken to ensure that the ferrite does not saturate under the applied magnetic field or else the SFR will be shifted. The chokes employed were confirmed to be unaffected at the applied magnetic field strength of 300 G used in this study. The maximum probe biasing voltage was adjusted to $V bias ≃ V p + 5$ V to avoid excessive electron current collection, which could damage the probe tip or modify its work function. Tip contamination effects were monitored by checking for hysteresis in the up and down sweeps of $V bias$⁠. Probe cleaning by ion bombardment was used when necessary.

The Langmuir probe is mounted on a movable doglegged eccentric dielectric shaft that can span the entire length of the plasma source tube.⁵³ The RFEA is mounted on a movable (rotation and translation) 6.35 mm diameter steel shaft introduced in the top port of the expansion chamber, with the probe orifice facing the plasma source exit to detect any possible directional ion beam. The LP and RFEA currents are deduced by measuring the voltage drop across a resistor with an operational amplifier. The output of the amplifier is then split and fed into both a data acquisition system (DAQ) for digitization and an analog differentiator that outputs the first and second time derivatives. The LP $V bias$ and the RFEA $V D$ are defined by several periods of a triangle wave function with a sweeping frequency of 10 Hz. The data acquisition is triggered by the output on the function generator to ensure signal synchronicity, and the voltages are digitized at 4000 samples per period. One LP measurement is the average of 100 $V bias$ sweeps, while one RFEA measurement averages 200 $V D$ sweeps.

The optimal evaluation of the ion and/or energy distribution functions through numerical methods requires a compromise in the choice of parameters to achieve sufficient noise reduction without causing significant distortion of the fitted I-V characteristic that could lead to large errors in the inferred plasma parameters and loss of information. In this section, the application of the different numerical methods on experimental data is presented as a reference example.

An analysis similar to that used by Fernández Palop et al.²⁸ and Magnus and Gudmundsson²⁹ is used to assess the numerical methods in this work. It is noted that in Refs. 28 and 29, the comparison of the different techniques was conducted on simulated I-V characteristics, while in this study, only experimental data were used. The following parameters, evaluated between the raw probe collected current $I p$ and the processed data $I f$⁠, were considered: the mean-squared error (MSE) and the Shapiro–Wilk (SW) test statistic W. This proposed optimization analysis cannot be employed for the AC superimposition method since its filtering and differentiating steps cannot be separated.

Specifically, the optimization of the IEDFs involves the careful computation of the first and second peaks of the energy distributions, while achieving enough noise reduction when evaluating $d I C / d V D$⁠. An excessive smoothing of the I-V characteristic would modify the voltage at which the I-V curves present the gradient changes that corresponds to the center of two ion populations present. In turn, this would yield an incorrect evaluation of the plasma and ion beam potentials, and an inaccurate derivation of the ion-energy distribution. The dependency of the choice of the numerical processing parameters on the quality of the ion-energy distribution function was studied for each numerical method. Figure 3 shows an example of the effect of the parameter choices on the first derivative for different values of polynomial degree N and number of knots k using the B-spline technique. The mean-squared error and the Shapiro–Wilk test are evaluated for each case, and the values are reported in Table I. When analyzing the dependence of the fitted polynomial degree N while k = 15 is kept constant, the MSE is seen to increase with N, W to decrease, and the floor ripple level to increase (the Runge phenomenon). Since the $I f$ obtained from a B-spline of degree N is N − 1 continuously differentiable, a quadratic B-spline does not provide a smooth first derivative [see Fig. 3(a)]. By extension, if evaluating the EEPF, N = 4 is the lowest B-spline degree that can be used to obtain a smooth second derivative. If the number of knots is too small, the evaluated IEDF is over-smoothed, and the features of the I-V characteristic are distorted (W < 0.996 for N = 3 and k = 5). In contrast, as shown in Fig. 3(f), too high a number of knots causes an insufficient smoothing of the noise, which could make it difficult to determine the location of the peaks. For the data analyzed in this study, a degree of N = 3 and a number of knots k = 15 are chosen for the B-spline polynomial since they result in the optimal IEDF, i.e., lowest curve distortion and sufficient noise reduction.

Similarly, optimal evaluation of the EEDFs requires simultaneously maximizing the energy resolution and the dynamic range (DR). The energy resolution is characterized by the low-energy depletion $δ ε e$⁠, which is defined as the $ε e$ value of the maximum of the EEDF. Because the bulk of the electrons is in the low-energy region of the EEDF, the larger the $δ ε e$⁠, the larger the errors in determining $n e$ and $T eff$ from Eqs. (6) and (7). An EEDF is considered to have sufficient energy resolution if $δ ε e ≤ ( 0.3 − 0.5 ) T e$ (Refs. 12 and 13). The dynamic range of the EEDF is defined as the ratio of its maximum to its minimum resolvable value. A dynamic range $≥ 60$ dB (equivalently of three orders of magnitude and above) is necessary to resolve high-energy electrons in the inelastic range of the EEDF, which is detrimental in investigating excitation, ionization, and transport processes.^12,13 An EEDF is often separated into its elastic and inelastic ranges, defined as $ε e < ε e *$ and $ε e > ε e *$⁠, respectively. $ε e *$ is the lowest excitation energy of the working gas, e.g., 11.55 eV for argon.⁶² Minimizing $δ ε e$ and maximizing DR are, therefore, paramount in the process of determining the optimal numerical parameters, together with minimizing the MSE and maximizing the value of W.

Figure 4 shows an example of the impact of varying the polynomial degree N and moving window size M of the Savitzky–Golay filter when applied to a compensated LP I-V characteristic acquired with the experimental apparatus. For each (N, M) pair, the corresponding MSE and Shapiro–Wilk test statistic W are given in Table II. In Figs. 4(a)–4(d)^, N = 2, the effect of increasing M is shown. As expected for an SG filter, increasing M is beneficial from the dynamic range point-of-view since it improves from 13 to 46 dB owing to the decreasing cutoff frequency. However, this is accompanied by distortions of the signal, visible from the worsening of the low-energy resolution $δ ε e$⁠, the MSE and W. Fixing M and increasing N to the next even integers (since N = 2 and N = 3 produce the same output), the signal distortion improves again at the cost of decreasing the DR, as shown in Figs. 4(d)–4(f) and Table II. Regarding the Shapiro–Wilk test, all the cases with M = 251 failed the normality test with W < 0.996. Overall, trading-off between $δ ε e$⁠, the DR, the MSE, and W, the best performing filter for these data is for $( N = 2 , M = 115 )$⁠. It is noted that none of the (N, M) pairs produced an EEPF fulfilling the aforementioned requirements, since the dynamic range never reaches three orders of magnitude. Different (N, M) pairs can be used whether the goal of obtaining the EEPF lies in the study of elastic or inelastic processes. Alternatively, the outputs of two SG filter applications could be combined into a high resolution and high dynamic range numerically derived EEPF, with the drawback of a more complicated post-processing.³⁰ 

Since this optimization analysis cannot be applied to the AC superimposition technique, the optimal application of this method was found by maximizing the energy distribution dynamic range when varying the AC perturbation signal amplitude u₀.

As mentioned in Sec. I, the ion-energy distribution function can be derived from the first derivative of the measurements acquired with a retarding field energy analyzer. The I-V characteristic presented in this study was taken in the expansion chamber with the RFEA orifice located 19 cm downstream of the magnetic nozzle (⁠ $z P = − 19$ cm). The experimental conditions were maintained fixed at an RF power of 250 W, a reflected power of $∼ 1 %$ or lower, an argon pressure in the expansion chamber of $0.5$ mTorr, and a maximum magnetic field on axis of $B z , max = 340$ G. Argon is injected from the inlet on the positive z end of the glass tube to reproduce the configuration from Charles and Boswell.² Under these conditions, the plasma density peaks axially under the solenoids at z = 8 cm reaching a value of $≈ 10 12 cm − 3$⁠, which then decreases to $∼ 10 10 cm − 3$ in the expansion chamber. The data processing of the collected current $I C ( V D )$ measured with the energy analyzer and the subsequent computation of the ion-energy distribution function followed the approach described in Sec. V A.

Table IV summarizes the results of the parametric analysis that provide the smoothest IEDF possible for the numerical methods analyzed, while keeping intact the peak features. Figure 5 shows the ion-energy distribution functions obtained with the different analog and numerical techniques from experimental data. The result of the optimal AC superimposition method implementation is shown in Fig. 6(a) for which $u 01 = 2.75$ and $u 02 = 2.85$ V were used. The IEDFs plotted are normalized by their peak value, and the position of the two peaks, i.e., the local plasma potential and the ion beam potential, is also reported. The raw I-V curve, shown in Fig. 5(a), was obtained by averaging 200 consecutive probe sweeps and has an S/N of approximately 66 dB. It is plotted with the fitted Gaussian deconvolution to illustrate one of the data processing methods. As seen in Figs. 5 and 6(a), the ion-energy distribution functions show a double-peak feature indicating the presence of a bi-population of ions: the first peak represents the background ion population (at zero energy) in which location can be considered the local plasma potential $V p$⁠, while the second peak corresponds to the accelerated ion beam population with potential $V B$ and energy $ε B = e ( V B − V p )$⁠. Comparing measurements taken with a source-facing and a radial-facing RFEA shows that the IEDF does not present a peak separation effect due to RF modulation, but it effectively detects the presence of a directional ion beam.⁵² It is also noted that the ion-neutral mean free path for the operating pressure of $0.5$ mTorr is shorter than the distance between the source exit and the probe position, i.e., $λ i ∼ 10$ and $z p = − 19$ cm. Thus, the amplitude of the accelerated ion beam component is reduced at the measurement location, while the local, background ion population is enhanced due to charge-exchange collisions. Furthermore, elastic collisions could affect the beam direction, further reducing the collected beam current.

Table III presents the key results obtained with the different data processing techniques, namely, the local plasma potential $V p$⁠, the ion beam potential $V B$⁠, the mean-squared error (MSE), and the Shapiro–Wilk test statistic W. The Gaussian deconvolution method required the fitting of three Gaussian curves to accurately model the two ion populations. Comparing the analog differentiation method and the numerical techniques, it is evident that both the plasma and ion beam potentials are approximately unchanged with a standard deviation for $V p$ and $V B$ of $± 0.22$ and $± 0.21$ V, respectively. For the collected data, all the concerned smoothing methods perform well in terms of mean-squared error and the SW test statistic. The MSEs are $≤ 9.25 × 10 − 4 μ A 2$ for all techniques, while the calculated W values are all higher than 0.998. The Saviztky–Golay filter yields the lowest MSE, while the Gaussian deconvolution method performs the best in the SW test, resulting in the highest W. It has to be noted that the Blackman filter showed an overshoot of the fitted polynomial at the edges caused by the zero padding in the convolution; this was addressed by disregarding the initial data points as the plasma characteristics in that region were not of interest. Finally, it is clear that while the AC superimposition method delivered plasma parameters within the standard deviation margins of the other methods, it has also produced a noisier IEDF. This performance is to be relativized by comparing it to the outcome of the central difference method used on the raw I-V data, shown as the background low contrast curve in Fig. 6(a). As expected, the latter method has been unable to resolve the IEDF owing to the amplification of the measurement noise.

The RF compensated LP data were acquired with the apparatus configured with the loop antenna placed at $z = 50$ cm and the movable Helmholtz solenoids at $z = 80$ cm. The RF power was set to 200 W, and the reflected power kept at $∼ 1 %$⁠. A maximum magnetic field on axis of $B z , m a x = 300$ G was used. An argon working pressure of $1$ mTorr was set by injecting the gas from the expansion chamber in order to minimize neutral pressure gradient effects along the glass tube.⁵³ These conditions create an inductively coupled plasma exhibiting a symmetrical single peaked axial plasma density gradient ranging from $10 12 cm − 3$ under the solenoids to $10 10 cm − 3$ at the extremities of the glass tube.^53,54,63 The LP was placed on axis halfway between the antenna and the solenoids, i.e., at $z = 65$ cm.

Figures 6(b) and 7 show the EEPFs obtained with the various methods. For the analog differentiation technique, the differentiator circuit output was scaled into an absolute value EEPF using the double integral method described in Sec. III A and Eq. (3) to produce the EEPF in Fig. 7(b). The processing of $I p ( V bias )$ to obtain the EEPF from the numerical methods was conducted following the steps described in Sec. V A. Table IV gives the resulting optimal numerical parameters. For the AC superimposition method, the optimal EEPF in Fig. 6(b) was obtained with u₀₁ varying from 1.5 to 8 V and u₀₂ from 1.25 to 7.75 V, between $ε e = 0$ and $ε e = 50$ eV. For reference, Fig. 6(b) also illustrates the inability of the central difference method to resolve any EEPF feature.

Figure 7(a) shows the unprocessed I-V LP characteristic together with the optimized Savitzky–Golay filter output from Fig. 4, i.e., for $( N = 2 , M = 115 )$⁠. As reference, the signal-to-noise ratio of the shown LP characteristic is 48 dB. In Figs. 6(b) and 7(b)–7(f), the EEPFs obtained from the various methods and the deduced plasma parameters are shown. $T eff$ and $n e$ are calculated by integrating the EEPF according to Eqs. (7) and (6), respectively, and $V p$ is obtained as the zero crossing of the second derivative. Their values are reported together with the dynamic range, energy resolution for all methods as well as the MSE and W for the relevant methods in Table V. The maximum variance in the obtained values of $V p$ is less than 2%. Likewise, for $T eff$ and $n e$⁠, their maximum variances are $∼ 13 %$ and $∼ 8 %$⁠, respectively. Thus, in terms of estimating the plasma parameters, all six methods have comparable performances, i.e., within the absolute error tolerances usually associated with LP measurements.^4,13

In terms of energy resolution and dynamic range, different methods perform best depending on whether the focus is on resolving low-energy or high-energy electrons dynamics and properties. The B-spline fitting in Fig. 7(d) resulted in the best energy resolution at $δ ε e ≃ 2.3$ eV and is the only method to satisfy the $δ ε e ≤ ( 0.3 − 0.5 ) T e$ condition. The second-best performing method is the analog differentiation with $δ ε e ≃ 3.2$ eV, i.e., 0.375 eV away from satisfying the energy resolution condition. Regardless of the employed method, the energy resolution is ultimately limited by the absence of internal resistance compensation in the present Langmuir probe implementation, in particular of the plasma-wall sheath impedance, or by the probe tip surface condition.¹³ Techniques to further improve the resolution are given in Refs. 13 and 30.

The analog differentiation method in Fig. 7(b) resulted in an EEPF with a significantly larger dynamic range of approximately three orders of magnitude compared to two orders of magnitude for the numerical smoothing methods and approximately one order of magnitude for the AC superimposed technique. The analog differentiation is, therefore, the only method presented, which can reliably resolve the inelastic domain and can provide confidence on the nature of the electron distribution function, i.e., closely matching a single-Maxwellian distribution for $ε e$ up to 30 eV. This is further highlighted in Fig. 7(b) by the linear regression of the EEPF, shown by a dashed red line. The same linear regression applied to the numerically derived EEPFs highlights that the reduced DR and nonphysical oscillations in the case of the B-spline would make it more difficult to reach a conclusion on the nature of the electron distribution function. The lower performance of the AC superimposition technique in terms of dynamic range compared to the example given in Bang and Chung³⁷ is to be attributed to the difference in respective noise levels. Indeed, the presently used LP trace has a signal-to-noise ratio of 48 dB compared to the 80 dB of the trace in Ref. 37.

In term of minimizing the distortion to the original I-V curve, the Savitzky–Golay and the B-spline are the best performing numerical methods for this dataset, as shown by their lower mean square errors and Shapiro–Wilk test statistics in Table V. The Gaussian and Blackman window filters both resulted in noticeable distortions of the data (W < 0.996, i.e., their residuals are not normally distributed), in particular around the plasma floating potential and plasma potential portions of the I-V curve.

Effects of the experimental conditions on the relative performances of the methods reviewed here have been evaluated by repeating the same methodology with the LP at a different location in Huia. With the probe located off-axis at $z = 70$ and $r = 1.8$ cm, the S/N was found to have degraded by $5$ dB. It was found that the same overall conclusions can be drawn with respect to the relative performances of the different methods. A shared effect of decreased dynamic range, which required a stronger filtering of the data with the numerical methods, was observed. Details of this analysis can be found in  Appendix.

Finally, it is worth noting that following the guidelines in Ref. 13, the low-energy gap $δ ε e$ impacts can be corrected by extrapolating the EEPF from its linear portion down to $ε e = 0$ eV when the right conditions are satisfied. Because electron–electron collisions are the dominant processes for thermalizing electrons and making the distribution Maxwellian, and since the frequency of these collisions is inversely proportional to $ε e 3 / 2$⁠, if the EEPF is Maxwellian for some $ε e$ greater than the energy of EEPF peak, then the EEPF is likely to be Maxwellian at lower energies.^13,64,65 This extrapolation can be achieved with the linear regressions shown in Fig. 7 as dashed red curves, starting the extrapolation from the lowest energy at which the dashed red curve contacts the EEPF. For example, when applying this process to the analog EEPF in Fig. 7(b), the effective electron temperature obtained with the corrected EEPF is $T eff ′ = 4.32$ eV and the corrected electron density $n e ′ = 2.57 × 10 11 cm − 3$⁠. These values are slightly different to the values in Table V. Indeed, the electron temperature obtained from the slope of the linear regression is $T e , fit = 4.2$ eV, providing evidence that the low-energy extrapolation has closely recovered the original EEPF shape. More details on the low-energy gap correction and its outcomes are discussed in  Appendix. It is relevant to note that the electron density measured here at $z = 65$ cm via the Druyvesteyn method is about half of the one found via the ion saturation method in Ref. 54 at the same location and under the same operating conditions. This illustrates another point highlighted in Godyak¹³ and Zhang et al.,⁶⁶ i.e., the Langmuir method and the ion saturation method often overestimate the plasma density. This is usually the result of electron distributions which depart from the presumed Maxwellian, and of practical difficulties to account for probe-sheath effects affecting the effective ion collection area.

This review outlines some of the most frequently employed data processing methods to obtain ion and electron-energy distribution functions from I-V characteristics measured with a retarding field energy analyzer and an RF compensated Langmuir probe, respectively. Design recommendations relevant to measurements in RF magnetized plasmas are presented for both instruments. After describing the different processing techniques, an example of a parametric analysis to obtain optimal IEDFs and EEPFs is presented and applied to experimental datasets collected in a magnetized RF plasma apparatus. The recommended parameters optimization involves checking the mean-squared error and the normal distribution of the residuals (i.e., the Shapiro–Wilk test), coupled with a visual assessment process to ensure sufficient noise reduction with minimal distortion of the curves.

With respect to the IEDF, the analysis focused on the evaluation of first derivatives that accurately represented a plasma containing two populations of ions. For the dataset under study, both the analog and numerical methods provide IEDFs with consistent values of the local plasma and the ion beam potentials. Regarding the numerical smoothing methods analyzed, they all perform well in terms of MSE and W; in particular, the Saviztky–Golay filter and the Gaussian deconvolution method delivered the best results in terms of obtaining accurate IEDFs with the least data distortion. While the AC superimposition method delivered an IEDF with a noise level higher than other methods, its features are till clearly defined.

The EEPFs obtained from each method are assessed in terms of maximizing the energy resolution and dynamic range while minimizing distortions of the Langmuir probe I-V curve. Among the numerical methods analyzed in the example application, the Savitzky–Golay filter and the B-spline fitting performed best in terms of distortion and energy resolution. The numerical methods delivered a low dynamic range of approximately two orders of magnitude and one order of magnitude for the AC superimposition method. This is to be related to the experimental dataset employed, having a low S/N of 48 dB. This could make it challenging to confidently identify the nature of the electron distribution function and to study inelastic electron processes. This is in contrast with the EEPF obtained from the same I-V curve using the analog differentiator which showed a dynamic range of approximately three orders of magnitude and can be clearly identified as being Maxwellian up to $∼ 30$ eV. When the intent of calculating the EEDF is to obtain the plasma parameters, all methods explored are found to deliver equivalent results for the particular dataset analyzed. It has to be stressed that the performances of the data processing techniques presented are dependent on the signal-to-noise ratio of the sample under analysis; thus, one has to assess the optimal method for their particular experimental apparatus, diagnostics, and dataset. The analog differentiator technique provides an appreciable reduction in the IEDF and EEDF processing time compared to the numerical methods since no dataset specific optimization is required.

See the supplementary material for ASCII tables containing the raw I-V characteristics measured with the RFEA and the Langmuir probe used in this study.

The authors would like to thank Alexander Bennet and Dimitrios Tsifakis for useful discussions on probe design and operation. This research was partly supported by the New Zealand Space Agency under Grant No. MBIE#00008060 as well as by the Asian Office of Aerospace Research and Development (AOARD), the international office of the Air Force Office of Scientific Research (AFOSR), under Grant No. FA2386-19-1-4012.

The authors have no conflicts to disclose.

Antonella Caldarelli and Félicien Filleul contributed equally to this work.

Antonella Caldarelli: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Resources (lead); Software (lead); Validation (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). Félicien Filleul: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Resources (lead); Software (lead); Validation (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). R. W. Boswell: Conceptualization (lead); Investigation (supporting); Methodology (lead); Resources (equal); Supervision (equal); Validation (lead); Visualization (lead); Writing – review & editing (supporting). Christine Charles: Conceptualization (lead); Funding acquisition (equal); Investigation (supporting); Methodology (lead); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – review & editing (supporting). Nicholas Rattenbury: Funding acquisition (equal); Project administration (equal); Resources (supporting); Supervision (equal). John Edward Cater: Funding acquisition (equal); Project administration (equal); Resources (supporting); Supervision (equal).

The data that support the findings of this study are available within the article and its supplementary material.

The respective performances of the methods compared in Sec. V C are expected to be influenced by experimental conditions specific to the locus of the measurement. Repeating the methodology at a different location in the present plasma discharge might lead to a change in the respective performances of the methods. To evaluate this, the RF compensated probe was placed off-axis at $z = 70$ and $r = 1.8$ cm with Huia operating under the same conditions as in Sec. V C. At this location, the S/N was measured to be about 44 dB, i.e., 5 dB lower than on axis at $z = 65$ cm. The numerical parameters used for this analysis are given in Table VII. An analysis of Fig. 8 and Table VI shows that the respective performances of the methods are mostly unaffected by the change in signal-to-noise ratio: the analog method still results in the highest dynamic range and second lowest low-energy gap after the B-spline method. The Gaussian and Blackman window filters are still the lowest performing methods in terms of signal distortion and dynamic range.

A noticeable difference with the previous dataset is the overall reduction in dynamic range, which is to be expected with the increased noise level. Also consistent with this expectation is the fact that the optimum numerical parameters (see Table VII) are ones which resulted in a more aggressive filtering/smoothing of the I-V characteristic. Compared to the performances of the methods at the different location, the low-energy gap $δ ε e$ is consistently larger. A plausible reason is the more aggressive filtering, which has resulted in a rounding-off of the second derivative peak. Table VIII shows the effective electron temperatures and densities obtained from each method after performing a linear extrapolation of the EEPFs to $0$ eV as a way to correct for the low-energy gaps, as described in Sec. V C. It can be seen that the variability in $T eff$ and $n e$ between the different methods, visible in Table VI, is greatly reduced once the low-energy gap is corrected. This illustrates the fact that the performances of a method in terms of determining the plasma parameters are a strong function of its ability to accurately resolve the low-energy section of an EEPF. This is to be expected since the great majority of the electrons in a distribution are those at low energy.

Finally, it is interesting to mention that the analog EEPF in Fig. 8(b) appears to be bi-Maxwellian with an inelastic electron population warmer than the bulk of the distribution. This feature being only resolved by the analog differentiator method is consistent with the fact that this method delivered the highest resolved electron energy at $ε e ≃ 40$ eV, owing to the method's higher dynamic range, compared to 20–25 eV for the numerical methods. The presence of a hotter high-energy population at this location is interpreted as the result of the transport of electrons from the antenna skin layer along the magnetic field line crossing this region, whose radial position is close to $1.8$ at $z = 70$ cm.

As a conclusion, different experimental conditions do influence the processing methods performances. Likewise, one can expect that employing differently constructed probes might also change the analysis outcomes. As such, it is recommended to the experimentalists to investigate which method performs best for their own experimental and plasma conditions, using the suggested methodology as a guideline.