The Langmuir probe system in the Wendelstein 7-X test divertor Open Access

Wendelstein 7-X (W7-X) is an optimized stellarator built to assess the fusion relevance of the stellarator concept by demonstrating steady state operation near reactor plasma pressures and confinement.^1–4 To achieve this objective, a good understanding of the plasma edge and wall interaction physics is crucial. Langmuir probes (LPs) play an important role here as a local system capable of simultaneously measuring multiple quantities. If the LPs are integrated in the target surface, the resulting perturbation of the plasma will be small. In LPs, a bias voltage V_bias is applied to a probe tip with respect to some ground, typically the plasma vessel wall. The current I that flows in response is measured, and from many data pairs, an I–V characteristic can be constructed, from which plasma properties are derived.

In W7-X, the angle of incidence of the magnetic field on the target plates is below 3° in the region of highest heat load (“strike line”). The analysis of LP characteristics is difficult under such shallow angles of incidence.⁵ It was thus important to choose an LP design with larger angles of incidence for all magnetic configurations while simultaneously not exposing the probe surface to excessive heat loads. This was achieved by a facetted probe surface. Even with this design, it proved important to select a model for the I–V characteristic that takes into account the shallow angle of incidence. In this paper, multiple models for this response characteristic are compared to establish which one is best suited to our case. We will show that the model for this response characteristic, for which many options exist in the literature, can have a large impact on the derived plasma parameters. For the comparison, we implemented the models in a generic and expandable way. Any adaptations of—and parameter choices for—the models will be described. If additional effects on the I–V characteristic are to be considered in the future, these can be easily implemented in this framework as further models. Special attention was paid to the determination of parameter uncertainties.

We will first introduce W7-X, its divertor, and the hardware of the LPs in Sec. II. In Sec. III, we will discuss the different model equations and their physical interpretation. Section IV describes the data pre-processing, the model comparison, and its results and implications. Limitations and observations so far not explained by the models are explained herein as well. The evaluation using the optimal models identified by the above comparison was implemented and performed by the Minerva framework⁶ using Bayesian probability and is discussed in Sec. V. The LP system performed reliably throughout the entire first divertor campaign of W7-X, providing a trove of data that should be evaluated. The automation necessary to accomplish this is described in Sec. VI. Finally, we highlight some key results and possibilities opened by this work in Sec. VII.

For W7-X, steady state is defined as all physics and machine operation processes reaching equilibrium.⁷ From this equilibrium requirement follow the choices of superconducting coils and electron cyclotron resonance heating (ECRH) using continuously operating gyrotrons as the main heating system,⁸ but more importantly the need for a capable divertor. Before installation of the water cooled high heat flux divertor (HHFD),⁹ the machine operated with a massive graphite divertor identical in size and shape, the test divertor unit (TDU).¹⁰ The most commonly used magnetic field configurations of W7-X have a shallow, i.e., low shear, profile of rotational transform iota that crosses a low rational value at the edge, forming large, stable magnetic islands.¹¹ These islands are used as inherent diverted sections of the plasma that like in tokamaks are intersected by divertor plates.¹² The divertor consists of ten discontinuous divertor segments, two in each of the five machine modules.¹³ One divertor module is shown in Fig. 1. The Langmuir system described in this paper is embedded into these divertor plates and diagnoses the plasma conditions there. Only the divertor plates in one of the modules are equipped with LPs, so this paper will only refer to the upper and lower divertors (UD and LD) of that module. The plasma facing surfaces of W7-X can be grouped by the heat load placed on them during an experiment. Weakly loaded areas are covered with stainless steel panels, more exposed ones are covered with graphite tiles, and the highest loads are borne by the graphite divertor. The design of the TDU is described by Peacock et al.,¹⁴ and first results of experiments are discussed by Pedersen et al.¹⁵ 

20 probes are installed in both the UD and the LD, arranged in two lines of ten. This is due to the construction of the divertor from poloidally arranged target elements (“fingers”). In the HHFD, individual carbon fiber reinforced carbon tiles are welded to the fingers with castellations between them to accommodate thermal expansion of the surface. In contrast, the TDU fingers are single pieces of fine-grained graphite but retain the castellations that define the poloidal spacing of the probes. To reduce the sources of metallic impurities, graphite¹⁷ was also the natural material of choice for the probe tips, each of which is machined from a single piece thereof. The probe centers are separated by 2.5 cm poloidally and 5 cm toroidally. The probe tips themselves are 15 mm long toroidally and 3 mm wide poloidally and protrude from the target less than 1 mm (see Fig. 2). The probes in finger 3 are offset in parallel to those in finger 4. Due to the shallow angle of incidence and to the direction of the magnetic field, each probe in finger 3, numbered 1–10, is displaced by less than 2 mm perpendicular to the magnetic field to a corresponding probe in finger 4, numbered 11–20. In the “standard” magnetic configuration, used in the experiments shown in this paper, each probe is shadowed by the target from one side such that they do not perturb the plasma and their projected area is unambiguous.

In the second half of the first divertor operation phase (OP) 1.2b, we swept the bias voltage of only the probes in finger 3 and used those in finger 4 to acquire negative current (I_sat) or voltage (V_f) measurements at constant bias (see also Sec. IV A).

To motivate the design of the LP surfaces in contact with the plasma, we made some general considerations on qualities that can be optimized by the design: impact of sheath effects, localization of the measurement, and robustness of the probe.

Designs are constrained mainly by the heat flux to the probe and its thermal contact, which should not cause the surface temperature to exceed 1800 °C¹⁸ to prevent electron emission or erosion. The high heat loads of up to 200 MW m⁻² parallel to the field require the probe surface to be inclined relative to the magnetic field direction, similar to the target surface itself. Shallow incidence of the field onto the probes however causes complications in the analysis as is well known from other experiments.⁵ 

Experiments with tilting probe arrangements^19,20 show a nonlinear effect of decreasing incidence angle γ (see Fig. 2) on the ion to electron current ratio, deviating from the expected dependence of sin γ. This is not accounted for in standard descriptions of the I–V characteristic and can distort the derived parameters. It is thus desirable to keep the incidence angle as large as feasible.

