On the applicability of three and four parameter fits for analysis of swept embedded Langmuir probes in magnetised plasma

The ion collection of a protruding LP corresponding to those in the divertor of COMPASS was simulated using the PIC method, which can properly resolve finite ion Larmor effects and non-quasineutral plasma within the Debye sheath. The geometry of the divertor probes is inherently three dimensional and so the preferred tool for its simulation is the 3D3V PIC code SPICE3 [7, 10, 11]. However, due to the size of the probe and the shallow angle of incidence of the magnetic field with respect to divertor tile surface, a full-scale calculation with realistic plasma parameters is unfortunately too computationally demanding. Instead, a subset of simulations with feasible plasma parameters was performed and compared with simulations of the 2D3V PIC code SPICE2 [12, 13] (2D simulations are in general computationally less demanding), which assumes infinite poloidal length of the probe. The differences between 2D and 3D simulations were characterised and finally the 2D simulations were performed for realistic plasma parameters and geometry.

The geometry of the simulations is depicted in figure 1. Particles are injected through the top boundary, the side boundaries are periodic and the bottom boundary represents the wall (for details of this scheme see [13]). The probe geometry is based on the LPs utilised in the combined LP and ball-pen probe array in COMPASS divertor [9]. The probe has a triangular cross-section in the toroidal direction, with base length L_p = 6.5 mm and height H_p = 1.5 mm. The poloidal cross-section is rectangular, with width W_p = 2.0 mm. The poloidal spacing between neighbouring probes of the array is D_p = 2.0 mm (note that in the context of the simulations, the poloidal direction coincides with the radial direction). The dimensions of the simulated area were chosen such that in the direction normal to the wall there was enough space above the probe to resolve the magnetic pre-sheath (L_z = H_p + 10r_Li [14], where r_Li is the ion Larmor radius). In the poloidal direction, the width of the simulation is simply L_x = W_p + D_p. The toroidal dimension was chosen to be long enough to capture the magnetic shadow cast by the probe (see later in this section). Given the fact that the boundaries at x = 0, x = L_x , y = 0 and y = L_y are periodic, this setup corresponds to a simulation of an infinite probe array or infinite probe matrix in case of 2D and 3D simulations respectively. During the course of the simulation, the particles were introduced through the top boundary and gradually filled an initially empty simulated volume. Once the simulation achieved saturation, the normalised potential of the probe Φ = (V_p − V_plasma)/kT_e was increased in a step-wise manner from −10 to +10 by increments of 0.1. The current during each step was averaged to form a data point on the I–V characteristics. First 50% of time steps during each voltage step were discarded in order to avoid possible interference of the recorded current by the reaction of the simulated particles to the sudden change of probe bias. Note that the duration of the sweep in the simulation is not representative of the experimental conditions: the entire simulation represents ∼1 μs, while the sweeping period in the experiment was ∼1 ms. The potential at the upper boundary was set to 0, which corresponded to the plasma potential. The potential of the surrounding divertor tiles (blue in figure 1) was fixed at −3 kT_e, which was an approximate value of the floating potential.

Zoom In Zoom Out Reset image size

Since our model involved only a small volume of plasma in the vicinity of the probe head, it could not capture long range effects which the probe may have induced in the plasma. In particular, when the probe had large positive bias (with respect to the floating potential), it may have drained electrons from the adjacent flux tube, which would result in a reduction of the electron probe current (so-called flux tube charging) [15, 16]. In order to avoid this effect, we did not consider electron current in this part of the I–V in further analysis.

The first series of simulations represented a scan in the angle of incidence of the magnetic field (ψ in figure 1) for computationally affordable plasma conditions. The angle was scanned in the range from 10° down to 1.25°, where the lower bound is representative of conditions in the COMPASS outer divertor, where the angle of incidence is typically between 3° (far SOL) and ∼0.5° (strike point). The choice of plasma parameters for this initial scan was not motivated by a specific discharge in COMPASS but by the desire to compare both 2D and 3D simulations, while capturing the scales which were expected to be important for probe current collection.

Initially, it was assumed that the driving mechanism of the sheath expansion effect was the polarisation drift in the magnetic pre-sheath (as described in [17] for the case of flush-mounted probes), which was a consequence of a strong gradient of E field in this region. As such, the magnitude of the drift was assumed to scale with the thickness of this pre-sheath (which in turn scales with r_Li [14]). The sheath expansion effect was expected to be more pronounced when the majority of the probe body was submerged in the magnetic pre-sheath, therefore we selected T_i = T_e = 40 eV, which is representative of plasma conditions in the vicinity of the outer strike point [18]. The density, on the other hand, was expected to have only minor impact on the sheath expansion as it influences the E field only in the Debye sheath via the Debye length ${\lambda }_{\text{D}}\sim {n}_{\text{e}}^{-1/2}$. As such, it was chosen to significantly reduce the density to speed up the calculations. Note that results presented in section 2.5 have contradicted this initial assumption and confirmed that in the particular case of this probe geometry and plasma conditions, the role of Debye sheath and plasma density is significant, which is consistent with the view presented in [19].

