Introduction

In the context of an increasing world energy production based on limited fossil fuel reserves and with a doubling time of 30–40 years in the last century, e.g. from 6000 million tons of oil equivalent (Mtoe) in 1973 to more than 13,000 Mtoe in 2013 [1], controlled thermonuclear fusion could be the ideal candidate to fulfil the future energy demand. Although magnetic confinement fusion based on tokamak is currently the most advanced option, many technological and physics challenges are still on the way to be solved before reaching the technological state of a commercial fusion reactor. In particular, dust and impurities originating from wall erosion can pollute the plasma core and decrease its stability and energy confinement by radiation [2, 3]. Besides, modern tokamaks like the International Thermonuclear Experimental Reactor (ITER) will not be able to use traditional carbon as a plasma-facing material due to its high tritium retention. Tungsten was retained instead for its low tritium retention associated with good mechanical properties, its resilience to intense heat fluxes and its low erosion rate [4]. Unfortunately, non-fully ionized impurities such as tungsten can represent a major threat for the tokamak performances because of intense radiation—especially due to line emission in the soft X-ray (SXR) range 0.1–20 keV [2]. Low concentrations of tungsten impurities in the order of 10–4 in the plasma core are sufficient to compromise the plasma performance of ITER and its goal of \(Q = P_{fusion} /P_{heating} \ge 10\) [5]. This implies that the development of X-ray diagnostic tools becomes crucial to monitor the local density of impurities in the plasma, to study their transport and to identify actuators that can prevent or mitigate their central accumulation [6, 7]. However, the correct estimation of the local impurity concentration relies on the proper characterization of the detector spectral response as well as the detailed emissivity features of the plasma, namely the filtered impurity cooling factor [5]. In this context, a synthetic diagnostic becomes a very valuable tool to study the tomographic reconstruction capabilities, to validate diagnostic design as well as to assess the error propagation during the reconstruction process [8]. Therefore, a complete synthetic diagnostic tool should consider the following steps:

  • Detector signal and/or spectrum reconstruction from a given geometry and plasma scenario,

  • Tomographic inversion to estimate the local emissivity from the synthetic measurements,

  • Impurity density reconstruction and subsequent transport analysis.

The goal of this contribution is to give some highlights on recent studies related to each of these three steps, for the development of SXR synthetic diagnostic tools in tokamak plasmas. In the second section, the potential of artificial neural network for fast tomographic inversions is briefly discussed. In the third section, a synthetic diagnostic is employed in order to predict the observed SXR spectrum in different tokamak scenarios. In the fourth section, the influence of perturbative noise on the reconstruction of impurity transport coefficients is investigated. Finally, some conclusions and perspectives are given in the last section.

Soft X-ray Tomography and the Potential of Neural Networks for Fast Inversions

As introduced in the previous section, reliable SXR diagnostic tools are essential to monitor the local impurity density, study W transport in the plasma core and identify actuators to avoid W central accumulation. Furthermore, the impurity density can exhibit poloidal asymmetries, for instance in the case of plasma rotation (induced by a neutral beam of fast particles) [9] or by the local modification of the electrostatic potential (e.g. during ion cyclotron resonance heating) [10]. It has been shown that these asymmetries can significantly impact the radial impurity transport [7]. Therefore, 2D tomography diagnostics become crucial to estimate the tungsten concentration profile in the plasma, quantify the 2D poloidal distribution and identify relevant impurity mitigation strategies [11].

Hereafter, the geometry of the SXR diagnostic of Tore Supra [12], consisting of two detector-pinhole cameras as depicted in Fig. 1, is considered. In the Line-of-Sight (LoS) approximation, the measurements \(f_{i}\) of the i-th detector (in W m−2) looking at the plasma is given by the line integral:

$$ f_{i} = \varphi_{i} /E_{i} = \mathop \int \limits_{LoS}^{ } \varepsilon^{\eta } \left( {x,y} \right)dr_{i} + \tilde{f}_{i} $$
(1)
Fig. 1
figure 1

Layout of the forward and inverse problem related to tokamak plasma tomography. The inverse problem is solved using an artificial neural network. The LoS geometry of the two SXR cameras of Tore Supra is plotted in white on the tomogram