A reduction of the relative importance of shallow incidence angle corrections can be achieved by increasing the probe size. The size of the probe is in our case constrained by the construction of the divertor, but also by the desire to have a localized measurement. Since the divertor in W7-X is not toroidally symmetric but segmented and has multiple toroidal interaction zones, a design like the rail probe system in Alcator C-Mod, intended to minimize sheath effects,²¹ is not usable.

To accommodate these design constraints, a facetted design of the probes was chosen, allowing for a larger angle of incidence than would have been achieved with flush probes, but lower heat flux than on a proud probe. A toroidally aligned wedge shape would have presented a perpendicular surface to the poloidal component of the field, necessitating a slope in the poloidal and toroidal directions (see surface normals in Fig. 2 and type 1 in Fig. 3).

Across the divertor and the extent of the probe array, the magnetic field incidence angle changes, passing through zero and forming a poloidal watershed at a location dependent on the magnetic configuration.²² Probes near this watershed may therefore receive flux from either the positive or the negative toroidal direction and were machined with four sloped surfaces (type 2 in Fig. 3).

The two or four faces of the probe tips are angled—modifying the ridges in the directions indicated by arrows in Fig. 3—such that 6° of incidence is not exceeded in any of the nine vacuum reference configurations.^23,24 This ensures acceptable heat loads and, at the same time, well defined values of the angle of incidence over the entire probe surface. Both properties would be less uniform if domed probes were used. The exact values of the facet angles and thus the detailed design differ for each probe position, depending on the local range of incidence angles.

The typical projected probe area A_proj is $O$(2 mm²). Because the probes are elongated and seen by the field at a shallow angle, the shape of their projection is a rounded sheared rectangle (see Fig. 2). While the area of this projection is strongly angle-dependent, its circumference or perimeter, important for the magnitude of sheath effects as will be discussed later (Sec. III B), varies only by <5%. For our probes, the projection circumference lies in the range of 12 mm–13 mm.

Errors in the probe area due to imprecise positioning would lead to errors in the determination of n_e as well as potentially problematic leading edges.²² To quantify errors, the distance by which each probe protrudes from the target was measured to be within 10 µm before the experimental campaign. These measurements were repeated for probes in one divertor after the campaign, showing little deviation. The perimeter and, therefore, the sheath growth estimation are dominated by the length of the probes and would not be affected by this.

Each probe is contacted with two wires to allow for independent measurement of current and voltage with the target as ground. The cables used inside the machine are needed to be vacuum compatible, resistant to ECRH stray radiation, and capable of carrying currents of several ampere to up to 200 V with minimal mutual influence (“crosstalk”) and good transmission characteristics up to 500 kHz (the chosen sampling frequency).

Below the divertor, we use a custom-built cable consisting of two twisted copper cores, each in an insulating ceramic mesh sleeve and together surrounded by a copper pipe. This pipe serves to reflect ECRH stray radiation and spreads heat, while the insulating material was chosen mainly for its heat resistance and because it does not release gas or deteriorate under vacuum conditions. Along the vessel wall inside the vacuum barrier, we use a shielded twisted-pair cable with Kapton insulation compatible with vessel temperatures of up to 150 °C¹⁸ during vacuum baking, and again in copper piping for stray-radiation protection. Beyond the vacuum barrier, a 50-core twisted-pair cable is used for the 70 m connection to the electronics outside the torus hall.

To have the flexibility of using alternative high sweeping frequency electronics, coaxial cables with better signal transmission properties and 50 Ω impedance are laid out in parallel to the conventional system outside the vacuum barrier. Switching a probe from one electronic to another is done by a patch field mounted on the vacuum vessel. This is also used to ground the return current of the probes. Potential differences therefore use this point as their reference, not the target directly next to the probes, which is however in good contact with the vessel. It was also necessary to ground the measurement cards in the electronics rack, which introduces a ground loop. The potential difference between the rack and the vessel was monitored and found to be small, V_{ground difference} < 1 V.

While the probe tips and cabling had to be specifically designed for the W7-X TDU, an existing system of measurement cards could be reused. It was previously operated with the static probe array on the Wendelstein Experiment in Greifswald für die Ausbildung (WEGA).^25,26 The system consists of eight measurement cards with eight inputs each, for which voltage or current measurement can be selected. The voltage measurement is realized with a 100/1 precision voltage divider and a factor two operational amplifier (OAmp) buffer. The current is directed across a 0.1% precision shunt resistor and the difference formed by a differential amplifier (DAmp).²⁷ Multiple shunt resistance values between 0.1 Ω and 100 Ω can be selected by a jumper on the card to improve the resolution of different probe current ranges. A fuse²⁸ rated for 1 A and 250 V is placed in-line to protect the measurement components, especially the DAmp. Data were recorded using analog-to-digital converters (ADCs)²⁹ with a sampling frequency of 500 kHz. Two issues arose with the measurement cards, related to the finite common mode rejection of the differential amplifiers and to the mutual influence of measurement channels on the same card. The procedures for correcting these are outlined in  Appendixes A and  B.

Because of the greater mobility of electrons, the current drawn for positive bias is much greater than I_sat (see Sec. III). To limit the current to the probe, the bias voltage range used was −180 V to 20 V. This range was sinusoidally swept at 500 Hz by bipolar power supplies (KEPCO 400-M) for two groups of ten probes on the upper and lower divertors, respectively. The attainable sweep frequency f_sweep and bias range were limited by the use of this power supply model. The maximal transient currents for a single probe recorded by the system were on the order of 1.5 A. Due to multiple probes being operated on the same power supply, total currents in excess of 4 A per power supply sometimes occurred during periods of strong plasma fluctuation. This caused a reduction in the provided bias voltage, which is taken into account for the I–V curves.

For the current and voltage measurements, we can determine an uncertainty by propagating the uncertainties in the individual components of the system. These are the imprecision of the shunt resistor and voltage divider, the resistance of the mode-selection jumpers, the temperature dependence of cables, the digitization error, and the error we make in our common mode and track resistance correction schemes. Of these, the resistance of the mode-selection jumper has the largest uncertainty, introducing a systematic error in the current of less than 3%. The random error and level of background noise in the measurements before each experiment begins and the plasma contacts the probes is less than 5 mA.