The cell size in a PIC grid is required to be smaller than λ_D, so a low density scenario (with large λ_D) allows to cover the simulated volume by a smaller number of cells. Therefore, the density n_e = 10¹⁷ m⁻³ was selected for the initial set of simulations, while the plasma density in the vicinity of the outer strike point in COMPASS is in the range of 10¹⁸ m⁻³.

As mentioned before, the aim of the first series of simulations was to compare the results of a simplified 2D simulation with the full 3D simulation as it was essential to cross-check whether the two codes agree in an identical simulation geometry. This was not just an exercise in algorithm benchmarking since SPICE3 calculates electric field in all three dimensions, which enables e.g. E × B drift in the y–z plane (which is not included in the 2D simulation). In order to achieve such a comparison, a dedicated 3D simulation was performed with W_p equal to L_x , which resulted in a simulation of an infinitely wide probe (just like in the 2D simulation). The resulting ion current (normalized to its value at Φ = −3) is shown in figure 2(A) for ψ = 10°. A very good agreement for the infinite width of the probe is achieved between 2D and 3D simulations, confirming that both codes are well compatible. However, as expected for the 3D finite width of the probe, there is a notable difference with the 2D simulation, which can be explained by the presence of magnetic pre-sheath in x direction (between the probes of the array in the poloidal direction). Regarding the electron current (figure 2(B)), all cases match (same slope), with the infinite probe having a greater current only due to its greater collection area.

Zoom In Zoom Out Reset image size

The scan was extended by further simulations in order to assess the effect of the mass ratio μ = m_i/m_e. In order to accelerate the simulations, it is common to employ a reduced value of μ = 200. This setup generally should not affect the potential profile in the sheath nor the ion current. In some of the previous simulations, it was however observed that it can affect the electron probe current. Simulations with realistic μ = 3670 (corresponding to deuterium plasma) were therefore performed. Fortunately, in the case of this particular setup the role of reduced μ is negligible both on ion (figure 2(A)) and electron (figure 2(B)) currents. Nevertheless, we will employ simulations with realistic μ for the analysis of the electron current whenever possible.

Another important check is the length of the simulation in the toroidal direction (y axis). At grazing angle of incidence ψ the probe casts a long magnetic shadow and if the simulation is not sufficiently large, the probe may start shadowing itself. As shown in figure 2(C), the ion current is insensitive to the length for L_y ⩾ 1.5 L_shadow, where L_shadow = H_p/tan ψ is the length of the magnetic shadow cast by the probe. All the following simulations were designed to fulfil this criteria.

2.3.1. Density, temperature and B field incidence

The variation of the ion current with the angle of incidence and density for both 2D and 3D simulations is shown in figure 3. The presence of the sheath expansion effect is clearly manifested by a slope on the ion current for V_p < V_float and is most pronounced for cases with low density and low angle of incidence. There is a notable difference between 2D and 3D simulations—the 3D simulations in general exhibit stronger sheath expansion (as it was already observed in section 2.2). Note that the variation of ion current is not limited to bias below the floating potential: for Φ > −3 the ion current is strongly reduced. In such case, Φ is higher than the potential of neighbouring tiles (which are set to −3 kT_e in the simulations) and it is the sheath expansion above these tiles which is responsible for this effect. The associated polarization drift is pushing ions onto the surface of these tiles, preventing them from reaching the probe. When performing analysis of the total probe current obtained in experiment, the reduction of the ion current produces the same effect as if the electron current was increased—therefore acting against the reduction of electron current due to e.g. flux tube charging. As such, the reduced ion current contributes to the formation of the local minimum of the T_e,fit(V_cut,upper) dependence (which is discussed in section 3.1).

Zoom In Zoom Out Reset image size

The values of the ion saturation current obtained in these simulations can be compared with its theoretical value calculated as

$$\begin{equation}{I}_{\text{sat,}\;\text{theor}}=e{n}_{\text{s.e}}{A}_{\text{eff}}\sqrt{\frac{{T}_{\text{e}}+{T}_{\text{i}}}{{m}_{\text{i}}}},\end{equation} \tag{ 3 }$$

where T_i is the ion temperature, n_s.e. is the electron density at the sheath entrance, m_i the ion mass and A_eff is the effective probe collecting area for ions, which was assumed in previous works to be equal to 2.8 mm², based on the geometrical cross-section of the probe [9] (the real probes have rounded edges, and as such slightly reduced cross-section). Note that in the simulations the nominal value of I_sat is assumed to be the value of ion current at Φ = −3, which is the approximate location of the floating potential. The variation of the ion saturation current (figure 4(A)) shows a sensitivity both on angle and density. The variation with angle is, however, reduced at small angles and, especially in case of 3D simulations, the absolute value is close to I_sat,theor (assuming T_e = T_i). These observations also confirm that the value of A_eff assumed in previous work is correct.

