Retarding field energy analyzers and Langmuir probes are routinely used to obtain ion and electron-energy distribution functions (IEDF and EEDF). These typically require knowledge of the first and second derivatives of the current–voltage characteristics, both of which can be obtained using analog and numerical techniques. A frequent problem with electric-probe plasma diagnostics is the noise from the plasma environment and measurement circuits. This poses challenges inherent to differentiating noisy signals, which often require prior filtering of the raw current–voltage data before evaluating the distribution functions. A review of commonly used filtering and differentiation techniques is presented. It covers analog differentiator circuits, polynomial fitting (Savitzky–Golay filter and B-spline fitting), window filtering (Gaussian and Blackman windows) methods as well as the AC superimposition and Gaussian deconvolution routines. The application of each method on experimental datasets with signal-to-noise ratios ranging from 44 to 66 dB is evaluated with regard to the dynamic range, energy resolution, and signal distortion of the obtained IEDF and EEDF as well as to the deduced plasma parameters.

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 rates7 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.14 

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.

The ion-energy distribution function has usually been computed by obtaining the first derivative of the collected current over the discriminator voltage.2,5,6,15–20 The total ion current measured by a retarding field energy analyzer is given by
I C = t 4 A e v 0 + v i f ( v i ) d v i ,
(1)
where t is the grid transmission factor (four grids in the case of the RFEA used in this work), A is the area of the probe orifice, v i is the component of the ion velocity perpendicular to the probe face, and f ( v i ) is the one-dimensional ion velocity distribution at the entrance of the probe.2, v 0 is the minimum velocity required for the particle with energy ε i to overcome the potential barrier in front of the probe and satisfy conservation of energy. It can be shown that, by assuming a one-dimensional velocity distribution,15 the ion-energy distribution function g i ( ε i ) is proportional to the negative derivative of the collected current I C with respect to the discriminator grid voltage V D,
g i ( ε i ) = 1 t 4 A o e 2 m i 2 e V d I C d V D ,
(2)
where m i is the ion mass. Thus, the IEDF can be calculated simply by differentiating the measured current with respect to the discriminator voltage. For a single Gaussian distribution, the voltage at which the IEDF peak occurs is defined as the local plasma potential, V p. If a population of accelerating ions is present, as in the case of a plasma expanding through a double-layer2,19,21 and for the considered gas pressures (i.e., ∼0.2–1 mTorr), the IEDF would show a second peak at a higher energy located at the ion beam potential, V B.
Druyvesteyn14 showed that for an isotropic plasma, the electron-energy probability function (EEPF) g ep is proportional to the second derivative of the Langmuir probe electron current I e with respect to the probe biasing voltage V bias as
g ep ( ε e ) = 2 2 m e A p e 3 d 2 I e d V bias 2 .
(3)
Here, ε e is the electron energy, A p is the probe area, e is the elementary charge, and m e is the electron mass. The EEDF relates to the EEPF as g ed ( ε e ) = g ep ( ε e ) ε e to give
g ed ( ε e ) = 2 2 m e e V A p e 3 d 2 I e d V bias 2 ,
(4)
where V = V p V bias and V p is the local plasma potential. For an isotropic electron gas, the EEDF and EEPF contain the same information as the electron distribution function, and as such the terms are often used interchangeably in the literature.12,13 Working with the EEPF rather than with the EEDF presents practical advantages. From the expression of a Maxwellian EEPF,7 
g ep ( ε e ) = n e 2 π T e 3 / 2 exp ( ε e T e ) ,
(5)
it can be seen that a Maxwellian g ep ( ε e ) is a straight line with a slope equal to 1 / T e when displayed on a semi-log plot. Displaying d 2 I e / d V bias 2 on a semi-log scale, therefore, allows to quickly appreciate whether a measured distribution is Maxwellian or not, and if so to easily obtain T e.
The choice of normalization of Eq. (4) is such as to retrieve the total electron number density from
n e = 0 g ed ( ε e ) d ε e .
(6)
Then, the effective electron temperature T eff can be calculated from the average electron energy ε e as
T eff = 2 3 ε e = 2 3 n e 1 0 ε e g ed ( ε e ) d ε e .
(7)
If the distribution is Maxwellian, T eff = T e. Finally, the local plasma potential V p is found from the zero crossing of d 2 I e / d V bias 2. From this, it is clear that the Druyvesteyn method is more general than the Langmuir method since once d 2 I e / d V bias 2 is measured, one can obtain the plasma parameters without a priori knowledge of the nature of the electron distribution. Owing to the limitations discussed by Godyak and Demidov12 regarding subtracting the ion current from the I-V characteristic, the full probe characteristic is used to compute the distribution function. It is, indeed, argued that doing otherwise will likely generate errors in the EEPF, and that the correction would only be significant if the dynamic range of the second derivative was to exceed five orders of magnitude.12,13

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.24 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 Fredriksen27 to smooth the I-V curve acquired with an RFEA and obtain the ion-energy distribution function. Fernández Palop et al.28 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 Gudmundsson29 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.30 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 MacGregor31 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 analogically31–34 or numerically.35 The analog AC superimposition method was also used to obtain IEDFs by Hatakeyama et al.36 Bang and Chung37 improved the numerical AC superimposition method of Jauberteau and Jauberteau35 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 Schoenberg38 has been often used to obtain EEPFs using appropriate electronics circuitry.3,39–42 Takahashi et al.40 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.18 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 Boswell2 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.41 and based on the work by Schoenberg.38 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.29 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.29 However, this leads to Runge's phenomenon and over-fitting,43 which is visible in the results presented by Magnus and Gudmundsson.29 Another shortcoming is the non-locality of the polynomial fit, i.e., distant data points directly impact the local fit.44 The SG filter and the B-spline methods both circumvent these problems by fitting low degree polynomials to subsets of the data.