Instead of measuring V_bias independently at each probe, we recorded the output voltage of the power supplies and calculated the probe voltage from the probe current and the known resistance in the current path. The inductive contribution from the changing current is neglected here as it is much smaller. On several probes, we verified this method by comparing the calculated probe voltage with direct measurements, using the second wire contacting the probe. The difference was zero on average with a spread of less than 3 V due to fluctuations. The current measurements were calibrated by using 1 kΩ test resistors in place of the probes. By using a larger shunt resistor, the small calibration currents thus generated were sufficient to confirm the linearity of the differential amplifiers in their entire input range.

Models of a LP make predictions of a characteristic that relates plasma parameters to the measured current given probe properties and V_bias. We determine these plasma parameters, which are usually the quantities of interest, by varying them until the model prediction matches the measured I–V curve. If the voltage is sinusoidally swept, then a full I–V curve covering all V_bias is recorded every half-period of the sweep. There is no general model of LPs, which holds true in all conditions; however, different ones have been proposed for different regimes and assumptions.

The derivation of the most well-known model and its assumptions are outlined by Hutchinson:³⁰ 

• No magnetic field or equivalently gyro-radii ρ much greater than the probe dimension a.

• No impediment to particle motion due to collision, implying the mean-free-path length λ_mfp similar to the plasma dimension A.

• A sheath surrounding the probe, with an extent defined by the location where the Bohm criterion c = c_s holds, with thickness on the order of the Debye length λ_D that is thin compared to a and relatively independent of V_bias.

• A Maxwellian electron energy distribution.

• V_bias sufficient to repel all electrons when measuring the ion saturation current I_sat such that their energy and momentum distributions are unchanged.

• Ions are drawn from a large area or reservoir and have low energy compared to the electrons and are completely absorbed by the probe if they reach the sheath.

• A density distribution of electrons governed by a Boltzmann factor such that n_e(x) = n_e,∞ exp[eV(x)/T_e], where V is the local sheath potential and n_e,∞ is the density far from the sheath.

These assumptions constitute the simple Langmuir model with the characteristic

c_s is the ion sound speed with γ_i = 3 accounting for the adiabatic cooling of the ions in the sheath, and Z is the ion charge. T_e and T_i represent energies of electrons and ions, respectively, corresponding to the temperature times the Boltzmann constant k_B. By n_e, we always refer to the density at the sheath edge, eliminating the necessity of a Boltzmann factor relating n_{e, surface} and n_{e, sheath edge}. A_coll is the effective collection area of the probe, and J_sat is thus the saturation current density.

The most obviously violated assumption is that of no magnetic field. If particles are magnetized, that is, their motion is heavily constrained across the magnetic field but essentially free along it, the collection area of the probe is no longer its entire surface but only the projection of its area onto a plane perpendicular to the field A_proj. This is of course related to the considerations about heat flux in Sec. II as the latter is proportional to the current or flux to the probe.

The second assumption we must relax is that the sheath is sufficiently thin to not contribute substantially to the collection area and its growth therefore can be ignored. Our measured I–V characteristics show a clear growth of the ion current in the negative bias region, which can be explained by an increasing collection area due to sheath growth. The sheath thickness can be approximated with the Child–Langmuir law for space-charge limited current

with c = 0.8 and λ_D being the Debye length.³¹ For special probe types and no magnetic field, there exist exact solutions such as those derived by Laframboise,³² which generally agree well with Eq. (4) within the typical measurement uncertainty of LPs. To account for the sheath growth effect in a magnetic field and for arbitrary probe shapes however, we introduce another model.

Tsui et al.³³ propose a model in which the sheath extends out perpendicular to the surface of the probe everywhere, but because of the projection along the magnetic field, only the perimeter of A_proj contributes.³³ This leads to the intuitive and generally applicable formulation

where P_proj is the perimeter of the probe projection and l_sheath is the sheath thickness. This thickness is calculated by Eq. (4) with

Using Λ = 2.5 here effectively neglects the magnetic sheath⁵ and c = 1 simplifies the Child–Langmuir expression and accounts for reported deviations from it.^34,35 We use this model with the calculated projection areas and circumferences of our probes, subtracting that part of the perimeter from which an extension perpendicular to the field would protrude into the divertor target surface, and refer to it as the PSE model.

The Bohm criterion requires ions to flow at c_s when they enter the sheath. In a magnetized plasma impinging on an inclined surface, this condition is required for the surface normal direction, as shown by Chodura.³⁶ This requires the formation of a pre-sheath with thickness on the scale of the ion gyro-radius ρ_i. Weinlich and Carlson derive a model that accounts for this sheath and describes its effect on flush mounted probes.³⁷ The sheath thickness predicted by this is given by

Using the notation of Ref. 37, l_Sheath,WC is the sheath thickness, ψ is the field incidence angle to the surface normal, Φ labels potentials, and prime denotes derivatives with respect to the distance from the wall or probe. Note that ψ in this notation corresponds to 90° − γ in ours. Dimensionless potentials are written as ϕ. Subscripts me, De, w, and pr refer to values at the magnetic sheath edge, Debye sheath, wall around the probe, and probe itself, respectively. ρ_s is the ion Larmor radius at sound speed. It is assumed that T_e = T_i, Z = 1, ψ ≥ 80°, and ions at the sheath edge are monoenergetic with E_i = T_i; ≈ in the equations signifies the use of these simplifications.

In the double probe model that will be introduced in Sec. III D, Φ_me can be determined by numerical iteration, but in practice, it is estimated to be close to V_p in Weinlich’s analysis code. We follow this approach in using Λ = 3 with Eq. (6b). This value of Λ is justified because we investigate primarily hydrogen plasmas (see Sec. IV A) and explicitly want to take the magnetic sheath into account. Φ_$w$ is zero as the vessel defines the ground and Φ_pr = V_bias. Substituting the conditions for which the Child–Langmuir law is derived, namely, an unmagnetized plasma, cold ions, and negligible electric field at the sheath edge, its $∝V34$ behavior is recovered by a power law expansion of the above expression. Our conditions are clearly different, which leads to a deviation of the sheath models. Sheath thickness predicted by the PSE and WSE models typically differs by a factor of 2–2.5, which is shown in panel (a) of Fig. 4, strongly influencing the predicted slope of the current in the ion dominated part of the characteristic. This can be understood when considering that the density is reduced in the magnetic pre-sheath by a factor dependent on the incidence angle, thereby increasing the local Debye length. Counteracting this, an angle dependent part of the total sheath potential drop occurs across the pre-sheath. Since the total sheath potential drop is constant,^5,36 the potential difference in the Debye sheath is reduced. The first effect scales as $1/cos(ψ)$ and the second as ln[cos(ψ)] such that the net effect increases the sheath thickness. Considering the magnetic sheath has two consequences illustrated by panel (b) of Fig. 4, the bias voltage for which a sheath can be formed shifts by δV_min. There is a constant offset of sheath thickness δl_min corresponding to the magnetic presheath and through the Debye length dependent on $1/ne$⁠.