Zoom In Zoom Out Reset image size

As a further step to characterise the sheath expansion, we introduce its normalised magnitude Δ_i as

$$\begin{equation}{{\Delta}}_{\text{i}}\enspace [{\text{kT}}_{\text{e}}^{-1}]=-\frac{1}{{I}_{\text{i}}({V}_{\text{p}}=-3\enspace {\text{kT}}_{\text{e}})\enspace [\text{A}]}\bar{\frac{\mathrm{d}{I}_{\text{i}}}{\mathrm{d}{\Phi}}}\enspace [\text{A}/{\text{kT}}_{\text{e}}],\end{equation} \tag{ 4 }$$

where $\bar{\frac{\mathrm{d}{I}_{\text{i}}}{\mathrm{d}{\Phi}}}$ is the slope of the ion current for Φ < −3 (obtained using linear fit of the current in this bias range, thus assuming β = 1 in equation (2)). This parameter essentially describes how large is the relative increment of the ion current per 1 kT_e when the bias decreases below the floating potential. Thus, it presents a convenient normalisation of the parameter α_i introduced in equation (2), which can be expressed as:

$$\begin{equation}{\alpha }_{\text{i}}\enspace [\text{A}/\text{V}]=-\frac{{{\Delta}}_{\text{i}}{I}_{\text{i}}({V}_{\text{p}}=-3\enspace {\text{kT}}_{\text{e}})}{{T}_{\text{e}}}.\end{equation} \tag{ 5 }$$

Figure 4(B) shows the variation of Δ_i in the performed simulations. Despite original expectations, the variation with plasma density is significant. However, the difference in the magnitude of Δ_i between 2D and 3D simulations can be reconciled using a fixed constant pre-factor 1.65, allowing us to employ 2D simulations to study the variation of Δ_i at high densities.

In addition to the simulations with varying density and angle of incidence, a scan in T_e was executed for n_e = 1 × 10¹⁸ m⁻³ and ψ = 3°. The resulting scaling (shown in figure 8(A)) is well characterised by a square root dependence.

2.3.2.  B field strength

2.3.3. Probe geometry

So far all the simulations targeted the geometry of the LPs in the combined probe array. For completeness, we have extended the study to include geometry of the swept probe array, which was installed in COMPASS already during its operation in Culham [20] and was retained after re-installation in Prague. Originally, the probes in the swept array had an approximately circular cross-sections protruding 2.0 mm above the divertor tile with a probe length equal to 7.0 mm and probe width 3.0 mm. It was suspected that after many years of operation, their geometry may have been affected by the repeated impact of plasma particles, reducing their height. This was confirmed by a post-mortem analysis using an optical 3D microscope (Leica DVM6, Leica Microsystems, Switzerland), revealing a range of probe heights between 0.5 mm (in the vicinity of the strike points) and 1.75 mm at the HFS, as shown in figure 5(A).

Zoom In Zoom Out Reset image size

In order to assess whether results obtained in section 2 are applicable also to these probes, an additional set of 2D simulations was performed for circular probes extending 2.0 mm, 1.5 mm and 0.5 mm above the divertor tile respectively for realistic plasma conditions (T_e = 40 eV, n_e = 5 × 10¹⁸ m⁻³), as shown in figure 6. The results for these plasma conditions show a strong insensitivity to the particular shape of the probe for fixed H_p = 1.5 mm. However, Δ_i is sensitive to the probe height, which was expected. In addition, the absolute values of the ion current show that A_eff of the swept probe array is 7.15, 4.68 and 2.49 mm², respectively (for H_p = 2.0, 1.5 and 0.5 mm), including the factor 1.4 derived from the comparison between 2D and 3D simulations in section 2.3.1. A linear interpolation was used to derive collecting areas for each LP, as shown in figure 5(B). Note that these collecting areas are for most of the probes lower than the value 6.52 mm² used in the experimental analysis [21], which has a direct impact on the deduced ion current density.

Zoom In Zoom Out Reset image size

So far we have assumed that the potential of the neighbouring divertor tiles is close to the floating potential—meaning that the SOL currents passing through the tiles can be neglected. However, such currents are routinely observed in COMPASS [22] and it is common for the floating potential of divertor LPs to be distinctly different from the machine ground (which is the potential of the divertor tiles). Since the sheath expansion of the probe is in competition with the sheath expansion of the neighbouring tiles, modification of tile potential could produce a visible effect on Δ_i. Additional 2D simulations for selected subset of plasma parameters were performed, in which the potential of the tiles V_T was scanned. Note that simulations with V_T > V_float can suffer from formation of an artificial sheath near the injection plane when large fraction of the impacting electron current is absorbed by the tile. For this reason the scan was limited to V_T ⩽ V_float + 1 kT_e, where the magnitude of such artificial sheath was insignificant (below 0.1 kT_e).