1. Savitzky–Golay filter

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 xi (forcing M to strictly be an odd number). Evaluating the resulting polynomial at xi gives the smoothed value Yi, 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 Golay22 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.22 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,22 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.46 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.46 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.46 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.47 Obtaining an optimal filtering, therefore, requires an iterative process and a trade-off between the values of N and M.

2. B-spline fitting

The piecewise polynomial approach used for data smoothing is the cardinal or uniform B-spline method.48 This polynomial fitting method avoids Runge's phenomenon by making use of a piecewise approach: the data are split in k + 1 pieces joined by k equidistant points called knots. Like other polynomial regressions, the fitting polynomial f(x), or spline function, of degree N can be expressed as a linear combination of coefficients βi and k + N + 1 basis functions (polynomials) B i , N + 1 ( x ),
f ( x ) = i β i B i , N + 1 ( x ) ,
(8)
where the ith B-spline basis function of degree N is defined for a given sequence of knots as
B i , 1 ( x ) = { 1 for x i x < x i + 1 , 0 otherwise
(9)
and
B i , N + 1 ( x ) = x x i x i + N x i B i , N ( x ) + x i + N + 1 x x i + N + 1 x i + 1 B i + 1 , N ( x ) ,
(10)
for N > 1. Spline functions can, therefore, be calculated on different bases, and the B-spline is a particular case for basis functions with local support, i.e., taking non-zero values only inside the interval on which they are defined, which spans N + 2 knots. As a result, this method results in a higher numerical stability compared to other splines.44 The spline function f is fitted to the data by using the least-square principle to calculate the k + N + 1 coefficients βi. Additionally, the continuity conditions need to be verified, such that the function is N − 1 times differentiable at the knots and that f = f = 0 at the endpoints. The larger the number of knots, the higher the accuracy of the data fitting. However, more knots may result in over-fitting and/or poor noise rejection. A small number of knots could also over-smooth the data. As for the SG filter, the first and second derivatives of f can be analytically derived.
Using window functions for noise filtering is common practice in signal processing. Applying a convolution operation between the signal and a window function constitutes one of the basis of so-called finite impulse response filters.49 The convolution of two continuous functions a and b evaluated at the variable t is defined as
( a b ) ( t ) = a ( τ ) b ( t τ ) d τ .
(11)
Equivalent definitions of the convolution and its properties can be given for discreet functions, which are used in practice for numerical data processing.49 An interesting result of complex analysis is the convolution theorem, which states that the convolution of two functions in the time-domain is equivalent to the product of their respective Fourier transforms in the frequency domain.49 Therefore, the convolution of a dataset with a given window function is equivalent to performing a discreet Fourier transform filtering. Owing to the identity ( a b ) = a b = a b , the signal to be filtered can be convolved with the first and second derivatives of the window function in order to differentiate the signal once and twice, respectively. This avoids resorting to the noise-amplifying central difference method. For treating plasma probes signals, the Gaussian filter and the Blackman window are the two most often encountered window functions.

1. Gaussian filter

The Gaussian filter was developed by Hayden50 who introduced a new smoothing routine involving the use of a Gaussian distribution function as a filter. Fernández Palop et al.28 and Magnus and Gudmundsson29 used this algorithm to smooth the I-V characteristic measured with a Langmuir probe and obtain the electron-energy probability function. The working principle of the Gaussian filter is based on the assumption that the response function of the instrument and its electronics can be approximated as a Gaussian distribution with a given standard deviation σ,
g ( x ) = 1 σ 2 π exp ( x 2 2 σ 2 ) .
(12)
Because the noise is not correlated with the instrument response function, g can be used for noise suppression. The filtering and the differentiated signals are evaluated by performing a convolution between the data and g(x) and its derivatives, respectively.

2. Blackman window filter

The Blackman window w is a standard moving average window, having zero value outside a chosen interval. It has been previously used to filter I-V characteristics in EEPF applications by Magnus and Gudmundsson29 and Roh et al.30 The Blackman window is defined by only one filtering parameter, i.e., the window size M,
w ( n ) = 0.42 0.5 cos ( 2 π n M ) + 0.08 cos ( 4 π n M ) , 0 n M .
(13)
It is characterized by the highest reduction in the sidelobe level (down to 58 dB) when compared to the other window functions (e.g., Hanning).51 However, it has a wider mainlobe resulting in a less sharp transition band in the passband of the filter that can be improved by choosing a wider window size. The first and second derivatives of Eq. (13) can be convolved with the measured data to obtain the differentiated filtered signals.
In the original analog implementation of this method, an alternating potential u ( t ) = u 0 sin ( 2 π ν t ) of amplitude u0 and frequency ν is superimposed onto the biasing voltage V of an electric probe.31 Provided that u0 is small with respect to V, the AC signal acts as a perturbation, and its resulting first and second harmonics contributions to the probe trace (noted I ν and I 2 ν) are proportional to the first and second derivatives of the I-V curve through the Taylor expansions neglecting high order terms,35,
I ν = u 0 d I ( V ) d V + u 0 3 8 d 3 I ( V ) d V 3 + , I 2 ν = u 0 2 4 d 2 I ( V ) d V 2 + u 0 4 48 d 4 I ( V ) d V 4 + ,
(14)
where I ν and I 2 ν can be isolated from the probe trace by filtering it in the frequency domain.33,34

The method employed herein is a numerical alternative proposed by Jauberteau and Jauberteau35 and improved by Bang and Chung.37 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 u0 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 u0. The larger the amplitude u0, the stronger the noise filtering, but at the cost of potentially distorting the real signal.35 To palliate to this limitation that particularly affects the computation of EEPFs, Bang and Chung37 proposed using an AC signal whose amplitude u0 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.

To correct for distortions in the inelastic range resulting from using high amplitude perturbations, it is necessary to take into account the higher order derivative terms. This is done by using two AC signals of different amplitudes u01 and u02 and performing the computation with each signal. In that case, the derivatives are37 
d I ( V ) d V = u 01 3 I ν , 2 u 02 3 I ν , 1 u 02 u 01 3 u 01 u 02 3 , d 2 I ( V ) d V 2 = 4 u 01 4 I 2 ν , 2 u 02 4 I 2 ν , 1 u 01 4 u 02 2 u 02 4 u 01 2 .
(15)
In what follows, the AC superimposition method with linearly varying u0 was used for computing the EEPF while a perturbation of constant amplitude was used for the IEDF since IEDFs have unobvious correlations between ion-energy and signal-to-noise ratio.

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.

For a single ion population, the ion-energy distribution function can be modeled with a Gaussian function,
g ( V D ) = a e ( V D V peak c ) 2 ,
(16)
where the parameters to be fit are the amplitude a, the curve width c, and the location of the Gaussian peak V peak. When a second ion population is present, as in the case of an accelerated ion beam in a collisional plasma, the IEDF requires the fitting of two or three Gaussian curves, which are summed to accurately model the measured I-V curve.2,19
The automatic fitting algorithm in the presence of an ion beam involves the initial estimates of the parameters a, c, and V peak in Eq. (16) for each Gaussian curve to be fitted (in this work, three Gaussians were necessary to accurately model the measured data presented in Sec. V B). The integral of the summed Gaussian curves G ( V D ) gives the reconstructed I-V curve. The automatic algorithm then searches for the parameters of the fitted Gaussian curves that minimize the difference between the measured I C V D curve and the integrated data, i.e.,
min V D F ( V D ) ,
(17)
where F ( V D ) is defined as
F ( V D ) = ( I C V D , min V D , max G ( V D ) d V D ) 2 .
(18)

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,52 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.52 

FIG. 1.

Schematic of the experimental setup (a), and plot of the magnetic field on axis generated by the fixed (solid line) and movable coils (dashed line).

FIG. 1.

Schematic of the experimental setup (a), and plot of the magnetic field on axis generated by the fixed (solid line) and movable coils (dashed line).

Close modal

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.2 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 observed2,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.12 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)].53 

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 Chen58 and Godyak and Demidov.12 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.