For the calculation of A_coll, we do not use the absolute values of sheath thickness shown in Fig. 4 but the difference between this value and that at the unbiased target [l_CL,WC(V_bias = 0)]. This reflects that the geometry of the system is only changed when the probe sheath extends further than the wall sheath. At positive bias, this difference is negative because the probe is in part shadowed by the sheath of the surrounding target. This is considered by Tsui et al.³³ but does not result in a negative correction there because due to their consideration of a proud probe, the probe perimeter effective for sheath expansion is larger than that shared with the wall or support. For flush mounted probes, both parts of the total perimeter are approximately equal.

In Weinlich and Carlson’s paper, A_coll of the probe is increased using l_Sheath,WC, but by expanding the sheath perpendicular to the probe surface rather than perpendicular to the magnetic field as in the PSE model. Since the difference for shallow incidence angles is small, we use the latter method of expanding the area around the perimeter but using Eq. (7a) for the sheath thickness and call this model WSE. There are some discussion and much simulation effort on the expansion of the sheath and thus current collecting area parallel to the target.^35,37–40 We follow the PSE model in extending the sheath outward from the perimeter, treating the target normal and parallel direction equally.

The calculation of the sheath thickness requires the angle of incidence onto the probe for ρ_s. This angle is not uniquely defined for the facetted probes, so we use the average, weighted for the contribution of each face’s area to the total.

One of the assumptions of the SL model was that ions and electrons could stream to the probe from a large volume. In a magnetized plasma, current should flow predominantly along the field lines. We would expect particles to reach the probe from a cone, aligned with the field and opening with an angle proportional to the relative strength of perpendicular transport. If the diffusion into this cone is insufficient to provide the current drawn by the probe, the remainder must come from a “counter electrode” where the field intersects the wall.³⁰ 

Drawxing a current from a region of the wall requires it to act as a probe itself, with the current limited by A_counterJ(V_{bias, counter}). Weinlich et al. therefore argue that we must treat a single probe in a magnetic field as a virtual asymmetric double probe^37,38 to which current is predominantly supplied from the wall. Using a double probe model allows us to fit the entire I–V characteristic and explains the low level of electron saturation current relative to the expectation of $Ie, satIi, sat≈mime$⁠. Février et al.⁴¹ consider a virtual asymmetric double probe model for this reason but describe the sheath growth by two free, not physically motivated parameters. Tsui et al.³³ also propose a double probe model using their description of sheath expansion but do not consider its applicability to single probes in magnetized plasmas. We use double probe models as an independent extension to other Langmuir characteristics, denoted by a letter D prefixed to the acronym.

For this, we introduce an additional fit parameter β describing the ratio of the areas of the probe electrodes or equivalently of the saturation currents and a general double probe characteristic (here with generic A_coll for all the above models),

The wall and target of W7-X are grounded, so the bias voltage of the return probe is zero.

Tests with a specially constructed checkerboard probe¹⁹ show that the current drawn by a probe can be accounted for when summing the return current to wall elements directly adjacent to the probe. Experiments in ASDEX Upgrade (AUG) and Alcator C-Mod both observe in accordance with this current being drawn across the field from elements of a flush mounted probe toward the biased probes.^38,42 This requires significant current to flow across the field rather than parallel to it. The observation is surprising but supported by the consideration of which resistance the current would encounter if it followed the field lines until they were intersected: It can be estimated by the parallel Spitzer resistivity that the plasma resistance alone would be responsible for the entire gradient of the I–V characteristic at V_f. If this were true, it would seriously draw into question the validity of Langmuir probe analysis assuming that the resistance is due to the sheath and contradicts other measurements of finite sheath temperature. The implications of current being drawn from the target surrounding the probe are nonetheless profound. Current to wall- or target-integrated probes depends heavily on cross-field transport, nearby probes potentially influence each other, and sheath size modifications happen in the immediate vicinity. We will discuss this further in Secs. IV A and IV D.

With the above options, we can formulate six different models, listed in Table I.

To provide inferred plasma parameters, we follow a two-step approach. In the first step, currents and voltages at the probe tips are calculated from the ADC readings by taking into account the assignment and settings of probes, measurement card channels, and ADC channels. The signals are corrected for the DAmp common mode and resistances and capacitance of cables and current paths on the amplifier cards, as described in Secs. II C and II D and in detail in  Appendixes A and  B. Corrected current and voltage signals for each probe tip are then made available in the data archive in addition to the raw ADC data. These pre-processed signals are the base for the further analysis.

The signals are split into segments containing at least one full sweep of the V_bias range (i.e., a half-period). As in the scrape-off layer (SOL) of many magnetic confinement devices, the fluctuations exhibit a broad frequency range, also covering the usual sweep frequency of 500 Hz (see the example shown in Fig. 5). This was observed in many discharges when either the bias voltage of a number of probes was not swept but kept constant at −180 V or probes were unbiased. These probes thus constantly recorded I_sat or V_f. We can therefore not expect the plasma quantities to be constant during a single (half) period of the V_bias sweep period. We shall nevertheless make the assumption that they are and further discuss the impact of fluctuations in Sec. IV B and specifically Sec. IV C.