Full size image

where \(\varphi_{i}\) is the impinging power on the detector (in W), \(E_{i}\) denotes the geometrical etendue (in m−2) of the detector-aperture system, \(\varepsilon^{\eta } \left( {x,y} \right)\) is the 2D emissivity field (in W m−3) filtered by the spectral response \(\eta \left( {hv} \right)\) of the detector and \(\tilde{f}_{i}\) the perturbative noise in the i-th channel. For tomographic purposes, the forward problem can be discretized as follows:

$$ f_{i} = \mathop \sum \limits_{j} T_{ij} { }\varepsilon_{j} + \tilde{f}_{i} $$
(2)

where \(\varepsilon_{j}\) is the plasma emissivity in the j-th pixel and the transfer matrix elements \(T_{ij}\) correspond to the response function of the detection system. Unfortunately, reconstructing the local SXR emissivity consists in an ill-posed inverse problem by nature, with in addition a sparse measurements set \(\left\{ {{\text{f}}_{{\text{i}}} } \right\}_{{{\text{i}} = 1..{\text{N}}}}\) due to the lack of available space in tokamaks. An adequate regularization procedure is thus required to perform the reconstruction.

The most common method relies on the Tikhonov regularization [13]. In this case, some a priori information is added in the reconstruction process through a regularization matrix \(H\) that imposes smoothness on the emissivity profiles, for instance by minimizing the Fisher information. It is then possible to extract a meaningful solution by finding the minimum of a functional:

$$ {\varvec{\varepsilon}}_{0} = \arg \mathop {\min }\limits_{{\varvec{\varepsilon}}} \left( {\left| {{\varvec{f}} - {\varvec{T}}.{\varvec{\varepsilon}}} \right|^{2} + \lambda {}_{\user2{ }}^{{\varvec{t}}} {\varvec{\varepsilon}}.{\varvec{H}}.{\varvec{\varepsilon}}} \right) = \left( {{}_{\user2{ }}^{{\varvec{t}}} {\varvec{T}}.{\varvec{T}} + \lambda {\varvec{H}}} \right)^{ - 1} .{}_{\user2{ }}^{{\varvec{t}}} {\varvec{T}}.{\varvec{f}} $$
(3)

where the regularization parameter \(\lambda\) acts as a balance between \(\chi^{2} = \left| {{\varvec{f}} - {\varvec{T}}.{\varvec{\varepsilon}}} \right|^{2}\) minimization and smoothness of the solution. The value of \(\lambda\) can be optimized based on several different methods [14]. Inversion methods alternative to the Tikhonov regularization can be found in the literature, e.g. based on Bayesian inference [15], a Monte-Carlo approach [16] or on a genetic algorithm [17]. However, most of these methods still require a significant computational time although significant effort is being performed towards real-time applications [18]. In this context, deep learning and artificial neural networks [19] represent the ideal candidate in the perspective of real-time impurity control. Neural networks are able to fit any multidimensional function in a given input range—with the adequate network complexity and training process—by adjusting properly the weights and biases of the neurons [20]. Besides tomography [19] and image recognition [21], they have been used in support of simulation with high computational cost, e.g. turbulent transport modelling [22], and to make empirical correlations from large experimental databases, e.g. disruption prediction [23].

A neural network with two fully connected layers of 30 sigmoid neurons has been recently developed in the Tore Supra geometry to achieve fast SXR inversions, as depicted in Fig. 1, and trained with a large database of synthetic profiles [24]. The input layer is composed of 69 neurons, corresponding to the number of SXR channels used for the inversion and of an output layer of 2500 neurons to produce 50 × 50 pixels tomograms. The training loss that assesses the network performances is defined by the mean square error cost function:

$$ {\text{C}}\left( {{\text{w}},{\text{b}}} \right) = \frac{1}{{2{\text{n}}}}\mathop \sum \limits_{{\text{x}}} \left[ {{\text{y}}\left( {\text{x}} \right) - {\text{a}}\left( {{\text{x}},{\text{w}},{\text{b}}} \right)} \right]^{2} $$
(4)