FIG. 2.

The attenuation spectrum of the RF compensated probe used in this study.

FIG. 2.

The attenuation spectrum of the RF compensated probe used in this study.

Close modal

The Langmuir probe is mounted on a movable doglegged eccentric dielectric shaft that can span the entire length of the plasma source tube.53 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.28 and Magnus and Gudmundsson29 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.

When comparing different regression models, a common mistake is to assume that an optimal processing technique would return the lowest mean-squared error. However, because the MSE value can range between zero and any larger number depending on the scale of the variables, it does not provide a standardized metric to assess how well the processing actually models the observed data. Moreover, the MSE on its own does not provide a way to highlight under-smoothing and over-fitting as these situations could return the lowest possible MSE. Thus, an analysis of the distribution of the residuals ( I p I f) is necessary to complement the MSE. In particular, verifying that the residuals are normally distributed ensures that the data processing does not significantly distort the I-V features.59 The Shapiro–Wilk test checks the null-hypothesis that the residuals are normally distributed by comparing the computed statistic 0 < W < 1 with tabulated values.60,61 The W statistic is defined as
W = ( i = 1 n a i y i ) 2 i = 1 n ( y i y ¯ ) 2
(19)
where yi is an ordered data sample of size n with mean y ¯ and ai is the tabulated coefficients. For a given significance level, if W is smaller than the tabulated value, then the null-hypothesis can be rejected, i.e., the residuals are not normally distributed. Since filtering and fitting always somewhat distort the data, one should aim to maximize W and reevaluate the data processing when the null-hypothesis is rejected. Taking a 0.01 significance level and given the samples size used in this study, the null-hypothesis is rejected if W < 0.996 (Ref. 61). Therefore, a parametric analysis is conducted together with a visual optimization of the IEDF and EEPF for determining the optimal parameters for each numerical method. The assessment requires a trade-off between minimizing the MSE, maximizing W, and minimizing alterations of the raw I-V characteristic features.

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.

FIG. 3.

Normalized first derivative of the I-V characteristic measured by the RFEA obtained with the B-spline polynomial for various values of degree N (a)–(c), and number of knots k (d)–(f).

FIG. 3.

Normalized first derivative of the I-V characteristic measured by the RFEA obtained with the B-spline polynomial for various values of degree N (a)–(c), and number of knots k (d)–(f).

Close modal
TABLE I.