An influence of the swept probes on the adjacent constant bias or floating probes is usually visible. This provides further evidence for the asymmetric double probe hypothesis of current being drawn across the field from the target surrounding a swept probe. An interval of the signal of a floating probe together with the voltage of the adjacent swept probe is shown in the upper panel of Fig. 5. The spectral density of this floating potential signal is shown in the lower panel of Fig. 5 together with the spectral density of a H_α signal.⁴³ The lines of sight of this signal penetrate the islands that are intersected by the target plates containing the Langmuir probes at a distance of 4 m and 6 m, respectively, along the magnetic field. No impact of the sweep frequency is visible in the H_α signal. This demonstrates that the sweeping of the probe bias does not cause a global perturbation of the SOL. We also note that the characteristics and inferred parameters of the swept probes show no systematic differences between the cases of adjacent biased or floating probes.

In the second step, we infer parameter values and uncertainties from a model. To find the most suitable model, we use the simple and well established method of $χν2$ (reduced Xi-squared) comparison.⁴⁴ This method is developed from the minimization of the squared residuals, commonly known as the least-square (LS) method, and both $χν2$ and LS make assumptions fully satisfied by linear models with uncorrelated parameters. Since our models are non-linear and, as we will see, have correlated parameters, we must turn to maximum likelihood estimation to correctly determine the values and especially the uncertainties of our parameters. Maximum likelihood is a special case of Bayesian inference if all measurement uncertainties are gaussian, which we assume to be the case. For generality and reasons we will further motivate in Sec. VI, we will nonetheless use the methods of Bayesian analysis for the parameter inference.

For both model comparison and inference, code was written such that only the physics model function changes and everything necessary to retrieve data or perform the parameter optimization is independent. In other words, the fitting code is agnostic of the model function. This allows the trivial extension to additional models. For both the LS method and the Bayesian inference, we used three or four free parameters for the single and double probe models, respectively. These were the electron temperature, floating potential, ion density, and, in the case of the virtual double probe models, β. Each parameter was set with bounds and an initial value, given in Table II.

The bounds for T_e, n_e, and V_f are set liberally to allow the solver freedom to explore the parameter space while preventing only obviously unphysical results. The lower bound for β reflects the assumed absence of negative cross-field flows; the upper bound is set to allow for no reduction of the electron saturation current in an impure hydrogen plasma. As we will see in Fig. 7, high values of β are inferred when in fact no electron current saturation is observable. The lower bounds for T_e and n_e were necessary for numerical stability.

For other unknowns of the characteristics, such as the ratio of the ion to the electron temperature and the ion charge, reasonable assumptions were made. Throughout this paper, we take T_i ≈ T_e as is expected for the machine edge⁵ and standard practice.^37,41 As for the ion charge Z, the main plasma species in OP 1.2b was almost exclusively hydrogen; however, there were significant levels of carbon impurities as well as oxygen prior to boronization.^45–47 Both line-integrated and profile measurements of Z_eff relevant for the level of bremsstrahlung were made but have not been fully evaluated yet.⁴⁸ Z_eff weighs impurities more heavily; it is not simply the ratio n_i/n_e that we would require. This being unavailable, we use Z = 1 everywhere for now and determine c_s according to Eq. (3).

To quantify the quality of the model fits, we will first present some examples and then introduce and compute $χν2$ for each fit. This requires finding the optimal parameters for each model to match the data, for which the least-square (LS) method is sufficient. To test the different models and adapt them for W7-X, the flexible python package symfit⁴⁹ was used. Equations, parameters, independent variables, and constants are logically separate and individual parameters can be defined with bounds and conditions, removing the need for checks. (Co-)Variances are calculated algebraically by evaluating the Hessian, approximated by the Jacobian around the solution. We report the LS parameter uncertainties here only for completeness and do not use them any further. The algorithm for minimization is the limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm for bound-constrained optimization (L-BFGS-B).⁵⁰ 

No weights were assigned to the values, and the higher density of points at the extremes of the bias range is thus not corrected for. The measurement uncertainty at each individual voltage and current measurement point is small compared to the fluctuations of the I–V curve (see Sec. IV D), so it was not considered in the optimization process.

A first example of LS fits is shown in Fig. 6. All panels show the same data in black, and model characteristics are drawn as colored lines. The top row [panels (a) and (b)] contains the full bias range, the middle row [(c) and (d)] a zoom on the ion current branch, and the bottom row [(e) and (f)] a zoom on the electron current branch. Characteristic functions for all models are overlaid, with single probe models in the left column [(a), (c), and (e)] and virtual double probe models in the right column [(b), (d), and (f)]. The necessity to account for sheath expansion to explain the non-saturation of the ion current can be clearly seen in panels (c) and (d), as well as the WSE model’s good match with the data. The bottom row supports the virtual double probe hypothesis, showing excellent agreement of data and characteristics. We conclude that, in this specific example, the DWSE model can best fit the I–V characteristic. This is supported quantitatively by the value of $χν2$ obtained for each model given in Table III.

The parameter values of the different models given in the caption vary significantly, especially the temperature. Trends in this variation clearly originate in the shape of the characteristics: Panel (c) of Fig. 6 shows how the transition point from electron current modeled by an exponential to the almost linear trend for the non-saturating ion current shifts further to positive V_bias going from SL to PSE to WSE. This increases the inferred growth rate of the exponential and therefore reduces the temperature. Simultaneously, it reduces the ion-saturation current at zero sheath expansion, reducing the inferred density.

The saturation of the electron current requires the single probe models to assume a shallower gradient near V_f, again reducing the growth rate and thus increasing the temperature.

For higher temperatures and lower densities, we see the differences between models much less clearly, as shown in a second example in Fig. 7. The horizontal extent of the characteristic increases with temperature, which in the fixed V_bias range means we do not reach saturated ion or electron currents. All models can fit the data equally well; $χν2$ lies in the range of 0.241–0.248. The inferred parameters however still vary strongly because the models interpret the data differently: again, the PSE and WSE models shift the transition from ion current to electron current dominated branches and can match the data assuming lower temperatures. This is illustrated by I_sat calculated for zero sheath expansion or A_coll = A_proj shown in the table for each model, which marks this transition. β cannot be determined from the data; the fitter therefore assumes the maximum value and correctly reports a large uncertainty. Note that this is matched by the small variation in the parameter values inferred by the single and double probe values. As we noted before, the parameter uncertainties of the LS fit (values in parentheses in the captions of Figs. 6 and 7) are only reported for completeness, and those for T_e quoted here are almost certainly underestimated. However, even a more accurate assessment of each model individually would not capture the uncertainty in the parameters due to the uncertainty associated with the model choice. An evaluation taking this into account is possible using Bayesian probability but beyond the scope of this paper. Experimentally, this could be resolved by a larger sweep range on the scale of multiple T_e, which might make model differences more evident.