where \(\left( {{\text{w}},{\text{b}}} \right)\) denotes the weights and biases of each neuron, \(\left( {{\text{x}},{\text{y}}\left( {\text{x}} \right)} \right)\) represent a training sample, i.e. measurements x and tomogram \({\text{y}}\left( {\text{x}} \right)\), with \({\text{a}}\left( {{\text{x}},{\text{w}},{\text{b}}} \right)\) the corresponding output of the network. A stochastic gradient-descent method [25] is used to iteratively updates the weights and biases in the direction of the negative gradient of the cost function, computed using a backpropagation algorithm [26]. The evolution of the training loss over the training process is shown in Fig. 2, where 90% of the data were used for training and 10% were kept outside of the training set for validation purposes.

Fig. 2
figure 2

Minimization of the cost function value of the neural network, as defined in Eq. (4), after # epochs of the training process

Full size image

Here, we test this neural network with synthetic emissivity patterns from the validation set, exhibiting poloidal asymmetries as presented in Fig. 3 and with 2% of Gaussian noise in the measurements to assess the neural network robustness. It can be seen that the network is able to properly recover the synthetic emissivity patterns. However, the emissivity error can locally peak from 15% to up to 40% as shown in Fig. 3g–i), which could be detrimental for the SXR analysis of impurities. The reconstruction error could be further reduced by optimizing the network parameters (number of hidden layers and neurons per layer) and the training process, e.g. using the dropout method. The relatively high level of error observed on Fig. 3 is the drawback of the simplicity of fully connected layers as compared with deeper de-convolutional networks, for which the structure is specifically adapted for image recognition and tomography tasks [19, 27].

Fig. 3
figure 3

Tomographic inversions performed using phantom 2D emissivity profiles with a high-field-side, b hollow, c low-field-side shapes. d–f Associated tomograms reconstructed by the neural network and g–i the corresponding error map between reconstructed and initial tomograms

Full size image

The main advantage of the neural network approach is that, once properly train with a synthetic or experimental database, the neural network can perform inversions in a time < 0.1 ms, two orders of magnitude faster than Tikhonov regularization and compatible with real-time control. Two options can be considered while applying a neural network to experimental data. The first option is to directly train the neural network on a database of experimental tomograms already obtained from another method, typically using Tikhonov regularization, like in [19]. This option allows performing reconstructions much faster with respect to the traditional regularization, at the cost of an additional error induced by the neural network, but requires the prior existence of a reliable and complete experimental database. The second option is to train the neural network on a synthetic database in a first step, and in a second step to apply it on experimental data. The advantage of such an approach is to be able to perform the tomographic task ab initio, for instance from the beginning of tokamak operations. However, the reliability of the latter approach depends largely on the capabilities of the synthetic diagnostics to describe accurately the reality of the experiments.

GEM Synthetic Diagnostic for Impurity Transport Studies on the WEST Tokamak

The Tore Supra tokamak was recently upgraded into WEST—for Tungsten (W) Environment in Steady-State Tokamak, to serve as a test-bed for ITER divertor components in long pulse operation [28]. It was also the opportunity to develop a new SXR detection system based on the Gas Electron Multiplier (GEM) technology [29, 30], after laboratory tests of a GEM prototype at CELIA [31] and on ASDEX-Upgrade [32]. Indeed, the semiconductor diodes traditionally used in tokamaks will not whistand the harsh ITER environment in the D-T operations and more robust technologies based for instance on a gaseous detection volume have to be considered [33]. A specificity of the GEM WEST detectors is that they will operate in photon-counting mode, allowing to classify incident photons in energy bands and thus recovering some spectral information on the SXR radiation [34]. This will open new possibilities to disentangle the different SXR contribution from W emissivity, i.e. line emission, radiative recombination and Bremsstrahlung radiation.

In this context, a synthetic diagnostic has been developed in order to investigate the reconstruction capabilities of the GEM WEST cameras [8]. The diagnostic is composed of two cameras placed in the same poloidal cross-section, as shown in Fig. 4. The synthetic diagnostic includes the SXR source from a given plasma scenario, the photoionization probability, electron cloud transport in the detection volume using the Magboltz code [35], and tomographic reconstruction of the radiation from the GEM signal in the WEST geometry [36]. In this section, we show a set of three simulations of the acquired SXR spectrum by the GEM cameras in the case of a standard WEST Physics basis scenario with heating power of 12 MW, Ip = 0.6 MA including W impurities [37]. The SXR emissivity is modelled as follows:

$$ \varepsilon^{SXR} = n_{e}^{2} \left[ {L_{H} \left( {T_{e} } \right) + c_{W} L_{W} \left( {T_{e} ,D,V} \right)} \right] $$
(5)
Fig. 4
figure 4

Poloidal cross-section of the WEST tokamak, with the geometry of the lines-of-sight of the new GEM diagnostic in red. The relative signal intensity on each pixel of the horizontal and vertical GEM cameras are represented for a WEST emissivity scenario in pure hydrogen plasma

Full size image

with \(n_{e}^{ }\), \(T_{e}\) the electron density and temperature and where \(L_{H}\) and \(L_{W}\) denote the Hydrogen and W cooling factors, respectively. The W spectral properties, as well as recombination and ionization coefficients of W ions are obtained from the open ADAS database [38]. The scenario includes the possibility to choose the diffusion and convective W transport coefficients (D,V) to modify the W ionization equilibrium.

The Fig. 5a reports the synthetic GEM spectrum obtained for the three scenarios considered: (1) pure plasma hydrogen plasma, i.e. a tungsten concentration \(c_{W} = 0\), (2) \(c_{W} = 10^{ - 4}\) with an intermediate transport case (D,V)x1 and (3) \(c_{W} = 10^{ - 4}\) with a strong transport case (D,V)x10, i.e. where the transport coefficients have been increased by one order of magnitude to highlight the transport effects. The transport coefficients values in the (D,V)x1 scenario correspond to D = 0.5 m2 s−1, V = − 1 m s−1 in the plasma core \(r/a < 0.6\) and D = 2 m2 s−1, V = − 5 m s−1 at the plasma edge. A significant increase of the SXR radiation is observed around 2 keV. This is expected since W line emission is dominant in this spectral region for the considered plasma core temperature Te ≈ 3–4 keV, while the Bremsstrahlung radiation becomes dominant for photon energies above 4 keV. The impact of transport intensity on the W SXR cooling factor LW can be seen on Fig. 5b. In particular, the enhanced transport reduces the W radiation in the plasma core, since high W ionization states are more easily transported from the core to the edge, decreasing the mean W effective charge in the core.

Fig. 5
figure 5

a GEM synthetic spectrum acquired in a WEST emissivity scenario (line-integrated spectrum summed over all GEM channels) with different assumptions on tungsten impurity concentration. b W SXR GEM filtered cooling factor in the two considered transport scenarios

Full size image

Reconstruction of Impurity Transport Coefficients and Impact of Perturbative Noise

In this section, we investigate the possibility of reconstructing the impurity transport coefficients (D,V) from given measurements of the local impurity density \(n_{S}\) in the presence of different levels of perturbative noise. The impurity density is generally estimated by reformulating Eq. (5):

$$ n_{S} \approx \frac{{\varepsilon^{SXR} - n_{e}^{2} L_{H} \left( {T_{e} } \right)}}{{n_{e}^{ } L_{S} \left( {T_{e} ,D,V} \right)}} $$
(6)

given that the species S (e.g. W) is the dominant impurity contributing to the emissivity. This can be performed experimentally by SXR tomography as seen in the previous sections, coupled with other diagnostics like Visible–UV spectroscopy to cover a wider radial range. However, in this section, the details of the detector response, impurity spectral features and tomographic reconstruction are not included in the analysis and the error in the impurity density reconstruction is rather described by a statistical noise with a Gaussian distribution and without systematic errors, for simplicity purposes. The issue of systematic errors in the geometry or camera calibration is thus let for further investigations. The (D,V) reconstruction relies on the Gradient-Flux (G-F) method [39, 40], which is applicable when the impurity density profile is evolving in time while the plasma conditions remain stable. Indeed, in this case the radial flux of impurities can be expressed as:

$$ \vec{\Gamma }_{S} \left( {r,t} \right) = - D\left( r \right)\vec{\nabla }_{r} n_{S} \left( {r,t} \right) + n_{S} \left( {r,t} \right)\vec{V}\left( r \right) $$
(7)