Mean-squared error and Shapiro–Wilk test for different values of degree N and number of knots k of the B-spline applied to the RFEA data in Fig. 3.

N k MSE / 10 4 ( μ A 2 ) W
45.62  0.9935 
15  9.25  0.9989 
35  8.92  0.9985 
80  8.87  0.9985 
15  9.42  0.9987 
15  9.75  0.9986 
N k MSE / 10 4 ( μ A 2 ) W
45.62  0.9935 
15  9.25  0.9989 
35  8.92  0.9985 
80  8.87  0.9985 
15  9.42  0.9987 
15  9.75  0.9986 

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.62 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.30 

FIG. 4.

EEPFs obtained from applying the Savitzky–Golay filter on an LP I-V curve for N = 2 and increasing window size M (a)–(d) and M = 251 with increasing value of N (e) and (f). The horizontal dashed lines mark the minimum resolvable values of the EEPFs.

FIG. 4.

EEPFs obtained from applying the Savitzky–Golay filter on an LP I-V curve for N = 2 and increasing window size M (a)–(d) and M = 251 with increasing value of N (e) and (f). The horizontal dashed lines mark the minimum resolvable values of the EEPFs.

Close modal
TABLE II.

Mean-squared error (MSE) and Shapiro–Wilk test (W) for various values of degree N and window size M of the Savitzky–Golay applied to the LP data shown in Fig. 4.

N M MSE / 10 4 ( mA 2 ) W
51  8.17  0.9973 
91  8.50  0.9961 
115  9.08  0.9964 
251  4.37  0.7135 
251  1.16  0.9132 
251  8.86  0.9937 
N M MSE / 10 4 ( mA 2 ) W
51  8.17  0.9973 
91  8.50  0.9961 
115  9.08  0.9964 
251  4.37  0.7135 
251  1.16  0.9132 
251  8.86  0.9937 

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 u0.

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.2 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.52 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.

FIG. 5.

Results of the analog and numerical differentiation and filtering methods from RFEA measurements showing the raw I-V characteristic with the integrated Gaussian deconvolution curve (a), and the ion-energy distribution functions against the discriminator voltage obtained with: the analog differentiator (black line) and the Gaussian deconvolution method (red line) (b), the Savitzky–Golay filter (c), the B-spline polynomial (d), the Gaussian filter (e), and the Blackman Window filter (f). The dashed curves in (b) represent the three fitted curves resulting from the Gaussian deconvolution method.

FIG. 5.

Results of the analog and numerical differentiation and filtering methods from RFEA measurements showing the raw I-V characteristic with the integrated Gaussian deconvolution curve (a), and the ion-energy distribution functions against the discriminator voltage obtained with: the analog differentiator (black line) and the Gaussian deconvolution method (red line) (b), the Savitzky–Golay filter (c), the B-spline polynomial (d), the Gaussian filter (e), and the Blackman Window filter (f). The dashed curves in (b) represent the three fitted curves resulting from the Gaussian deconvolution method.

Close modal
FIG. 6.

Result of applying the AC superimposition method to the RFEA (a) and to the RF compensated LP (b) I-V traces. The background low contrast curves are the result of attempting to compute the distributions functions using the central difference method. The output of the central difference method for the EEPF was scaled down by a factor of 100 to fit it within the plotted range.

FIG. 6.

Result of applying the AC superimposition method to the RFEA (a) and to the RF compensated LP (b) I-V traces. The background low contrast curves are the result of attempting to compute the distributions functions using the central difference method. The output of the central difference method for the EEPF was scaled down by a factor of 100 to fit it within the plotted range.

Close modal

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.

TABLE III.

Results of the different IEDF processing methods in Fig. 5(a).

Method V p (V) V B (V) MSE / 10 4 ( μ A 2 ) W
Analog  21.03  39.02  ⋯  ⋯ 
Gaussian deconv.  20.94  39.07  9.16  0.9990 
SG  20.90  39.07  8.97  0.9986 
B-spline  20.58  39.03  9.25  0.9989 
GF  21.16  38.61  9.08  0.9984 
BW  21.16  38.61  9.15  0.9984 
AC sup.  21.24  38.9  ⋯  ⋯ 
Method V p (V) V B (V) MSE / 10 4 ( μ A 2 ) W
Analog  21.03  39.02  ⋯  ⋯ 
Gaussian deconv.  20.94  39.07  9.16  0.9990 
SG  20.90  39.07  8.97  0.9986 
B-spline  20.58  39.03  9.25  0.9989 
GF  21.16  38.61  9.08  0.9984 
BW  21.16  38.61  9.15  0.9984 
AC sup.  21.24  38.9  ⋯  ⋯ 

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.53 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 u01 varying from 1.5 to 8 V and u02 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.

FIG. 7.

Comparison of the analog and numerical methods for obtaining the EEPF, showing the raw I-V characteristic (a), the EEPF obtained from the analog differentiator (b), and from the different numerical methods (c)–(f). The red dashed lines are the linear fit used to extract the electron temperature T e , fit and to perform the extrapolation to 0 eV.

FIG. 7.

Comparison of the analog and numerical methods for obtaining the EEPF, showing the raw I-V characteristic (a), the EEPF obtained from the analog differentiator (b), and from the different numerical methods (c)–(f). The red dashed lines are the linear fit used to extract the electron temperature T e , fit and to perform the extrapolation to 0 eV.

Close modal
TABLE IV.

Set of parameters for the processing used of the I-V characteristics in Figs. 5 and 7.