A decision on the most suited model is clearly important given the large discrepancy in inferred parameters. Based on the result for a hot plasma alone, we cannot make this choice by a $χν2$ comparison. Occam’s razor suggests we should choose the model with the smallest amount of parameters and assumptions, which would be the SL model. The more complex models do not improve the quality of the fit. The cold plasma examples however show that the DWSE model describes the data best. Since the same physics governs both situations, the same model should be used everywhere. We will base our decision only on those cases where a difference between models is clearly visible and argue that the choice is valid also for those conditions where it is not.

To quantify which model best fits our characteristics and should be used to evaluate the LP results in the future, we have been comparing the reduced Xi-squared (⁠$χν2$⁠) of the fits. χ² is a scaling of the squared residuals by the uncertainty of the data and therefore a measure of the size of observed deviations relative to expected deviations. This is normalized by the number of degrees of freedom to allow for comparison of fits to different numbers of data points and penalize complex models prone to overfitting.⁴⁴ The residuals are generally composed of random measurement error, model error attributable to a systematically incorrect model function, and parameter choice error due to sub-optimal parameters for that model. This last term is minimized by a successful run of the fitting algorithm.

In our case, two physical phenomena affect the collected current data: the sweep of V_bias generating the I–V curves from an otherwise constant plasma on its timescale and fast fluctuations for which we do not have a model. To decide which model to use, we are interested in only the first part, so we must either separate the effects by their timescale or find a model for the fluctuations. A separation is not normally possible for our data. As mentioned before, Fig. 5 shows that there is no minimum in spectral power, which we could choose as cutoff for a low-pass filter. The frequency band between half and approximately four times f_sweep is problematic because fluctuations in this band can be neither resolved nor removed by low-pass filtering without modifying the characteristics.

To calculate a comparable $χν2$ that emphasizes differences in the model error, we therefore consider fluctuations as part of the expected deviations. The fluctuation amplitude varies with the applied bias voltage, increasing with the magnitude of current. There is also a slower time dependence on the discharge conditions, with the fluctuation amplitude increasing with heating power and reducing significantly in phases when the divertor is not loaded (strong edge radiation or detachment). While we do have an independent measurement of the fluctuations without sweeps by the probes in divertor finger 4 (see Fig. 1), these only give us information on the amplitude and power spectrum at the ion saturation and floating potential bias level (−180 V and 0 V, respectively). We would need at least one more measurement point at positive bias to be able to estimate the underlying fluctuation amplitude on the swept signals.

Lacking this, the following procedure is used to estimate the expected deviation σ: For each time point, we detrend the current values in a bin centered at this time across adjacent voltages and adjacent sweeps. This kernel covers a range of 10 V and five sweeps. We assign the standard deviation of these data to the central measurement as σ. The kernel dimension in sweeps takes account of the variation of the fluctuation level during a discharge, while the voltage dimension does the same for the aforementioned difference between ion and electron current branches.

To get a quantitative answer how well the different models perform in general rather than for a single characteristic, $χν2$ is calculated for each model and sweep of a discharge. We use experiment 20180814.007 during which the plasma density was continuously increased. This led to a transition from attached divertor conditions to the entire power being dissipated by radiation in the edge and ultimately power starvation of the plasma. It must be distinguished from divertor detachment involving neutral compression, which was not observed in this discharge. An overview of this discharge is shown in Fig. 8.

Figure 9 shows the $χν2$ distributions for ∼7000 I–V curves from that experiment, separated into two panels for the first and second phases of the discharge. This corresponds approximately to the transition to the radiated power fraction $PradPECRH>0.5$⁠. As shown before in Figs. 6 and 7, a difference in the $χν2$ values of different models is only visible in colder plasmas. Here, the DWSE model systematically performs best, supporting our conclusion from above.

It should be noted that $χν2$ itself is uncertain:^56,57 For models with non-orthogonal or covariant parameters, the effective number of degrees of freedom K is less than K = N − P used by us, where N is the number of samples and P is the number of parameters. Furthermore, what we are implicitly comparing is the expectation value of the χ² distribution. To distinguish models at the above noise level, we must take into account its standard deviation of $σ=2/K$ in the gaussian approximation. In our case however, N = 500 and P < 4 such that neither effect significantly impacts the comparison in the cold plasma case.

Note as well that due to our use of the fluctuation distribution for σ, we do not expect $χν2$ to be exactly unity for an ideal model. The values and distributions are comparable relative to each other, but absolute values of $χν2<1$ do not necessarily imply overfitting.

The Langmuir data in Fig. 8 were obtained with Minerva inferring values using the DWSE model. This framework and the automation of the analysis using it will be described in Sec. VI.

As is visible in Fig. 5, fluctuations of V_f (and, presumably, of the other plasma quantities) frequently exist over a broad range of frequencies extending from below the sweep frequency of 500 Hz to and beyond several kHz. The plasma parameters can therefore vary significantly between adjacent sweeps as well as during a single sweep (see also Fig. 6). If the aim is to measure plasma density, temperature, and electric potential at the location of each probe, the best solution would obviously be to sweep the characteristic faster than fluctuations with a significant amplitude, which would require sweep frequencies in the range of several 100 kHz,^26,58 or the use of equivalent techniques (such as analysis of harmonics,⁵⁹ use of a “time-domain triple probe,”⁶⁰ or use of a “mirror Langmuir probe”⁶¹), which imply similarly fast variation of the probe bias. The application of these methods has already been investigated or is planned, but up to now, they have not been used regularly for the W7-X divertor Langmuir probes. This is fundamentally due to restrictions of torus hall space and accordingly long cable lengths.

To obtain the evolution of the plasma quantities on timescales significantly below the sweep frequencies, several approaches are possible:

1. Analysis of every half sweep, averaging of the obtained model parameters p over N half sweeps, denoted $⟨p(1/2)⟩$_N.

