Data Processing Techniques for Ion and Electron Energy Distribution Functions

The data used in this study has been acquired using Moa Caldarelli et al. (2022), a radio-frequency plasma reactor that is an expansion of Huia Filleul et al. (2021, 2022). It consists of a $150\leavevmode\nobreak\$cm long, $9\leavevmode\nobreak\$cm inner diameter borosilicate glass tube plasma source connected to a $70\leavevmode\nobreak\$cm long, $50\leavevmode\nobreak\$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)\leavevmode\nobreak\$cm. A base pressure of ${\sim}10^{-7}\leavevmode\nobreak\$Torr is routinely achieved, while the argon working pressure ranges between $0.5\times 10^{-3}\leavevmode\nobreak\$Torr and $1\times 10^{-3}\leavevmode\nobreak\$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\leavevmode\nobreak\$kW RF generator working at $27.12\leavevmode\nobreak\$MHz. A static magnetic field is applied with a pair of movable Helmholtz solenoids, placed concentrically with the glass tube. The coils position along the plasma source length can be adjusted to obtain the desired magnetic field topology.

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\leavevmode\nobreak\$mm, and a $0.1\leavevmode\nobreak\$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 Charles et al. (2000); Charles and Boswell (2004); Cox et al. (2010); Bennet, Charles, and Boswell (2018), and a detailed description of the RFEA and its electric circuit used herein are presented in Ref. Caldarelli et al., 2021. In order to map the energy of the ions, the voltage of the discriminator grid is swept from $0\leavevmode\nobreak\$V to $80\leavevmode\nobreak\$V to obtain the $I$-$V$ characteristics. 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 $\pm 40^{\circ}$ are obtained Charles and Boswell (2004). 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 Conway, Perry, and Boswell (1998); Charles et al. (2000); Charles and Boswell (2004) 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. $\lambda_{\rm i}>>s,d$ (where $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 $\lambda_{\rm e}$ Godyak and Demidov (2011). This ensures that the probe only disturbs the plasma in ways that can be accounted for by theory. Following recommendations in Ref. Godyak, 2021, Langmuir’s theory is satisfied if the probe tip radius $r_{\rm p}$, the tip length $l_{\rm p}$, the tip holder radius $r_{\rm h}$ and the local Debye length $\lambda_{\rm D}$ all satisfy the condition $r_{\rm p}\ln{\left(\pi l_{\rm p}/4r_{\rm p}\right)},r_{\rm h},\lambda_{\rm D}\ll\lambda_{\rm 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 be much smaller than the electron Larmor radius Godyak and Demidov (2011). To meet these conditions, the LP used in this work is made from a tungsten tip with $r_{\rm p}=25\leavevmode\nobreak\ \mu$m and $l_{\rm p}=3\leavevmode\nobreak\$mm, mounted at the extremity of a glass pipette tapering down to $r_{\rm h}=0.75\leavevmode\nobreak\$mm.

Measuring an EEDF requires biasing the LP above the local plasma potential and can typically draw electron currents in 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. Godyak and Demidov, 2011 for more details). For this purpose a grounded aluminium sleeve of adequate surface area was inserted inside the dielectric plasma source tube at the end opposite to the diffusion chamber Filleul et al. (2021).

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 Sudit and Chen (1994); Godyak and Demidov (2011). While this effect is not easily visible in the probe $I$-$V$ characteristic, it becomes evident with the second derivative and, without proper mitigation, it severly compromises EEDF studies. Because the root-mean-square of $V_{\rm p}$ at the first and second harmonics are larger than $T_{\rm e}/2$ in the present apparatus, RF compensation was added to the Langmuir probe, following design steps described in detail in Sudit and Chen (1994) and Godyak and Demidov (2011). 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\leavevmode\nobreak\$MHz and $54.24\leavevmode\nobreak\$MHz with shunt capacitors of approriate 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. 1, 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\leavevmode\nobreak\$dB at the first and second harmonics provides sufficient filter impedance to suppress the $V_{\rm 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 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. The maximum probe biasing voltage was adjusted to $V_{\rm bias}\simeq V_{\rm p}+5\leavevmode\nobreak\$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_{\rm bias}$. Probe cleaning by ion bombardment was used when necessary.

The Langmuir probe is mounted on a movable doglegged eccentric dielectric shaft which can span the entire length of the plasma source tube Filleul et al. (2021). 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 current is deduced by measuring the voltage drop across a resistor with an isolation amplifier. The output of the amplifier is then split and fed into both a data acquisition system (DAQ) for digitization, and into an analog differentiator which outputs the first and second time derivatives. The LP $V_{\rm bias}$ and the RFEA $V_{\rm D}$ are defined by several periods of a triangle wave function. The data acquisition is triggered by the output on the function generator to ensure signal synchronicity and the voltages are digitized at $4,000$ samples per period. One LP measurement is the average of $100$ $V_{\rm bias}$ sweeps, while one RFEA measurement averages $200$ $V_{\rm D}$ sweeps.

The use of analog differentiator circuits predates numerical differentiation methods, and they are still regularly employed for obtaining EEDFs and occasionally for IEDFs Godyak, Piejak, and Alexandrovich (1992b); Takahashi et al. (2009); Takahashi, Shida, and Fujiwara (2010); Boswell et al. (2015); Yadav et al. (2018). The active analog differentiator used in this study, designed to work with a sweeping frequency of $10\leavevmode\nobreak\$Hz, follows the design described in Takahashi et al. (2010). Two cascading operational amplifiers (Renesas CA3140) are each tuned to resonate at ${\small{\sim}}800\leavevmode\nobreak\$Hz by adjusting the resistors and capacitors values. This results in a linear gain frequency increase of $20\leavevmode\nobreak\$dB per decade up to the resonance frequency, i.e. input signals of frequencies up to ${\small{\sim}}800\leavevmode\nobreak\$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 ${\small{\sim}}80\leavevmode\nobreak\$kHz. Therefore the analog differentiator configuration also works as a band-pass filter, attenuating high frequency noise. Since real-world circuits deviate from the ideal-case, the linearity was checked to be strictly valid up to at least ${\small{\sim}}150\leavevmode\nobreak\$Hz. Each stage of the analog differentiator gives an output proportional to the first and second time derivative of the probe $I(t)$ signal. Knowing the biasing/discriminator voltage sweep period and voltage swing, the differentiator output can be multiplied by $\mathop{}\!\mathrm{d}t/\mathop{}\!\mathrm{d}V$ to retrieve the $\mathop{}\!\mathrm{d}I_{\rm p}/\mathop{}\!\mathrm{d}V$ derivative. Differentiator circuits have the advantage of producing repeatable outputs since no tuning of parameters for different experimental conditions are 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 Magnus and Gudmundsson (2008). Since the $I$-$V$ characteristics 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$ the polynomial degree Magnus and Gudmundsson (2008). This leads to Runge’s phenomenon and over-fitting Celant and Broniatowski (7 62), which is visible in the results presented by Magnus and Gudmundsson (2008). Another shortcoming is the non-locality of the polynomial fit, i.e. distant data points directly impact the local fit Perperoglou et al. (2019). The SG filter and the B-spline methods both circumvent these problems by fitting low degree polynomials to subsets of the data.

The Gaussian deconvolution method, or Gaussian fitting, has been used to obtain the ion-energy distribution function from the $I$-$V$ characteristics measured by a retarding field energy analyzer Charles and Boswell (2004); Cox et al. (2008); Imai and Takahashi (2022). 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 minimises 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 curve

where the parameters to be fit are: the amplitude $a$, the curve width $c$, and the location of the Gaussian peak $V_{\rm 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 Charles and Boswell (2004); Bennet, Charles, and Boswell (2018). 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 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$ characteristics that could lead to large errors in the inferred plasma parameters and loss of information.

An analysis similar to that used by Fernández Palop et al. (1995) and Magnus and Gudmundsson (2008) is used to assess the numerical methods in this work. It is noted that in Ref. Fernández Palop et al., 1995 and Ref. Magnus and Gudmundsson, 2008, the comparison of the different techniques was conducted on simulated $I$-$V$ characteristics, while in this study only experimental data are used. The following parameters, evaluated between the raw probe collected current $I_{\rm p}$ and the processed data $I_{\rm f}$, were considered: the mean squared error (MSE) and the Shapiro-Wilk (SW) test statistic $W$.

When comparing different regression models, 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_{\rm p}-I_{\rm 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 Sen and Sen (1990). 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 Shapiro and Wilk (1965); Rahman and Govindarajulu (1997). 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 distorts the data, one should aim to maximize $W$ and re-evaluate 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$ Rahman and Govindarajulu (1997). 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 peak of the energy distributions, while achieving enough noise reduction when evaluating $\mathop{}\!\mathrm{d}I_{\rm C}/\mathop{}\!\mathrm{d}V_{\rm D}$. An excessive smoothing of the $I$-$V$ characteristics 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 2 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 1. 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_{\rm 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. 2 (a)). By extension, $N=4$ is the lowest B-spline degree that can be used to obtain a smooth second derivative for the EEPF. 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. 2 (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 $\delta_{\varepsilon_{\rm e}}$ which is defined as the $\varepsilon_{\rm e}$ value of the maximum of the EEDF. Because the bulk of the electrons are in the low-energy region of the EEPF, the larger $\delta_{\varepsilon_{\rm e}}$, the larger the errors in determining $n_{\rm e}$ and $T_{\rm eff}$ from Eq. 3 and Eq. 4. An EEDF is considered to have sufficient energy resolution if $\delta_{\varepsilon_{\rm e}}\leq(0.3-0.5)T_{\rm e}$ Godyak and Demidov (2011); Godyak (2021). The dynamic range of the EEPF is defined as the ratio of its maximum to its minimum resolvable value. A dynamic range $\geq 60\leavevmode\nobreak\$dB (equivalently of 3 orders of magnitude and above) is necessary to resolve high-energy electrons in the inelastic range of the EEPF, which are detrimental in investigating excitation, ionization, and transport processes Godyak and Demidov (2011); Godyak (2021). An EEDF is often separated into its elastic and inelastic ranges, defined as $\varepsilon_{\rm e}<\varepsilon_{\rm e}^{*}$ and $\varepsilon_{\rm e}>\varepsilon_{\rm e}^{*}$, respectively. $\varepsilon_{\rm e}^{*}$ is the lowest excitation energy of the working gas, e.g. 11.55  eV for argon Puech and Torchin (1986). Minimizing $\delta_{\varepsilon_{\rm e}}$ and maximizing the DR is therefore paramount in the process of determining the optimal numerical parameters, together with minimizing the MSE and maximizing the value of $W$.

Figure 3 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 2. In Fig. 3 (a)-(d) $N=2$, and 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  dB to 46  dB owing to the decreasing cut-off frequency. However, this is accompanied by distortions of the signal, visible from the worsening of the low-energy resolution $\delta_{\varepsilon_{\rm 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 Fig. 3 (d)-(f) and Table 2. Regarding the Shapiro-Wilk test, all the cases with $M=251$ failed the normality test with $W<0.996$. Overall, trading-off between $\delta_{\varepsilon_{\rm e}}$, the DR, the MSE and $W$, the best performing filter for this 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 3 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 Roh et al. (2015).

As mentioned in Section I, the ion-energy distribution function can be derived from the first derivative of the measurements acquired with a retarding field energy analyzer. The respective $I$-$V$ characteristics was taken in the expansion chamber with the RFEA orifice located 19  cm downstream of the magnetic nozzle ($z_{\rm P}=19$ cm). The experimental conditions were maintained at an RF power of $250\leavevmode\nobreak\$W, argon pressure in the chamber of $0.5\leavevmode\nobreak\$mTorr, and a maximum magnetic field on axis of $B_{z\rm,max}=330\leavevmode\nobreak\$G. Under these conditions, the plasma density peaks axially under the solenoids reaching a value of $\approx\leavevmode\nobreak\ 10^{12}\leavevmode\nobreak\ \rm cm^{-3}$ which then decreases to ${\sim}\leavevmode\nobreak\ 10^{10}\leavevmode\nobreak\ \rm cm^{-3}$ in the expansion chamber. The data processing of the collected current $I_{\rm C}(V_{\rm D})$ measured with the energy analyzer and the subsequent computation of the ion-energy distribution function followed the approach described in Section IV.1.

Table 4 summarizes the results of the parametric analysis that provide the smoothest IEDF possible for the numerical smoothing methods analyzed, while keeping intact the peak features. Figure 4 shows the ion-energy distribution functions obtained with the different analog and numerical techniques from experimental data. 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 data in Fig. 4 has a signal-to-noise ratio of approximately $66\leavevmode\nobreak\$dB. The raw $I$-$V$ curve, shown in Fig. 4 (a), was obtained by averaging 200 consecutive probe sweeps and it is plotted with the fitted Gaussian deconvolution to illustrate one of the data processing methods. As seen in Fig. 4, 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) which location can be considered the local plasma potential $V_{\rm p}$, while the second peak corresponds to the accelerated ion beam population with potential $V_{\rm B}$ end energy $\varepsilon_{\rm B}=e(V_{\rm B}-V_{\rm 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 Caldarelli et al. (2022). It is also noted that the ion-neutral mean free path for the operating pressure of $0.5\leavevmode\nobreak\$mTorr is shorter than the distance between the source exit and the probe position i.e., $\lambda_{\rm i}{\sim}10$ cm and $z_{\rm p}=19$ cm. Thus, the accelerated ion beam component is decreased 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 3 presents the key results obtained with the different data processing techniques, namely: the local plasma potential $V_{\rm p}$, the ion beam potential $V_{\rm 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 deviations for $V_{\rm p}$ and $V_{\rm B}$ of $\pm 0.21\,$V and $\pm 0.23\,$V, respectively. For the collected data, all the smoothing methods perform well in terms of mean squared error and the SW test statistic. The MSEs are $\leq 9.25\times 10^{-4}\leavevmode\nobreak\ \mu 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 methods 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.

The RF compensated LP data were acquired with the apparatus operated at $200\leavevmode\nobreak\$W of RF power, $1\leavevmode\nobreak\$mTorr of argon, $B_{z,\rm max}=300\leavevmode\nobreak\$G and with the Helmholtz solenoids placed 30  cm away from the loop antenna. These conditions create an inductively coupled plasma exhibiting a symmetrical single peaked axial plasma density gradient ranging from $10^{12}\leavevmode\nobreak\ \rm cm^{-3}$ under the solenoids to $10^{10}\leavevmode\nobreak\ \rm cm^{-3}$ at the extremities of the glass tube Filleul et al. (2021, 2022); Bennet, Charles, and Boswell (2019). The LP was placed on axis halfway between the antenna and the solenoids. The processing of $I_{\rm p}(V_{\rm bias})$ to obtain the EEPF with the analog and numerical methods was conducted following the steps described in Section IV.1. Table 4 gives the resulting optimal numerical parameters.

Figure 5 (a) shows the raw $I$-$V$ LP characteristic together with the optimized Savitkzy-Golay filter output from Fig. 3, i.e. for $(N=2,M=115)$. As reference, the signal-to-noise ratio of the LP characteristic is $48\leavevmode\nobreak\$dB. In Fig. 5 (b)-(f), the EEPFs obtained from the various methods and the deduced plasma parameters are shown. $T_{\rm eff}$ and $n_{\rm e}$ are calculated by integrating the EEPF according to Eq. 4 and Eq. 3, respectively. Their values are reported together with $V_{\rm p}$, the MSE, $W$, dynamic range and energy resolution for each method in Table 5. The maximum variance in the obtained values of $V_{\rm p}$ is less than 2$\%$. Likewise, for $T_{\rm eff}$ and $n_{\rm e}$, their maximum variances are ${\sim}13\%$ and ${\sim}8\%$, respectively. Thus, in terms of estimating the plasma parameters, all 5 methods have comparable performance, i.e. within the absolute error tolerances usually associated with LP measurements Boswell et al. (2015); Godyak (2021).

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. 5 (d) resulted in the best energy resolution at $\delta_{\varepsilon_{\rm e}}\simeq 2.3\leavevmode\nobreak\$eV and is the only method to satisfy the $\delta_{\varepsilon_{\rm e}}\leq(0.3-0.5)T_{\rm e}$ condition. The second-best performing method is the analog differentiation with $\delta_{\varepsilon_{\rm e}}\simeq 3.2\leavevmode\nobreak\$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 LP implementation, in particular of the plasma-wall sheath impedance, or by the probe tip surface condition Godyak (2021). Techniques to further improve the resolution are available in Ref. Roh et al., 2015; Godyak, 2021.

The analog differentiation method in Fig. 5 (b) resulted in an EEPF with a significantly larger dynamic range of approximately 3 orders of magnitude, compared to 2 orders of magnitude for the numerical methods. The analog differentiation is therefore the only method presented which can reliably resolve inelastic processes and can provide confidence on the nature of the electron distribution function, i.e. a single-Maxwellian for $\varepsilon_{\rm e}$ up to 30  eV in the present case. This is further highlighted in Fig. 5 (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 artificial oscillations in the case of the B-spline would make it difficult to reach conclusions on the nature of the electron distribution function. Among the numerical methods, the Savitzky-Golay and the B-spline are the best performing for this data set, with the lowest distortion of the original $I$-$V$ curve, as shown by their lower mean-square errors and Shapiro-Wilk test statistics in Table 5. 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.

Finally, it is interesting to note that following the guidelines in Ref. Godyak, 2021, the low-energy gap $\delta_{\varepsilon_{\rm e}}$ impacts can be corrected by extrapolating the EEPF from its linear portion down to $\varepsilon_{\rm e}=0\leavevmode\nobreak\$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 $\varepsilon_{\rm e}^{3/2}$, if the EEPF is Maxwellian for some $\varepsilon_{\rm e}$ greater than the energy of EEPF peak, then the EEPF is likely to be Maxwellian at lower energies Tsendin (2009); Kaganovich et al. (2009); Godyak (2021). This extrapolation can be achieved with the linear regressions shown in Fig. 5 as dashed red curves. Applying this process to the analog EEPF in Fig. 5 (b) for example, the effective electron temperature obtained with the corrected EEPF is $T_{\rm eff^{\prime}}=4.32\leavevmode\nobreak\$eV and the corrected electron density $n_{\rm e^{\prime}}=2.57\times 10^{11}\leavevmode\nobreak\ \rm cm^{-3}$, to compare with the values in Table 5. The electron temperature obtained from the slope of the linear regression is $T_{\rm e,fit}=4.2\leavevmode\nobreak\$eV, providing evidence that the low-energy extrapolation has closely recovered the original EEPF shape.