Figure 7 summarises the results. In general, when the potential of the tiles is above the floating potential (as measured by the probe), Δ_i is reduced with respect to the ambipolar case. This is mainly due to an increase of ion current to the probe, which is used as normalisation parameter in equation (4). In the first order approximation, this trend can be characterised by a linear function:

$$\begin{align}\hfill {G}_{\text{nonamb}}({V}_{\text{T}}-{V}_{\text{float}})& ={{\Delta}}_{\text{i}}/{{\Delta}}_{\text{i}}({V}_{\text{T}}={V}_{\text{float}})\hfill \\ \hfill & =-1.68({V}_{\text{T}}-{V}_{\text{float}})+2.25.\hfill \end{align} \tag{ 6 }$$

Zoom In Zoom Out Reset image size

In case of a divertor tile which is ∼1 kT_e above the local floating potential (common case in the vicinity of the outer strike point), Δ_i is reduced by $\sim 20\%$.

Combining all the simulations corresponding to the geometry of LPs in the combined probe array together, we derived an empirical scaling for Δ_i:

$$\begin{align}\hfill {{\Delta}}_{\text{i,}\;\text{model}}& =0.0041\times {G}_{\text{nonamb}}\times {\left({T}_{\text{e}}\enspace [\text{eV}]\right)}^{1/2}\hfill \\ \hfill & \quad \times {\left({n}_{\text{e}}[\times 1{0}^{18}\enspace {\text{m}}^{-3}]\right)}^{-1/3}{\mathrm{sin}}^{-1/3}\,\psi .\hfill \end{align} \tag{ 7 }$$

This scaling provides a decent fit to the data obtained in simulations, as demonstrated in figure 8(B) (only simulations with V_T = −3 kT_e shown here). It predicts an increase of sheath expansion magnitude for small angles of incidence, however the overall magnitude of Δ_i is suppressed for high density and/or low temperature.

Zoom In Zoom Out Reset image size

Note that unlike in earlier works on this topic (e.g. [4, 19]) scaling parameters and their functional dependencies in equation (7) were not guided by simplified models of the sheath but derived solely empirically, using the output of simulations. This approach allowed us to avoid possible biases due to imperfections of such simplified models and was largely enabled thanks to the advancements of high performance computing, which allowed us to perform series of simulations for realistic plasma conditions. Nevertheless e.g. the scaling with angle of incidence ∼ sin^−1/3  ψ is reminiscent of the dependence sin^−1/2  ψ, which appeared also in model for sheath expansion parameter for the flush-mounted probe by Bergmann [4], which was confirmed by Gunn [17]. More detailed discussion on this parameter can be found in [7]. The square root dependence on T_e is a surprising result since this means that α_i is practically independent on T_e (since ${\alpha }_{\text{i}}\sim 1/\sqrt{{T}_{\text{e}}}$). The derived scaling is far from being completely generic, for example the pre-factor 0.0041 was derived for the particular shape of the probe used in COMPASS. It was outside the scope of this work to characterise the influence of all possible probe shapes which are used in tokamaks.

The scaling of Δ_i with T_e and n_e is somewhat unfortunate, because it depends on quantities which should be obtained by the fit of I–V characteristic itself. A solution of this problem can be the use of tangent estimate method, which is described in section 3.2.

Equation (7) suggests that Δ_i should be rather small for realistic plasma conditions and B field orientation in the COMPASS divertor (⩽0.1 per 1 kT_e). In order to verify whether the sheath expansion is small enough to be completely neglected in the analysis of swept I–V characteristics, it was important to assess how does its value influence the final value of T_e obtained by a fit to the total probe current (consisting of contributions of both ions and electrons). This is achieved by fitting the probe current for V_p < V_float +1 kT_e using the three-parameter fit. Unlike the analysis of the experimental data, in case of simulations the deviation of T_e,fit can be easily quantified, since the temperature of electrons which are injected in the simulation T_e,orig is known. This comparison is presented in figure 9 and confirms the initial findings by Matthews [3]—the stronger sheath expansion (larger Δ_i), the large over-estimation of temperature. The dependence can be characterised by a linear dependence and suggests that, for realistic plasma conditions, the expected over-estimation is ⩽25%.

Zoom In Zoom Out Reset image size