2. Simultaneous analysis of N adjacent half sweeps, assuming constant model parameters for all these sweeps, denoted p(N/2).

3. Averaging over the raw data of N half sweeps—this can be executed by defining voltage intervals of the sweep voltage and averaging over all current values for each voltage interval; the probe model is then applied to the resulting characteristic. Denoted p(⁠$⟨N/2⟩$⁠).

We compare the effect of these three averaging methods with N = 10 for one fixed probe model (DWSE) in Fig. 10. The comparison is performed for three time intervals of 0.2 s within W7-X discharge 20180814.007. The three time intervals differ significantly in plasma density and temperature in front of the divertor and in the level of fluctuations. For the first averaging method, we exclude results of half sweeps with the uncertainty of T_e above 50 eV or T_e > 250 eV, indicating too strong fluctuations during the half sweep. We use LS fitting for all three methods. We note from this comparison that the resulting plasma parameters on the timescale of 5 times the sweep period do not differ significantly between the three averaging methods. Small systematic differences are visible for some plasma parameters, however not in the same direction under all discharge conditions. In the following, we prefer to analyze the data with the full half sweep time resolution of 1 kHz, bearing in mind that the fluctuations of the resulting plasma parameters at this time resolution may be misleading but give some indication of the level of fluctuations, whereas the evolution on the 100 Hz scale is reliable.

Effects not included in the forward models can introduce systematic errors in the derived quantities. We will outline these and discuss their importance as well as our measures to limit their impact.

The slope of the ion current is still slightly greater than that predicted even by the DWSE model. By scaling the sheath thickness with an additional free parameter, it was found that this results in a systematic overestimation of T_e and n_e within the error bars of the Bayesian analysis. This matches Weinlich and Carlson’s results³⁷ who obtain their best results with a correction factor of α = 1.2. We have identified a possible source of this discrepancy and aim to correct it in the next version of this analysis: The sheath thicknesses calculated in Sec. III were obtained for each “element” of the probe separately—actively driven probe, unbiased wall, and counter electrode area. This can only be a first order approximation. Realistically, the transition from the wall sheath to the probe sheath would be smooth and not a sudden step as we assume in order to derive a sheath thickness difference for components with different V_bias, as shown in Fig. 4. This would reduce the shadowing of the probe by the wall sheath and thus increase the predicted non-saturation of the ion current. Additionally for the double probe models, it might be necessary to consider the modification of the sheath potential in those areas from which the current is drawn. Since this counter electrode area is only typically larger than that of the driven probe by the factor $β=O(5)$⁠, the current density at the counter electrode can become non-negligible. To fully account for this, a better understanding of the virtual double probe mechanism or direct numerical simulation is likely required, possibly by extending the work by Bergmann et al.⁴⁰ 

When the data are separated into half-periods for analysis, attention must be paid to systematic differences between profiles recorded with the increasing or decreasing bias voltage. Plotting the inferred parameters of the up- and downsweeps separately, we observe small systematic deviations of the parameters within the error bars. At high sweep frequencies f_sweep ≳ 500 kHz, such hysteresis effects are expected in multiple regions of the characteristic due to different mechanisms.⁶² Using a sweep rate of only 500 Hz, we would not expect to see such effects in our system. The reason is potentially the delay of the current signal relative to the voltage signal when measured on the power supplies due to the propagation time of 0.4 µs. This is shorter than the 2 µs sampling interval of the ADCs, so we have not corrected for it yet. As we cannot exclude the presence of some physical effect, we choose not to fit data from adjacent sweeps together. This would only obfuscate this phenomenon without giving insight into its origins.

Bayesian analysis is a more general approach to our problem that is not limited to simple constraints on parameters and can give us meaningful parameter uncertainties. The equivalent confidence intervals of the LS method are not strictly valid for non-linear models and cannot capture the full information Bayesian or maximum likelihood inference provides. For further physics analysis, we require correctly determined uncertainties; therefore, we will use Bayesian analysis to routinely evaluate our measurements.

Unlike for the LS evaluation where the bounds on the parameter values were only set for practical reasons to help the fitting algorithm, prior distributions are mathematically necessary in the Bayesian approach. As its name suggests, it encodes prior knowledge of the parameter values. We use the same values as in Table II for uniform priors, except β that we approximate with a truncated gaussian distribution (⁠$β=N(σ=4,s=4)β∈[1,50]$⁠). This restricts it more strongly to the physically plausible and theoretically expected range.³⁸ 

The Bayesian analysis is performed in a framework we will introduce in Sec. VI that uses acyclic directed graphs to show the conditional dependencies of each model. The Langmuir model graph is shown in Fig. 11. The prediction node contains the implementation of the formulas of the models to calculate the predicted current. The only external information required is that of the magnetic field calculated by the variational moments equilibrium code (VMEC).⁶³Bayes’s formula,⁶⁴ 

is used to infer the posterior distribution of the parameters F given the data D.

The Hooke and Jeeves algorithm⁶⁵ is used to find the approximate maximum of this posterior distribution, i.e., the maximum a posteriori (MAP). Because our parameter space has only few dimensions, this process is very quick, typically requiring $O(10)$ steps to find the MAP. The full distribution is then reconstructed by Markov chain Monte Carlo (MCMC) sampling. We initialize MCMC at the MAP and let it run for $O(104)$ iterations to explore the posterior and generate samples. The result of one such optimization can be seen in Fig. 12. By integrating over parameters (marginalization), we may obtain projections of the posterior distribution showing parameter value probabilities and covariances. Equivalently, we count the samples on a grid of two parameters, thereby creating two-dimensional histograms of the posterior probability density for different parameter combinations, such as the one shown in Fig. 13. These reveal a strong covariance of the parameters and asymmetric deviations from the bivariate gaussian distribution. Marginalizing an additional parameter and applying a kernel density estimate (KDE)^66,67 approach result in a one-dimensional (1D) probability density function (PDF) for each parameter, shown in Fig. 14. The width assigned to each sample in the KDE is determined by Scott’s window.⁶⁷ For our model, the posterior is typically unimodal, with the MAP quite close to the weighted mean of the distribution.

The use of Bayesian analysis for Langmuir probe evaluation has so far only been reported for cold plasma applications outside of fusion physics.^68,69