where the radial flux and gradient of impurity are time dependent while the diffusive and convective coefficients remain constant. The conditions can be fulfilled during the recovery phase of a sawtooth (ST) crash or just after a controlled impurity injection at a trace level. The method allows extracting experimentally the transport coefficient thanks to a linear fit:

$$ \frac{{\Gamma_{S} }}{{n_{S} }}\left( {r,t} \right) = - D\left( r \right)\frac{{\nabla_{r} n_{S} }}{{n_{S} }}\left( {r,t} \right) + V\left( r \right) $$
(8)

where \(- D\left( r \right)\) is the slope and \(V\left( r \right)\) is the offset. The impurity flux can be reconstructed from:

$$ \Gamma_{S} \left( {r,t} \right) = - \frac{1}{r}\mathop \int \limits_{0}^{r} \frac{{\partial n_{S} }}{\partial t}\left( {r^{\prime},t} \right)r^{\prime}dr^{\prime} $$
(9)

The G-F linear fit in Eq. (8) is subject to experimental uncertainties, in particular in the plasma region where the impurity density gradient is low (i.e. flat radial profile). It is thus valuable to develop a synthetic diagnostic tool to assess the reconstruction capabilities for different levels of noise in the reconstructed impurity density.

Hereafter, we will consider an ideal tokamak in circular geometry with a major radius of R = 2.5 m and a minor radius of a = 0.5 m. The assumed (D,V) coefficients are shown in Fig. 6a, together with the steady-state impurity density profile obtained using Eq. (8) when \(\Gamma_{S} \left( {r,t} \right) = 0\):

$$ \frac{{\nabla_{r} n_{S} \left( r \right)}}{{n_{S} \left( r \right)}} = \frac{V\left( r \right)}{{D\left( r \right)}} $$
(10)
Fig. 6
figure 6

a W transport coefficient profiles considered in the simulation. b W density profile resulting from the chosen transport coefficients

Full size image

Equation (10) leads to the following impurity density profile:

$$ n_{S} \left( r \right) = n_{S,0} \exp \left( {\mathop \int \limits_{0}^{r} \frac{{V\left( {r^{\prime}} \right)}}{{D\left( {r^{\prime}} \right)}}dr^{\prime}} \right) $$
(11)

The simulation is simply performed by solving numerically the continuity equation for the impurities:

$$ \frac{{\partial n_{S} }}{\partial t} + \overrightarrow {{\nabla_{r} }} .\vec{\Gamma }_{S} = S_{ext} $$
(12)

The obtained impurity 2D map is presented in Fig. 7a. The time evolution of the impurity density for a few radial locations is also displayed in Fig. 7b. ST crashes are simply implemented by transiently increasing the diffusive coefficient up to D = 30 m2 s−1 in the region \(r/a < 0.3\), as presented in Fig. 8a. An impurity injection is also added by putting a non-zero source term \(S_{ext} > 0\) at t = 0.11 s at the plasma edge, i.e. for \(r/a > 0.95\).

Fig. 7
figure 7

Simulated time evolution of the impurity density for different radial locations. a Impurity density 2D map (logscale), where the time window used for the reconstruction is indicated by red dash lines. b Times traces for different radial locations, where the impurity injection is indicated in red (LBO—for Laser-blow off trigger)

Full size image
Fig. 8
figure 8

Reconstructed transport coefficients using the ST recovery phase, i.e. excluding the effect of the controlled impurity injection. a Diffusion coefficient profile during a ST crash. b–d Reconstruction of the transport coefficients with a level of Gaussian noise of b 0.1%, c 0.5% and d 2%, where \(D_{ini}\), \(V_{ini}\) are the transport coefficients used in the simulation, \(D_{rec}\), \(V_{rec}\) are the reconstructed ones with error bars corresponding to \(\pm 2\sigma\) in the G-F linear fit—not plotted for large errorbars, i.e. \(r/a > 0.3\) in b and c due to visibility issues. The black dashed line “ST” denotes the sawtooth inversion radius at \(r/a = 0.3\)

Full size image