Method Window size, M Degree, N Sigma, σ Knots, k
EEPF IEDF EEPF IEDF EEPF IEDF EEPF IEDF
Savitzky–Golay filter  115  301  ⋯  ⋯  ⋯  ⋯ 
B-spline fitting  ⋯  ⋯  ⋯  ⋯  42  15 
Gaussian filter  ⋯  ⋯  ⋯  ⋯  25  24  ⋯  ⋯ 
Blackman window  156  140  ⋯  ⋯  ⋯  ⋯  ⋯  ⋯ 
Method Window size, M Degree, N Sigma, σ Knots, k
EEPF IEDF EEPF IEDF EEPF IEDF EEPF IEDF
Savitzky–Golay filter  115  301  ⋯  ⋯  ⋯  ⋯ 
B-spline fitting  ⋯  ⋯  ⋯  ⋯  42  15 
Gaussian filter  ⋯  ⋯  ⋯  ⋯  25  24  ⋯  ⋯ 
Blackman window  156  140  ⋯  ⋯  ⋯  ⋯  ⋯  ⋯ 

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

TABLE V.

Results of the different EEPF processing methods in Fig. 7.

Method V p(V) T eff(eV) n e / 10 11 ( cm 3 ) MSE / 10 3 ( mA 2 ) W DR (dB) δ ε e(eV)
Analog  68.08  5.65  1.87  ⋯  ⋯  55  3.2 
SG  67.14  5.52  1.85  0.91  0.9964  32  3.4 
B-spline  66.95  5.21  1.76  0.84  0.997  36  2.3 
GF  67.01  5.97  1.85  3.3  0.8015  33  3.4 
BW  67.08  5.38  1.91  3.1  0.8537  33  3.5 
AC sup.  66.95  5.43  1.83  ⋯  ⋯  23  3.3 
Method V p(V) T eff(eV) n e / 10 11 ( cm 3 ) MSE / 10 3 ( mA 2 ) W DR (dB) δ ε e(eV)
Analog  68.08  5.65  1.87  ⋯  ⋯  55  3.2 
SG  67.14  5.52  1.85  0.91  0.9964  32  3.4 
B-spline  66.95  5.21  1.76  0.84  0.997  36  2.3 
GF  67.01  5.97  1.85  3.3  0.8015  33  3.4 
BW  67.08  5.38  1.91  3.1  0.8537  33  3.5 
AC sup.  66.95  5.43  1.83  ⋯  ⋯  23  3.3 

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.13 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 Chung37 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 Godyak13 and Zhang et al.,66 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.

FIG. 8.

Analog and numerical methods for obtaining the EEPF, showing the raw I-V characteristic (a), the EEPF obtained from the analog differentiator (b), and from the different numerical methods (c)–(f). The data were acquired at z = 70 and r = 1.8 cm.

FIG. 8.

Analog and numerical methods for obtaining the EEPF, showing the raw I-V characteristic (a), the EEPF obtained from the analog differentiator (b), and from the different numerical methods (c)–(f). The data were acquired at z = 70 and r = 1.8 cm.

Close modal
TABLE VI.

Results of the different EEPF processing methods in Fig. 8.

Method V p(V) T eff(eV) n e / 10 11 ( cm 3 ) MSE / 10 3 ( mA 2 ) W DR (dB) δ ε e(eV)
Analog  70.1  6.86  1.92  ⋯  ⋯  50  4.3 
SG  69.17  7.03  2.8  0.91  0.999  33  5.2 
B-spline  69.1  6.03  2.7  0.9989  0.997  37  4.2 
GF  68.96  6.61  2.29  6.2  0.9141  32  4.8 
BW  68.96  6.9  1.67  4.0  0.9801  30  4.8 
Method V p(V) T eff(eV) n e / 10 11 ( cm 3 ) MSE / 10 3 ( mA 2 ) W DR (dB) δ ε e(eV)
Analog  70.1  6.86  1.92  ⋯  ⋯  50  4.3 
SG  69.17  7.03  2.8  0.91  0.999  33  5.2 
B-spline  69.1  6.03  2.7  0.9989  0.997  37  4.2 
GF  68.96  6.61  2.29  6.2  0.9141  32  4.8 
BW  68.96  6.9  1.67  4.0  0.9801  30  4.8 

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.

TABLE VII.

Set of parameters for the processing used of the I-V characteristic in Fig. 8.

Method Window size, M Degree, N Sigma, σ Knots, k
Savitzky–Golay filter  155  ⋯  ⋯ 
B-spline fitting  ⋯  ⋯  17 
Gaussian filter  ⋯  ⋯  30  ⋯ 
Blackman window  168  ⋯  ⋯  ⋯ 
Method Window size, M Degree, N Sigma, σ Knots, k
Savitzky–Golay filter  155  ⋯  ⋯ 
B-spline fitting  ⋯  ⋯  17 
Gaussian filter  ⋯  ⋯  30  ⋯ 
Blackman window  168  ⋯  ⋯  ⋯ 
TABLE VIII.

Effective electron temperatures and densities of the different methods after correction for the low-energy gap.

Method T eff (eV) n e / 10 11 ( cm 3 )
Analog  5.3  2.57 
Savitzky–Golay filter  5.41  2.44 
B-spline fitting  5.36  2.27 
Gaussian filter  5.36  2.27 
Blackman window  5.36  2.27 
Method T eff (eV) n e / 10 11 ( cm 3 )
Analog  5.3  2.57 
Savitzky–Golay filter  5.41  2.44 
B-spline fitting  5.36  2.27 
Gaussian filter  5.36  2.27 
Blackman window  5.36  2.27 

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.