The inferred plasma parameters of all diagnostics on W7-X should be available to every member of the team for higher-level analysis. It is therefore crucial to automate and integrate the analysis within the project’s infrastructure. We achieve this using the Minerva software framework for Bayesian analysis developed by Svensson and Werner,⁶ which among other benefits provides this capability. Its idea is the following: In the field of fusion physics especially, there are many diagnostics with different operation principles measuring the same quantities. Since the plasma must have one true value of, for instance, density at each point and time, it makes sense to use evidence from multiple sources constrained by a common prior and find the most likely result. The directed graph shown in Fig. 11 is crucial to this approach: It can be expanded to include multiple models making observations of a shared set of priors. Implementing our models in this framework thus allows for a later inclusion of the LPs in a joint analysis with other diagnostics. It furthermore opens the possibility of using the developed forward models at Joint European Torus (JET) or Mega-Ampere Spherical Tokamak (MAST) where Minerva is also used,^70,71 standardizing analysis and enabling objective comparison of results.

To store the inferred values and uncertainties, we use the center of mass and standard deviation of the 1D KDE of each parameter. This is effectively the same as making the approximation of an independent gaussian distribution of each parameter. This last step may seem to negate the benefit of using a more general Bayesian approach but makes sense when considering two different use cases: For analysis based on the Langmuir data alone, the symmetric and independent confidence regions are precise enough and widely understood. For more careful analysis, for example, joint with other diagnostics, the posterior would be recalculated for the entire graph. This is in our case much quicker than storing and loading a representative number of samples or a description of the posterior with additional statistical moments. The procedure is included in a scheduler system such that analyses can be run automatically after each discharge.

In this paper, we described the Langmuir diagnostic system of the probes used in the first divertor campaign of W7-X. A special facetted design of the tips was chosen to mitigate the problems of analyzing flush mounted probes while ensuring the probes withstand the high heat flux. Despite the facetted probe shape, the shallow incidence angle of the magnetic field required going beyond the textbook Langmuir model for the analysis of the characteristics. We compared six models, generalizing and then adapting them to our special use case in a workflow that can easily be expanded to include variations of assumptions or additional models. The I–V characteristics recorded are best matched by a model including the effect of a finite counter electrode, important to describe the electron region at lower T_e, and sheath growth over the ion current region, accounting for the magnetic field. We could show that significantly different plasma parameter values (e.g., up to a factor two lower temperature) resulted when using this model, which best described the characteristics in low temperature cases. Notably, this also holds true at higher temperature plasma conditions, where the WSE/DWSE models indicate significantly reduced T_e values, although the fit quality of the different models is similar. This information is crucial for the LPs to achieve their goal in aiding in the understanding of detachment in the W7-X island divertor. Finally, we describe the automation of our evaluation to make plasma parameters available to every user of the W7-X data archive. Our evaluation is the first of LPs in fusion plasmas using Bayesian probability. This permits the rigorous treatment of uncertainties, revealing significant correlation between inferred plasma parameters. This is done in the analysis framework Minerva, opening the path to a joint analysis of multiple diagnostics.

Raw data were generated at the Wendelstein 7-X large scale facility. Derived data supporting the findings of this study are available from the corresponding author upon reasonable request.

The authors would like to thank S. Freundt, S. Klose, D. Rondeshagen, and J. Kügler for their work and support in building and maintaining the Langmuir system and D. Brida and T. Sunn Pedersen for their advice and valuable discussions.

This work was carried out within the framework of the EUROfusion Consortium and has received funding from Euratom research and training programs 2014–2018 and 2019–2020 under Grant No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

The current flowing to the Langmuir probes is determined by measuring the voltage drop across shunt resistors. The shunt resistance must be small to limit the voltage drop to within the output voltage range of the differential amplifier, in this case ±10 V. This implies measuring a small difference on a large baseline, up to −200 V in the ion saturation mode. The amplifier (Analog Devices AD629) was therefore chosen for its high common mode rejection ratio of 95 dB at the sweep frequency of 500 Hz. In practice, this ratio is reduced by the periphery of the DAmp. We correct for the common mode by determining from the pre-plasma phase a time offset and factor that must be applied to V_bias to match the observed current. This correction current is then calculated for the entire discharge duration and subtracted. This procedure corrects for both the common mode and capacitance at the same time since both are a linear function of the voltage. The capacitive current $Icap=CdVdt$ is dominated by the cable capacitances that were measured individually and are typically (20.0 ± 0.5) nF.

After the subtraction, the pre-plasma current signal only fluctuates about 0 A with an amplitude of ≈5 mA. This remaining background is still periodic with 500 Hz but no longer a linear function of voltage.

The voltage drop proportional to current is not measured directly at the terminals of the shunt resistor. Instead, one DAmp contact connects to ground via a common track for all channels on the measurement card, as shown in Fig. 15. Unfortunately, this means that each DAmp also sees a voltage drop proportional to the other probe currents and resistance of relevant track segments. This effect was considered negligible before and indeed is if the shunt resistance is large compared to that of the track. For the operation in WEGA, and OP1.1 in W7-X, this was mostly true as the currents were low and the 10 Ω shunt resistor could be used. With the higher densities and larger probes of OP1.2, higher currents were recorded, which necessitated the use of smaller shunt resistor values. This made the “common track resistance” error comparable to other effects. We determined the value of the track resistances through direct measurement and careful calibration of the measurement cards with known currents. From these, we derive a formula to calculate the recorded voltages given the currents,

Here, $I8$ is an 8 × 8 identity matrix, I_i and V_Ai are, respectively, the current and voltage on the channel i, R_Si is the shunt resistance used on that channel, and R_ΣTi is the total resistance from that channel to ground (⁠$RΣn=∑i=ni=8RTi$⁠). Additionally, R_Si must be corrected by the resistance of the jumper R_J selecting the shunt value that contributes another ≈(110 ± 30) mΩ. By inverting this matrix, we can calculate the currents from the recorded voltages. The values of the resistances and errors are shown in Table V.⁷² To minimize the mutual influence of the currents upon each other, we typically used channels 1–3 for voltage measurements and selected channels for attributed current measurements considering the expected signal strength, such that the probe with the lowest expected current was connected to channel 4 and the highest signal was recorded on channel 8.