Dedicated PIC simulations have demonstrated that for the particular geometry of the divertor LPs and the plasma conditions relevant to the COMPASS divertor, the sheath expansion effects should not have large influence on the analysis of experimental data. While there appears to be a systematic difference in the ion current obtained in 2D and 3D simulations, this variation can be reconciled by a simple constant pre-factor, meaning that relatively inexpensive 2D simulations can be ran instead of more computationally costly 3D runs. In addition, in case of the combined probe array, for realistic angles and plasma parameters the effective ion collecting area A_eff matches quite remarkably the value assumed in experimental analysis (i.e. 2.8 mm²), although this was previously deduced using only a simple geometrical consideration. For probes in the swept probe array, the A_eff used for experimental analysis appears to be overestimated by 50% or more, depending on the degree of probe erosion.

We have developed a fully automatised routine to analyse measurements of the swept LPs by adopting the minimum T_e method [23, 24] and the stationary bootstrap [25], which is described in detail in this section. The workflow of data analysis is illustrated here for a single Ohmic low confinement discharge which is run at the beginning of each experimental day on COMPASS (I_p = −180 kA, B_T = −1.15 T, n_e = 4 × 10¹⁹ m⁻³). Based on previous observations, the area of interest was the vicinity of the outer strike point, where largest T_e and smallest angle of incidence were expected (which should result in largest sheath expansion effects, as shown in the previous section). Therefore, probe #LPA36 (R = 510 mm, located 33 mm outside the outer strike point) was used. The measurements during the discharge flat-top were analysed using the three- and four-parameter fits (assuming β = 1 in equation (2)).

For each I–V characteristic, a series of fits was performed with a varying upper boundary (V_cut,high) of the fitted voltage range: from ${V}_{\text{float}}+\frac{1}{4}\enspace {\text{kT}}_{\text{e,}\;\text{tangent}}$ up to ${V}_{\text{float}}+\frac{3}{2}\enspace {\text{kT}}_{\text{e,}\;\text{tangent}}$, where T_e,tangent was an estimate of T_e obtained by the tangent estimate method described in section 3.2. From this series of fits, the one with the minimum resulting T_e was selected, as illustrated in figure 10. This method [23] allows to maximise the precision of the fit (by including some data at V_p > V_float) and it avoids over-estimation of T_e due to distortion of electron current in the magnetised plasmas [15]. However, due to fluctuations of the probe current, some of the fits may under-estimate the real temperature and by choosing the minimum T_e one could introduce a systematic bias into the analysis. In order to mitigate this effect and also to achieve realistic errorbars of the measurement, the stationary bootstrap method [25] was incorporated into the analysis toolchain. This method is based on the construction of a large number of virtual datasets based on random portions of the original experimental data (in our case each data point consist of a pair of current and voltage measurements), fitting these datasets and obtaining the distribution of resulting T_e, I_sat and V_float. If these distributions have Gaussian shape, the precision of the fit can be approximated as the standard deviation of the distribution and the final value as its mean.

Zoom In Zoom Out Reset image size

In the stationary bootstrap method, the new datasets are constructed by choosing blocks of data at random location in the original dataset. The length of the block is drafted from an exponential distribution with the decay parameter being equal to the autocorrelation time. This time, as deduced from the measurements of ion current at negative probe voltage, is taken as ∼5 μs. New blocks are selected until the length of the new virtual dataset is equal to the length of the original dataset. Since the blocks are selected randomly, some data points may be selected multiple times, while other data points may be completely missing. This allows to assess whether some outlying data points (e.g. due to impacting filaments) are affecting the resulting T_e (especially when the minimum T_e method is applied). Typically, ∼100 datasets are generated for each I–V characteristic, while ∼50 fits are needed for the minimum T_e method. In total, ∼5000 fits are required to obtain a reliable result from a single I–V characteristic.

Figure 10 shows an example of I–V and the scan for minimum T_e using a three-parameter fit. As can be seen, sheath expansion effect is not pronounced, which is in agreement with PIC modelling for this range of plasma parameters. The semi-empirical formula introduced in equation (7) predicts Δ_i = 0.07 per kT_e (I_sat = 60 mA corresponds to density n_e = 2 × 10¹⁸ m⁻³, angle of incidence ψ ∼ 1.5°). Note that no binning of the experimental data was applied for the analysis, instead all available data points were used for the fit. Binning is sometimes used to speed up fit convergence and for visual inspection of the fidelity of the fit, however, it can introduce bias in the data in situations where fluctuations of the probe current are of non-Gaussian nature, which is a common observation with probe current [26]. In such a case, the mean value is not a good representation of the underlying distribution of measurements.