1.
V.
Godyak
,
R.
Piejak
, and
B.
Alexandrovich
, “
Probe diagnostics of non-Maxwellian plasmas
,”
J. Appl. Phys.
73
,
3657
3663
(
1993
).
2.
C.
Charles
and
R.
Boswell
, “
Laboratory evidence of a supersonic ion beam generated by a current-free ‘helicon’ double-layer
,”
Phys. Plasmas
11
,
1706
1714
(
2004
).
3.
K.
Takahashi
,
C.
Charles
,
R.
Boswell
,
W.
Cox
, and
R.
Hatakeyama
, “
Transport of energetic electrons in a magnetically expanding helicon double layer plasma
,”
Appl. Phys. Lett.
94
,
191503
(
2009
).
4.
R. W.
Boswell
,
K.
Takahashi
,
C.
Charles
, and
I. D.
Kaganovich
, “
Non-local electron energy probability function in a plasma expanding along a magnetic nozzle
,”
Front. Phys.
3
,
14
(
2015
).
5.
S.
Ingram
and
N. S. J.
Braithwaite
, “
Ion and electron energy analysis at a surface in an RF discharge
,”
J. Phys. D
21
,
1496
(
1988
).
6.
L.
Habl
,
D.
Rafalskyi
, and
T.
Lafleur
, “
Ion beam diagnostic for the assessment of miniaturized electric propulsion systems
,”
Rev. Sci. Instrum.
91
,
093501
(
2020
).
7.
M. A.
Lieberman
and
A. J.
Lichtenberg
,
Principles of Plasma Discharges and Materials Processing
(
John Wiley & Sons
,
New York
,
2005
), pp.
165
206
.
8.
V.
Godyak
,
R.
Piejak
, and
B.
Alexandrovich
, “
Evolution of the electron-energy-distribution function during RF discharge transition to the high-voltage mode
,”
Phys. Rev. Lett.
68
,
40
(
1992
).
9.
K.
Takahashi
,
C.
Charles
,
R. W.
Boswell
, and
A.
Ando
, “
Thermodynamic analogy for electrons interacting with a magnetic nozzle
,”
Phys. Rev. Lett.
125
,
165001
(
2020
).
10.
H. M.
Mott-Smith
and
I.
Langmuir
, “
The theory of collectors in gaseous discharges
,”
Phys. Rev.
28
,
727
(
1926
).
11.
R. L.
Merlino
, “
Understanding Langmuir probe current-voltage characteristics
,”
Am. J. Phys.
75
,
1078
1085
(
2007
).
12.
V.
Godyak
and
V.
Demidov
, “
Probe measurements of electron-energy distributions in plasmas: What can we measure and how can we achieve reliable results?
,”
J. Phys. D
44
,
233001
(
2011
).
13.
V.
Godyak
, “
RF discharge diagnostics: Some problems and their resolution
,”
J. Appl. Phys.
129
,
041101
(
2021
).
14.
M. J.
Druyvesteyn
, “
Der Niedervoltbogen
,”
Z. Phys.
64
,
781
798
(
1930
).
15.
C.
Böhm
and
J.
Perrin
, “
Retarding–field analyzer for measurements of ion energy distributions and secondary electron emission coefficients in low–pressure radio frequency discharges
,”
Rev. Sci. Instrum.
64
,
31
44
(
1993
).
16.
C.
Charles
,
A. W.
Degeling
,
T. E.
Sheridan
,
J. H.
Harris
,
M. A.
Lieberman
, and
R. W.
Boswell
, “
Absolute measurements and modeling of radio frequency electric fields using a retarding field energy analyzer
,”
Phys. Plasmas
7
,
5232
5241
(
2000
).
17.
W.
Cox
,
C.
Charles
,
R. W.
Boswell
, and
R.
Hawkins
, “
Spatial retarding field energy analyzer measurements downstream of a helicon double layer plasma
,”
Appl. Phys. Lett.
93
,
071505
(
2008
).
18.
K.
Takahashi
,
Y.
Shida
, and
T.
Fujiwara
, “
Magnetic-field-induced enhancement of ion beam energy in a magnetically expanding plasma using permanent magnets
,”
Plasma Sources Sci. Technol.
19
,
025004
(
2010
).
19.
A.
Bennet
,
C.
Charles
, and
R.
Boswell
, “
In situ electrostatic characterisation of ion beams in the region of ion acceleration
,”
Phys. Plasmas
25
,
023516
(
2018
).
20.
R.
Imai
and
K.
Takahashi
, “
Deflections of dynamic momentum flux and electron diamagnetic thrust in a magnetically steered RF plasma thruster
,”
J. Phys. D
55
,
135201
(
2022
).
21.
A. M.
Keesee
,
E. E.
Scime
,
C.
Charles
,
A.
Meige
, and
R.
Boswell
, “
The ion velocity distribution function in a current-free double layer
,”
Phys. Plasmas
12
,
093502
(
2005
).
22.
A.
Savitzky
and
M. J. E.
Golay
, “
Smoothing and differentiation of data by simplified least squares procedures
,”
Anal. Chem.
36
,
1627
1639
(
1964
).
23.
I. D.
Sudit
and
R. C.
Woods
, “
A workstation based Langmuir probe system for low-pressure DC plasmas
,”
Rev. Sci. Instrum.
64
,
2440
2448
(
1993
).
24.
S.
Bétchu
,
A.
Soum-Glaude
,
A.
Bès
,
A.
Lacoste
,
P.
Svarnas
,
S.
Aleiferis
,
A.
Ivanov
, Jr.
, and
M.
Bacal
, “
Multi-dipolar microwave plasmas and their application to negative ion production
,”
Phys. Plasmas
20
,
101601
(
2013
).
25.
L.
Conde
,
J. L.
Domenech-Garret
,
J. M.
Donoso
,
J.
Damba
,
S. P.
Tierno
,
E.
Alamillo-Gamboa
, and
M. A.
Castillo
, “
Supersonic plasma beams with controlled speed generated by the alternative low power hybrid ion engine (ALPHIE) for space propulsion
,”
Phys. Plasmas
24
,
123514
(
2017
).
26.
J.
Damba
,
P.
Argente
,
P.
Maldonado
,
A.
Cervone
,
J.-L.
Domenech-Garret
, and
L.
Conde
, “
Multiprobe characterization of plasma flows for space propulsion
,”
J. Phys.: Conf. Ser.
958
,
012002
(
2018
).
27.
N.
Gulbrandsen
and
Å.
Fredriksen
, “
RFEA measurements of high-energy electrons in a helicon plasma device with expanding magnetic field
,”
Front. Phys.
5
,
2
(
2017
).
28.
J.
Fernández Palop
,
J.
Ballesteros
,
V.
Colomer
, and
M.
Hernández
, “
A new smoothing method for obtaining the electron energy distribution function in plasmas by the numerical differentiation of the I-V probe characteristic
,”
Rev. Sci. Instrum.
66
,
4625
4636
(
1995
).
29.
F.
Magnus
and
J. T.
Gudmundsson
, “
Digital smoothing of the Langmuir probe I-V characteristic
,”
Rev. Sci. Instrum.
79
,
073503
(
2008
).
30.
H.-J.
Roh
,
N.-K.
Kim
,
S.
Ryu
,
S.
Park
,
S.-H.
Lee
,
S.-R.
Huh
, and
G.-H.
Kim
, “
Determination of electron energy probability function in low-temperature plasmas from current - Voltage characteristics of two Langmuir probes filtered by Savitzky–Golay and Blackman window methods
,”
Curr. Appl. Phys.
15
,
1173
1183
(
2015
).
31.
R.
Sloane
and
E.
MacGregor
, “
An alternating current method for collector analysis of discharge-tubes
,”
London, Edinburgh, Dublin Philos. Mag. J. Sci.
18
,
193
207
(
1934
).
32.
C.
Chung
,
S.
Kim
, and
H.-Y.
Chang
, “
Electron cyclotron resonance in a weakly magnetized radio-frequency inductive discharge
,”
Phys. Rev. Lett.
88
,
095002
(
2002
).
33.
B.
Crowley
,
D.
Homfray
,
S.
Cox
,
D.
Boilson
,
H.
de Esch
, and
R.
Hemsworth
, “
Measurement of the electron energy distribution function by a Langmuir probe in an ITER-like hydrogen negative ion source
,”
Nucl. Fusion
46
,
S307
(
2006
).
34.
J. Y.
Bang
,
A.
Kim
, and
C. W.
Chung
, “
Improved measurement method for electron energy distribution functions with high accuracy and reliability
,”
Phys. Plasmas
17
,
064502
(
2010
).
35.
J.-L.
Jauberteau
and
I.
Jauberteau
, “
Determination of the electron energy distribution function in the plasma by means of numerical simulations of multiple harmonic components on a Langmuir probe characteristic–measurements in expanding microwave plasma
,”
Meas. Sci. Technol.
18
,
1235
(
2007
).
36.
R.
Hatakeyama
,
N.
Sato
,
Y.
Tsunoda
,
H.
Sugai
, and
Y.
Hatta
, “
Ion-energy distribution in a plasma under a diverging magnetic field
,”
J. Appl. Phys.
45
,
85
88
(
1974
).
37.
J.-Y.
Bang
and
C.-W.
Chung
, “
A numerical method for determining highly precise electron energy distribution functions from Langmuir probe characteristics
,”
Phys. Plasmas
17
,
123506
(
2010
).
38.
K. F.
Schoenberg
, “
Pulsed electrostatic probes as a diagnostic for transient plasmas
,”
Rev. Sci. Instrum.
49
,
1377
1383
(
1978
).
39.
V.
Godyak
,
R.
Piejak
, and
B.
Alexandrovich
, “
Measurement of electron energy distribution in low-pressure RF discharges
,”
Plasma Sources Sci. Technol.
1
,
36
(
1992
).
40.
K.
Takahashi
,
C.
Charles
,
R.
Boswell
,
T.
Kaneko
, and
R.
Hatakeyama
, “
Measurement of the energy distribution of trapped and free electrons in a current-free double layer
,”
Phys. Plasmas
14
,
114503
(
2007
).
41.
K.
Takahashi
,
C.
Charles
,
R.
Boswell
,
M. A.
Lieberman
, and
R.
Hatakeyama
, “
Characterization of the temperature of free electrons diffusing from a magnetically expanding current-free double layer plasma
,”
J. Phys. D
43
,
162001
(
2010
).
42.
S.
Yadav
,
S.
Ghosh
,
S.
Bose
,
K.
Barada
,
R.
Pal
, and
P.
Chattopadhyay
, “
Role of ion magnetization in formation of radial density profile in magnetically expanding plasma produced by helicon antenna
,”
Phys. Plasmas
25
,
043518
(
2018
).
43.
G.
Celant
and
M.
Broniatowski
,
Interpolation and Extrapolation Optimal Designs 1: Polynomial Regression and Approximation Theory
(
John Wiley & Sons
,
New York
,
2016
), pp.
57
62
.
44.
A.
Perperoglou
,
W.
Sauerbrei
,
M.
Abrahamowicz
, and
M.
Schmid
, “
A review of spline function procedures in R
,”
BMC Med. Res. Methodol.
19
,
1
16
(
2019
).
45.
K. G.
Xu
and
M. L.
Walker
, “
Plume characterization of an ion-focusing hall thruster
,”
J. Propul. Power
28
,
1105
1115
(
2012
).
46.
J.
Luo
,
K.
Ying
,
P.
He
, and
J.
Bai
, “
Properties of Savitzky–Golay digital differentiators
,”
Digital Signal Process.
15
,
122
136
(
2005
).
47.
W. H.
Press
and
S. A.
Teukolsky
, “
Savitzky–Golay smoothing filters
,”
Comput. Phys.
4
,
669
672
(
1990
).
48.
C.
De Boor
and
C.
De Boor
,
A Practical Guide to Splines
(
Springer-Verlag
,
New York
,
2001
), Vol.
27
.
49.
R. G.
Lyons
,
Understanding Digital Signal Processing
,
3rd ed.
(
Pearson Education
,
2011
).
50.
H. C.
Hayden
, “
Data smoothing routine
,”
Comput. Phys.
1
,
74
75
(
1987
).
51.
R.
Chassaing
,
Digital Signal Processing and Applications with the C6713 and C6416 DSK
(
Wiley-Interscience
,
Hoboken, NJ
,
2005
).
52.
A.
Caldarelli
,
F.
Filleul
,
C.
Charles
,
R.
Boswell
,
N.
Rattenbury
, and
J.
Cater
, “
Radial characterization of an ion beam in a deflected magnetic nozzle
,”
J. Electric Propul.
1
,
10
(
2022
).
53.
F.
Filleul
,
A.
Caldarelli
,
C.
Charles
,
R.
Boswell
,
N.
Rattenbury
, and
J.
Cater
, “
Characterization of a new variable magnetic field linear plasma device
,”
Phys. Plasmas
28
,
123514
(
2021
).
54.
F.
Filleul
,
A.
Caldarelli
,
R.
Boswell
,
C.
Charles
,
N.
Rattenbury
, and
J.
Cater
, “
The role of ion magnetization on plasma generation in a magnetic nozzle RF device
,”
J. Electric Propul.
1
,
20
(
2022
).
55.
W.
Cox
,
C.
Charles
,
R. W.
Boswell
,
R.
Laine
, and
M.
Perren
, “
Magnetic ion beam deflection in the helicon double-layer thruster
,”
J. Propul. Power
26
,
1045
1052
(
2010
).
56.
A.
Caldarelli
,
F.
Filleul
,
C.
Charles
,
N.
Rattenbury
, and
J.
Cater
, “
Preliminary measurements of a magnetic steering system for RF plasma thruster applications
,” AIAA Paper No. 2021-3401,
2021
.
57.
G. D.
Conway
,
A. J.
Perry
, and
R. W.
Boswell
, “
Evolution of ion and electron energy distributions in pulsed helicon plasma discharges
,”
Plasma Sources Sci. Technol.
7
,
337
347
(
1998
).
58.
I. D.
Sudit
and
F. F.
Chen
, “
RF compensated probes for high-density discharges
,”
Plasma Sources Sci. Technol.
3
,
162
(
1994
).
59.
A. K.
Sen
and
A.
Sen
,
Regression Analysis: Theory, Methods and Applications
(
Springer-Verlag
,
New York
,
1990
).
60.
S. S.
Shapiro
and
M. B.
Wilk
, “
An analysis of variance test for normality (complete samples)
,”
Biometrika
52
,
591
611
(
1965
).
61.
M. M.
Rahman
and
Z.
Govindarajulu
, “
A modification of the test of Shapiro and Wilk for normality
,”
J. Appl. Stat.
24
,
219
236
(
1997
).
62.
V.
Puech
and
L.
Torchin
, “
Collision cross sections and electron swarm parameters in argon
,”
J. Phys. D
19
,
2309
(
1986
).
63.
A.
Bennet
,
C.
Charles
, and
R.
Boswell
, “
Non-local plasma generation in a magnetic nozzle
,”
Phys. Plasmas
26
,
072107
(
2019
).
64.
L.
Tsendin
, “
Electron kinetics in glows–from Langmuir to the present
,”
Plasma Sources Sci. Technol.
18
,
014020
(
2009
).
65.
I.
Kaganovich
,
V.
Demidov
,
S.
Adams
, and
Y.
Raitses
, “
Non-local collisionless and collisional electron transport in low-temperature plasma
,”
Plasma Phys. Controlled Fusion
51
,
124003
(
2009
).
66.
Y.
Zhang
,
C.
Charles
, and
R. W.
Boswell
, “
Density measurements in low pressure, weakly magnetized, RF plasmas: Experimental verification of the sheath expansion effect
,”
Front. Phys.
5
,
27
(
2017
).

Supplementary Material