An attempt of reconstructing the transport coefficients in the ST recovery phase, excluding the impurity injection phase i.e. for the time window \(t = \left[ {0.033s;0.099s} \right]\), is shown in Fig. 8 for 0.1%, 0.5% and 2% of perturbative noise. The impurity density is assumed to be reconstructed over the full radial range including a statistical error with a Gaussian distribution varying from 0.1 to 10%. The impurity flux is estimated using the obtained noisy density in Eq. (10) and the reconstructed coefficients (\(D_{rec}\), \(V_{rec}\)) are estimated from the G-F linear fit described by Eq. (9). The error bars in the linear fit correspond to \(\pm 2{\upsigma }\), i.e. 95% of values around the mean for a normal distribution. As a result, it is observed that the reconstruction is only reliable inside the ST inversion radius \(r/a < 0.3\). This is expected since it is the region of significant time evolution of the impurity flux during the ST recovery phase. While the reconstruction is still robust for 0.5% of perturbative noise and \(r/a < 0.3\), see Fig. 8c, a Gaussian noise of \( \gtrsim 2\%\) is enough to significantly perturb the reconstruction over the full radial range, as shown on Fig. 8d.

Finally, the GF reconstruction method is applied just after the impurity injection, where the time window used for the reconstruction is indicated by the red dash lines on Fig. 7a. The reconstruction results are presented in Fig. 9a–c) for different levels of Gaussian noise—0.1%, 2% and 10%—together with the linear fit of Eq. (8) in Fig. 9d–f. Overall, it can be seen that the transport coefficient profiles are recovered over the full radial range with a better robustness against statistical noise than during the ST recovery phase alone. Indeed, a strong noise of 10% is necessary to significantly perturb the (\(D_{rec}\), \(V_{rec}\)) profiles. However, it is observed that the reconstruction is less reliable in the limit regions \(r/a < 0.1\) and \(r/a > 0.8\) where the coefficients can be underestimated (in absolute value), even for 0.1% and 2% of Gaussian noise, as shown in Fig. 9a, b. This can be interpreted by the fact that both the radial gradient and flux of impurities are lower in these regions, decreasing the quality of the G-F linear fit. As a conclusion, the reconstructions during the ST recovery phase and during the impurity injection phase showed that the G-F method is mostly robust in the regions with a significant radial impurity gradient and for strong time evolution of the impurity profile.

Fig. 9
figure 9

Reconstructed transport coefficients using the time evolution of the tungsten density during the impurity injection phase, assuming different level of Gaussian noise in the impurity density reconstruction: a 0.1%, b 2% and c 10%, where \(D_{ini}\), \(V_{ini}\) denote the coefficients used in the simulation and \(D_{rec}\), \(V_{rec}\) the reconstructed ones with error bars corresponding to \(\pm 2\sigma\) in the G-F linear fit. d–f Associated G-F linear fit of the normalized impurity fluxes and radial gradients, based on Eq. (8)

Full size image

Conclusion

In this paper, some recent studies related to the development of X-ray synthetic diagnostic tools have been presented, in the prospect of impurity transport monitoring. First, the potential of neural networks to perform tomographic inversions with real-time capabilities was demonstrated based on a synthetic database in the geometry of the Tore Supra SXR diagnostic. A good capability of the network to learn and reproduced the emissivity patterns was found with 2% of Gaussian noise in the measurements, although the local reconstruction error can peak from 15% to up to 40%. Secondly, a GEM synthetic diagnostic was applied in a WEST scenario to illustrate the influence of plasma parameters, in particular the tungsten concentration level and transport effects, on the acquired spectrum and the possibility to disentangle the different contributions to the SXR spectrum. Finally, the robustness of the reconstruction of the impurity transport coefficients (D,V) against 0.1–10% of perturbative Gaussian noise was investigated with a using the G-F method, during a ST recovery phase and during an impurity injection phase. The study allowed to reveal the plasma regions in which the (D,V) reconstruction is the most sensitive to experimental uncertainties. More specifically, the (D,V) reconstruction is no more valid outside of the ST inversion radius, here for \(r/a > 0.3\), during the ST recovery phase. In the impurity injection phase, the reconstruction is valid over the full profile up to 10% of noise, although less reliable in the limit regions \(r/a < 0.1\) and \(r/a < 0.8\) of lower gradients and flux of impurities for the studied scenario.