The distributions of T_e, I_sat, V_float and the upper boundary of the fit V_cut,max are shown in figure 11. Except for V_cut,max all of them have approximately Gaussian distribution of values, which allows to use their mean and standard deviation as representations of the most probable value and confidence interval. The distribution of the upper voltage boundary of the fit shows multiple peaks, which is caused by the presence of probe current fluctuations in the vicinity of the upper boundary (forming multiple local minima in figure 10(B)). The scan in number of samples constructed for bootstrapping shows that estimation of mean values and standard deviations is possible already for 100 samples, however at least 1000 samples is required for analysis of the exact shape of the distributions. Figure 12 demonstrates correlations between the fitted parameters. There appears to be a significant correlation between the obtained T_e and I_sat, which is not surprising since ${I}_{\text{sat}}\sim \sqrt{{T}_{\text{e}}}$. An important result is the weak dependence only of the fitted T_e on the choice of the upper voltage boundary for the fit. Note that the distribution of V_float in figure 11 does exhibit a bimodal distribution, however the scatter of the values is generally quite small.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the purpose of further analysis of swept probes, we introduce a new method of T_e estimation, which uses only data in the vicinity of the floating potential, which we call the tangent method. The method works as follows:

• (a)  
  The I–V data points are sorted according to voltage and then smoothed using running mean filter. The smoothing window size N_window is calculated as N_window = S_w N₁₀₀, where S_w is typically equal to 0.4 and N₁₀₀ is the number of data points with −100 $< {V}_{\text{p}}<$ 0 V. The optimum value of S_w depends on the level and character of noise of current measurements and as such needs to be tailored to a specific data acquisition system.
• (b)  
  An estimate of V_float is found by iterating through the I–V characteristic and looking for voltage where current of smoothed I–V characteristic crosses 0. In order to make this procedure robust, the algorithm looks for interception of zero from left and right side independently. The estimate of V_float is then obtained as an average of these two values. This helps to mitigate problems in case of distorted I–Vs with multiple crossings of zero.
• (c)  
  A tangent to the I–V characteristic is found at V_p = V_float. The slope is obtained by linear fit of the smoothed I–V characteristic in the interval of V_float ± 5 V. Note that the slope S of this tangent is equal to −I_sat/T_e (when sheath expansion effects are neglected).
• (d)  
  The tangent is extended towards negative bias as I_tangent. For each datapoint of the smoothed I–V (I_j , V_j ) having V_j < V_float, the following ratio is calculated:

  $$\begin{equation}r={I}_{\text{tangent}}({V}_{j})/{I}_{j}.\end{equation} \tag{ 8 }$$

  Such point satisfies the following equation

  $$\begin{equation}-({V}_{j}-{V}_{\text{float}})/{T}_{\text{e}}=r(1-{e}^{({V}_{j}-{V}_{\text{float}})/{T}_{\text{e}}}).\end{equation} \tag{ 9 }$$

  This equation is transcendental but a solution for T_e can be expressed using the Lambert function W [28]:

  $$\begin{equation}{T}_{\text{e},j}=-({V}_{j}-{V}_{\text{float}})/(-r-W(-r/{e}^{r})).\end{equation} \tag{ 10 }$$

  In practice the calculation of Lambert function is performed using a fast algorithm developed by Fukushima [29]. For example, for a data point with corresponding tangent current being twice the value of the smoothed measured current (r = 2) the solution reads T_e = −(V_j − V_float)/1.59. The final T_e,tangent estimate is obtained as a mean value of all obtained values of T_e with 0.6 < r < 3. Extending the procedure towards higher values of r would bring risk of influence of the sheath-expansion effects (which are small in the vicinity of the floating potential since these scale as ∼(V − V_float)).

The method is illustrated in figure 13(A). The magenta dashed line shows the tangent to the smoothed I–V, the vertical dashed red line shows obtained value of V_j . The reliability of this method was tested on the dataset of swept LPs in the combined probe array as shown in figure 13(B) (including measurements with large sweeping range). While there is noticeable scatter in the data, the method works reasonably well for a wide range of temperatures, especially for T_e < 40 eV. As the tangent method does not require any iterative fitting, it is notably faster than the method described in section 3.1 and appears to be a good candidate for real-time applications. Note that the tangent method is distinctly different from the double probe technique, which can also make use of the first derivative of the I–V characteristic [2].

Zoom In Zoom Out Reset image size

It is instructive to use synthetic probe data to inspect this method. One important question is its sensitivity to sheath expansion—since this method does not use current measurements at large negative voltage, it should be insensitive to Δ_i. We applied two approaches when dealing with synthetic probe data: (i) we generated smooth probe current for a pre-set of constant values of T_e, I_sat, V_float and α_i (labelled synthetic) and (ii) we used fast measurements of the combined probe array and employed time-dependent values of T_e, I_sat, V_float (using data from probe triplet #35 in #20474). The value of Δ_i was calculated using equation (7), where ψ was artificially manipulated to obtain a scan in sheath expansion. In both cases, the probe current was generated for a synthetic voltage waveform, which had the same range and frequency of sweeping as in the experiment.

The resulting dependence of T_e,tangent on Δ_i is shown in figure 14(A). Both approaches yield very similar trends, while using the time-dependent values of plasma parameters resulted in significantly larger and asymmetric errorbars (here represented by 10% and 90% quantiles of dataset measured over a flat top of #20474). The method is reasonably stable for Δ_i < 0.2. As such, it provides an independent mean to cross-check the validity of three- or four-parameter fits (which is presented in next sections).

Zoom In Zoom Out Reset image size

Another important check is to evaluate the stability of the resulting T_e,tangent with respect to the small errors in determination of the slope of the tangent S = dI_p/dV_p(V_p = V_float) and V_float, which are shown in figures 14(B) and (C). In both cases, we have introduced a random error in the process of its determination (with Gaussian probability function) and observed the effect on T_e,tangent. This to some extent mimics the real-world case where the experimental data are acquired with some level of noise and fluctuations. In case of the slope S, the error in its determination gradually increases the scatter of T_e,tangent. However, it does not introduce a significant bias in the results. Error in 10% of determination of S yields $\sim 20\%$ scatter in T_e,tangent. In case of V_float the situation is more complicated, and having δV_float/T_e > 0.2 leads to an over-estimation of the obtained temperature estimate. However, note that the scale of the scan is rather extreme compared to the properties of the experimental data. Fitted I–V characteristics typically show match of V_float with precision of ∼1 V, which for T_e ∼ 40 eV yields δV_float/T_e ∼ 0.05.

3.3.1. Scan in V_cut,low

The generous range of swept voltage applied to the probe allowed us to assess the role of the voltage range V_range = −(V_min − V_float)/T_e on precision of fits of the probe data, especially with respect to the measured T_e. This can be done easily by setting different limits of the lower voltage boundary V_cut,low in the analysis, thus excluding parts of the I–V characteristics from the fitting procedure. Figure 15 shows an example of results obtained with both three-parameter and four-parameter fits from the I–V characteristic presented in figure 10. It can be seen that for the largest available voltage range the results for three- and four-parameter fits are very similar. As the voltage range is reduced, the three-parameter fit shows good stability in terms of T_e,fit for −320 $< \ {V}_{\text{cut,}\;\text{low}}\ <$ 200 V, which suggests that the fitting is valid. On the other hand, the four-parameter fit starts to yield lower and lower temperatures. For V_cut,low ⩾ −250 V (V_range ∼ 5 kT_e) the resulting temperatures are in strong disagreement with the values obtained with maximum available voltage range. This bias towards low temperatures can be explained by an interference between the two fitting parameters: T_e and α_i. In vicinity of V_float, the exponential term in equation (2) can be linearized and for β ∼ 1 the fit becomes ill-defined:

$$\begin{equation}{I}_{\text{p}}\sim -\frac{{I}_{\text{sat,}+}}{{T}_{\text{e}}}({V}_{\text{p}}-{V}_{\text{float}})+{\alpha }_{\text{i}}({V}_{\text{p}}-{V}_{\text{float}}).\end{equation} \tag{ 11 }$$

Zoom In Zoom Out Reset image size

Both the second and the third term scales linearly with V_p − V_float, it is impossible to reliably determine both T_e and α_i. This problem is demonstrated in figure 15(B), which shows the dependence of α_i on V_cut,low (and thus on available voltage range). In principle, its value should be constant (and close to zero in this respective case), however, there is a clear trend of increased magnitude of α_i with reduced voltage range.

As a possible way to circumvent this problem, the four-parameter fit can be applied in two steps: (i) the slope of the linear term is determined using data obtained at the lowest probe bias and (ii) the extra current due to sheath expansion is subtracted from the I–V characteristic, which can be then treated using the three-parameter fit. We have attempted to employ this two-step technique (green points in figure 15) but it did not yield more reliable results. On the contrary, this method is prone to fluctuations in the probe current (which can hamper determination of the linear term), which results in larger scatter of the resulting parameters. These fluctuations can be caused e.g. by MHD modes with characteristic frequency close to the probe sweeping frequency.

The two-step fitting technique can be significantly improved if the coefficient α_i is not determined from the slope of the I–V characteristic but if it is calculated using equations (7) and (5). This is demonstrated by the yellow points in figure 15. For the full range of swept voltage, the prediction of the semi-empirical law yields larger α_i than obtained by the four-parameter fit, which leads to a certain under-estimation of T_e in this situation (26.5 eV vs 37 eV given by the three-parameter fit). However, the agreement improves when the voltage range is reduced and the match is generally better than for a regular four-parameter fit. Since the prediction of α_i is based on the tangent method, it is rather insensitive to V_range (as shown in figure 15(B)).

In case of the three-parameter fit, the reduction of obtained temperatures is much less pronounced, however, having V_cut,low > −200 V (V_range ∼ 4 kT_e) also leads to a reduction of the fitted temperature. This effect may be caused by the I–V characteristic not strictly following an exponential dependence, as predicted e.g. by Rozhanski [27].

A practical counter-measure against using the four-parameter fit in case of insufficient voltage range is to rely on the value of V_range and apply a four-parameter fit only if V_range is above a certain threshold [7]. However, the above analysis shows that this technique will fail if one uses T_e originating from the four-parameter fit itself. When the range is not sufficient, T_e will be under-estimated, leading to an over-estimation of V_range. The same problem, albeit less pronounced, can also occur for the three-parameter fit (e.g. V_range based on T_e from three-parameter fit in figure 15(C) never drops below 4). The tangent method allows to circumvent this problem, since it is not using the part of I–V characteristic most affected by the sheath expansion.

3.3.2. Criteria for application of the fits

The boundaries for safe application of I–V characteristic fitting presented in figure 15 may be to some extent accidental, demonstrated only for a single particular case. In order to achieve more robust criteria, we have repeated the scan in V_range over a larger data set measured during diverted phases of selected L-mode and ELM-free H-mode discharges (each data point represents a value obtained using a 10 ms time window). In each case we noted at which value of V_range does T_e drop by 20 and 50 percent, respectively, with respect to the value obtained at full voltage range. The resulting histograms are presented in figure 16. Note that we have used T_e,tangent to determine V_range, as it provides an independent reference point which only weakly depends on the actual voltage range. Vertical dashed lines indicate median values of the distributions. It can be seen that indeed the four-parameter fit requires at least 4 kT_e of voltage range to achieve a systematic error better than 20%. If the voltage range is less then 3 kT_e, then the resulting temperature is under-estimated by factor of two or more. For the three-parameter fit, these boundaries are 2.9 and 1.4 kT_e, respectively. Figure 16(C) demonstrates that when prediction of the sheath expansion based on PIC simulations (equation (7)) is used in the two step technique, the reliability of the fitting procedure is, to large extent, recovered: V_range of at least 2.5 kT_e is needed to achieve precision of the fit of 20% and for V_range smaller than 1.6 kT_e the resulting temperature is under-estimated by at least 50%.

Zoom In Zoom Out Reset image size

Finally, figure 16(D) shows the results of the same analysis for the tangent method, which requires V_range ⩾ 1.7 kT_e and ⩾1.1 kT_e to achieve precision of 20% and 50% respectively. The good stability of this method with respect to reduced V_range is caused by the fact that it uses data points in a close vicinity of the floating potential.

The usefulness of V_range derived using the tangent method is further demonstrated when α_i obtained by four-parameter fits are compared to the semi-empirical scaling formula equation (7). In particular, the sheath expansion parameter α_i obtained using the four-parameter fit can be confronted with the prediction of Δ_i presented in equation (7) (here the parameters needed for prediction of Δ_i were obtained using the tangent method). Such comparison (as shown in figure 17(A)) demonstrates that the semi-empirical scaling serves as good estimate of experimental Δ_i for most of the measurements. However, measurements with V_range < 4 kT_e show strong deviation from this trend. The same group of points exhibits disagreement between four-parameter fit and the tangent method, as shown in figure 17(B). This indicates that for small voltage range, the four-parameter fit yields incorrect results—over-estimated α_i and under-estimated T_e. This is the result of interference between fitting parameters described earlier.

Zoom In Zoom Out Reset image size

During regular operation, the combined probe array determines T_e using the difference in floating potentials of LP and ball-pen probe measurements as

$$\begin{equation}{T}_{\text{e,BPP}}=\frac{{V}_{\text{float,BPP}}-{V}_{\text{float,LP}}}{\alpha },\end{equation} \tag{ 12 }$$

where α is assumed to have a fixed value of 1.4 [9]. This measurement remained available even in the experiments reported here, which provided a straightforward way to cross-check the results obtained using swept probes.

Figure 18 shows results of such a comparison for our dataset of swept probe measurements. Data points with V_range > 3 kT_e are plotted in colour, measurements with insufficient V_range are plotted as grey points (with grey outline).

Zoom In Zoom Out Reset image size

The agreement between the two methods is satisfactory in most cases. Peak temperatures at the outer target in the range of 50 eV are consistently observed by both methods (see an example of such measurements in figure  19(c)), which also agrees with older measurements performed at COMPASS-D [20]. This indicates that the LFS SOL is in the sheath-limited regime with almost constant T_e along the magnetic field lines. Further analysis of the outliers revealed that when the floating potential was positive—either in L-mode (discharges #12416 and #12401, represented by stars) or following a transition into an ELM-free H-mode (#20504 and #20505, represented by diamonds), the LP + BPP method yields very low (figure 19 (b)) or even negative (figure 19 (a)) temperatures, while the temperatures derived from the swept LP remain above 20 eV. The detailed cause of this behaviour is not understood at this moment and will be a subject of future investigations. Note that it is often observed that swept LPs fail to deliver temperatures below a certain threshold (∼5 eV) even if other diagnostics suggest colder plasmas [30, 31]. Unfortunately, we have not yet devised a method to determine this threshold value of T_e.