Implementing an X-ray tomography method for fusion devices

In fusion devices, the X-ray plasma emissivity contains essential information on the magnetohydrodynamic activity, the magnetic equilibrium and on the transport of impurities, in particular for tokamaks in the soft X-ray (SXR) energy range of 0.1–20 keV. In this context, tomography diagnostics are a key method to estimate the local plasma emissivity from a given set of line-integrated measurements. Unfortunately, the reconstruction problem is mathematically ill-posed, due to very sparse and noisy measurements, requiring an adequate regularization procedure. The goal of this paper is to introduce, with a didactic approach, some methodology and tools to develop an X-ray tomography algorithm. Based on a simple 1D tomography problem, the Tikhonov regularization is described in detail with a study of the optimal reconstruction parameters, such as the choice of the emissivity spatial resolution and the regularization parameter. A methodology is proposed to perform an in situ sensitivity and position cross-calibration of the detectors with an iterative procedure, by using the information redundancy and data variability in a given set of reconstructed profiles. Finally, the basic steps to build a synthetic tomography diagnostics in a more realistic tokamak environment are introduced, together with some tools to assess the capabilities of the 2D tomography algorithm.


Introduction
In fusion devices, the X-ray plasma emissivity contains essential information on the magnetohydrodynamic activity, the magnetic equilibrium and on the transport of impurities, in particular for tokamaks in the soft X-ray (SXR) energy range of 0.1-20 keV. In this context, tomography diagnostics are a key method to estimate the local plasma emissivity from a given set of line-integrated measurements. Unfortunately, the reconstruction problem is mathematically ill-posed, due to very sparse and noisy measurements, requiring an adequate regularization procedure. In practice, adapted reconstruction methods have to be chosen in order to perform the inversion task and select meaningful solutions. Though a great variety of different inversion techniques have been developed for plasma tomography in different Focus Point on Modern Developments in Applications for Plasma Physics. a e-mail: axel.jardin@ifj.edu.pl (corresponding author) energy ranges such as VUV, SXR, hard X-ray or gamma-rays [1][2][3][4][5][6][7][8][9], this paper is intended to focus on one of the most common methods used for SXR tomography in tokamaks, based on the Tikhonov approach [10].
The crucial issue of heavy impurity radiation has been raised by modern tokamaks like ITER-International Thermonuclear Experimental Reactor, which selected tungsten (W) instead of traditional carbon as the main plasma facing material for its divertor, to minimize tritium retention in the walls. In future reactors, W core concentration should be kept below 0.01% to prevent an unmanageable cooling of the plasma by radiation [11]. Besides, it has been shown that the poloidal distribution of impurities can affect their central accumulation [12], highlighting the fact that 2D tomographic tools are crucial to estimate the W radial profile in the plasma, quantify its poloidal distribution and identify relevant impurity mitigation strategies.
This paper is intended for any scientist interested in the field of plasma tomography or any student wishing to develop his/her own tomography algorithm. The goal is to introduce some methodology and tools to develop a tomography algorithm for fusion devices. In the second section, the tomography problem is introduced with a focus on the Tikhonov regularization. In the third section, we study the optimal values of the reconstruction parameters, such as the emissivity spatial resolution and the regularization parameter, based on a simple 1D tomography problem. A methodology is also proposed to perform an in situ sensitivity crosscalibration and position correction of the detectors with an iterative procedure, by using the information redundancy and data variability in a given set of reconstructed profiles. Finally, in the last section, the basic steps to build a synthetic soft X-ray tomography diagnostics in a more realistic tokamak environment are introduced, together with some tools to assess the capabilities of the 2D tomography algorithm.

Generalities and main assumptions
The incident power ϕ (in W ) measured by a detector-pinhole camera looking at the plasma can be expressed in general as: ϕ X Y Z K (x, y, z) ε η (x, y, z)dx dy dz (1) with ε η the plasma emissivity (in W·m −3 ) filtered by the spectral response η(hν) of the camera-where hν is the energy of the incident radiation and K (x, y, z) the geometric function of the detector-pinhole system in the (X, Y, Z ) space. As pictured in Fig. 1, the incident power ϕ onto the detector D is the sum of the contributions from the infinitesimal emissivity volumes dV p in the plasma volume V p : where dV p S e f f D /r 2 is the solid angle under which D is viewed by the volume dV p through the pinhole P, with r the distance between dV p and D. Considering that the spatial Fig. 1 Representation of a detector-pinhole system in a vertical cross section of the plasma variation of the emissivity ε η is small in each section S p of the plasma perpendicular to the pinhole-detector axis (and contained in the volume-of-sight), it is possible to write: where E is called the geometrical etendue (in m 2 ) of the detector-pinhole system. This simplification is referred to as the Line-of-Sight (LoS) approximation. Furthermore, it is possible to show [13] that the geometrical etendue is conserved along the LoS, since the plasma section S p seen by the system increases as r 2 while the solid angle of the detection system as seen by the plasma decreases as 1/r 2 , such that: Provided that the pinhole is small as compared to the detector-pinhole distance L, it is possible to obtain a simplified formula for the geometrical etendue: E S P S D / 4π L 2 .

The inverse problem
Based on Eq. (4), the line-integrated measurements f i of the ith detector, also referred to as brightness (in W·m −2 ), can be given, in the line-of-sight (LoS) approximation, by: where ϕ i denotes the incident power on the ith channel (in W), E i is the geometrical etendue (in m 2 ) of the detector-aperture system, ε η (x, y) is the 2D filtered emissivity field (in W·m −3 ) andf i represents any perturbative noise in the ith measurements (e.g., electronic noise). Discretization of the space containing plasma in N p emissivity pixels leads to the reformulation of Eq. (5) as follows: with ε j the emissivity (hereafter, the upper index η associated with the detector spectral response will be omitted for simplicity) in the jth plasma pixel. The transfer matrix elements  2 Layout of a line-of-sight in a tomographic system, with pixel discretization of the plasma region T i j correspond to the response function of the detection system. In our case, T i j will simply denote the length of the ith LoS in the jth pixel, as depicted in Fig. 2.
It is worth noting that several tokamak plasma tomography methods were first developed using global basis functions (e.g., Bessel and Fourier expansion of the plasma emissivity profile), see for instance the use of the Cormack technique on Alcator C [2], rather than using local basis functions such as square pixels. However, the latter approach gave superior results and has become dominant nowadays in this topic.
For ideal measurements (i.e., without any systematic or statistical error), denoted f i,0 , the previous equation can be rewritten as: In experimental applications, the real measurements f i will result from the experimental deviation δ f i from ideal measurements f i,0 such that f i f i,0 + δ f i . In the same way, the transfer matrix can be expressed as: in the small perturbation approximation, where C i 1 + δC i depicts the ith detector imperfect calibration (with δC i 1) andT i j results from detector and pinhole positioning errors. Therefore, Eqs. (6) and (7) lead to: where the second-order term δC iTi j has been neglected here (it should be kept if δC i 1 is not fulfilled). It is clear from Eq. (8) that the total measurement error δ f i results from: -the detector miscalibration δC i , -the detector-pinhole mispositioningT i j , -the statistical (e.g., electronic) noisef i .
We will assess in the next sections the effect of measurement errors on the tomographic reconstructions, and the possibility and limitations of correcting in situ, thanks to some tomography procedure, the calibration and positioning errors.

Tikhonov regularization
The intuitive approach of obtaining a solution of the tomographic problem would be to find a solution ε rec which minimizes the measurement error χ 2 between the retrofits f rec T ·ε rec of the reconstruction and the actual measurements f meas : Unfortunately, reconstructing the local emissivity is not a straightforward task as it is an ill-posed inverse problem by nature. As we have seen previously, the presence of statistical and systematic errors in the measurements affects the reconstruction. Besides, the measurements set {f i } i 1...N los is usually quite sparse in tokamaks, due to the lack of available space. This implies that the number of plasma pixels N p , necessary to obtain a satisfactory spatial resolution is usually much greater than the number of measurements, i.e., N p N los , making the problem underdetermined. Therefore, a solution with satisfying specific properties should be selected among the infinity of possible solutions, by adding some a priori information in the reconstruction. Many different inversion methods have been developed to solve the inversion problem for tokamak plasmas. Some of them are based on global basis functions like the Fourier-Bessel decomposition [2,3], and more recently on local basis functions (e.g., pixels). The traditional approach is to minimize (or maximize) a functional that includes some a priori knowledge on the plasma emissivity shape, such as in the maximum entropy [4,5], maximum-likelihood [6] and other Bayesian methods [8]. Other approaches involve genetic algorithms [14] or Monte Carlo methods [15]. In the prospect of real-time applications, machine learning and neural networks are also of great interest for performing tomographic reconstructions at low computational cost [9,16].
Nevertheless, in this paper we focus on one of the most common approaches for SXR tomography in tokamak plasmas, with the central emphasis placed on the smoothness of the solution. A. N. Tikhonov introduced such regularization procedure [10], for which the a priori information is added through a regularization matrix H that imposes smoothness on the emissivity profiles. For example, it can aim at minimizing the local spatial gradient H t ∇ · ∇ (the superscript t denotes the matrix transpose operation), with the gradient in the matrix form (e.g., forward difference) along the X direction: with X the spatial step size. A meaningful solution can then be obtained by finding the minimum of a functional φ: where the regularization parameter λ allows to obtain a compromise between the minimization of χ 2 and of the regularization term R t ε · H · ε. The problem aims at finding a solution ε rec that satisfies: The explicit solution ε rec of the minimization problem written in Eq. (12) can be expressed as: where the value of λ is a parameter of the reconstruction and has to be chosen carefully. An empirical approach that can be found in [17] is to choose the ratio of matrix traces in Eq. (13): This method has the advantage of being simple and fast, although it can be far from the optimal value as we will see in the next sections. A second approach, so-called L-curve corner, is to perform a scan of the regularization parameter and to plot the regularization term R as a function of the measurement error χ 2 . Such plot has usually a L-shape with an apparent inflection point or corner which is chosen as the optimal λ value. However, this method is computationally costly since it requires to perform a scan over a wide range of the regularization parameter. Besides, the L-curve corner is not always clearly defined in the case of noisy and sparse measurements, as we will see in the next sections. Another method relies on the discrepancy principle. The goal is to obtain a measurement error χ 2 N equivalent to the expected noise level σ 2 noise (assuming a Gaussian distribution): The purpose is to smooth out emissivity structures that are below the noise level, and keep those of higher intensity-considered as valuable information-in the reconstruction.

Second-order Philips-Tikhonov regularization
For benchmark purposes, two Tikhonov regularization methods are introduced and used in this paper. The first one is a second-order Philips-Tikhonov regularization (PTR) and aims at minimizing the second-order derivative (the Laplacian) of the emissivity. Such a method was for instance developed at JET for neutron tomography [18]. In this case, the regularization operator H PTR is expressed as: Such an operator imposes constraint on the spatial variation of the emissivity gradient and not on the emissivity gradient itself. It thus prevents the appearance of highly fluctuating emissivity structures in the plasma core and sharp profiles at the edge.

Minimum Fisher information
The second reconstruction technique introduced here, so-called Minimum Fisher information (MFI), is one of the most widely applied methods nowadays for SXR tomography in tokamaks, for instance in JET [19], ASDEX Upgrade [20], TCV [1], COMPASS Upgrade [21] or Tore Supra [22] and WEST [23]. For a given 2D emissivity distribution ε(x, y), the associated Fisher information is defined as: Fig. 3 Workflow of the algorithm used to find the solution that minimizes functional of the MFI method, with a convergence criterion of 5%. Symbols in red are quantities susceptible to change after each iteration where the prime denotes the spatial derivative. Numerical transcription of the Fisher information in the regularization operator H MFI gives: aiming at selecting the solution with the least curvature, with W a ponderation matrix defined as: where δ i j denotes the Kronecker delta and ε min > 0 is a lower bound used to avoid singularities for the zero emissivity regions. The ponderation matrix W imposes a constraint of lower gradients in the low emissivity regions, e.g., plasma edge (where W i j value is high). The constraint is relaxed in the plasma core (where W i j value is low), allowing the reconstruction of spatially fluctuating structures in that region. This can be useful, e.g., to monitor MHD activities or impurity asymmetric distributions [24]. It can be noted that, unlike PTR, the prior knowledge of the emissivity field is necessary to estimate the ponderation matrix W. This issue can be overcome with an iterative approach. A first guess, here with a flat emissivity ε (0) 1, is chosen for the first estimation of the ponderation matrix elements W i j using Eq. (19). From this, a new emissivity ε (1) can be calculated using Eq. (13). The procedure can then be repeated k times as depicted in Fig. 3, until convergence ε (k+1) ε (k) is reached with some tolerance of typically 5-10% (usually 3-5 iterations are needed).

Definition of the geometry
We first implement the tomography method for a simple symmetric case, where the plasma consists of several nested homogeneous layers, as represented in Fig. 4. We have thus a 1D problem where the plasma emissivity depends only on the distance from the plasma center ε ε(r ). This symmetry implies that the detector LoS is only defined by its closest distance from the plasma center z LoS , where all LoS are represented horizontally on the figure. Four synthetic emissivity profiles of different shapes are used to assess the capabilities of the tomography method in various plasma conditions: allowing to obtain in the dataset, after normalization to 1, a so-called flat profile (ε 1 ), a peaked profile (ε 2 ), a regular profile (ε 3 ) and a hollow profile (ε 4 ), as shown in Fig. 5.

Onion-peeling method
An intuitive reconstruction method would be to choose N p N los N in order to simply retrieve the emissivity by performing "onion-peeling" of the plasma emissivity. Indeed, since the last plasma layer would only be crossed by one line-of-sight, the last emissivity element could be retrieved using the fact that ε N f N /T N N and afterwards all the other emissivity values, iteratively: Although simple and smart at first sight, this method is giving very poor results in practice. First, it limits drastically the spatial resolution of the reconstructed emissivity distribution. Second, the method shows a high sensitivity to perturbative noise, as depicted in Fig. 6 with a regular profile. Random noise is added here according to Eq. (6), with a zero-mean Gaussian distribution and a standard deviation equal to 2% of the signal intensity. The onion-peeling method implies that the reconstruction errors obtained at the plasma layers δε k+1:N accumulate iteratively to the inner layer δε k , making the reconstruction less and less reliable as we move to the plasma center. Thus, the use of a regularization procedure like Tikhonov regularization is required for plasma tomography. It is worth noting that despite poor reconstruction results, the retrofits match well the line-integrated measurements, see Fig. 6 (right). This is referred to as measurements overfitting.
A tomographic reconstruction using MFI is displayed in Fig. 7, for different values of the regularization parameter λ, in order to illustrate the impact of λ value on the smoothness of the solution. A full λ scan is performed in order to investigate the optimal solution, as shown in Fig. 8. The figure of merit to estimate the quality of the reconstruction is based on the root-mean-square emissivity error R M S em : Eur. Phys. J. Plus (2021) 136:706 It can be seen that there is one optimal value of λ which minimizes RMS em , noted λ opt . It is also visible that the discrepancy principle allows to find a solution λ discr very close to the optimum, while the trace ratio solution λ tr leads here to an overfitting of the measurements. Figure 9 shows the L-curve for two different set of reconstruction parameters. It can be seen on the left that finding the L-curve corner can lead to a solution close to λ opt when the corner is well-defined. Unfortunately, it is not always the case, as shown in Fig. 9 (right), in situations with significant noise and low number of measurements. In such situations, the L-curve corner method is not applicable. For these reasons, the discrepancy principle will be used hereafter for λ optimization. Since the number of plasma layers N p is a free parameter of the reconstruction, it is valuable to study how it affects the overall quality of the reconstruction. Such N p scan is shown in Fig. 10 for the four profiles, together with the associated computation time (normalized to the minimum one). In this specific case, no significant decrease of RMS em is observed for N p > 50, while the computation time is increasing exponentially. In practice, it is therefore recommended to limit wisely the plasma resolution to avoid unmanageable computation cost.

Tomography tests
We will now assess the general capabilities of the PTR and MFI methods for the four different synthetic emissivity profiles: regular, peaked, flat and hollow, as defined in the previous sections. A set of 75 tomographic reconstructions for each profile is performed after adding 1% of random noise with a zero-mean Gaussian distribution. The average reconstruction for each profile is plotted in Fig. 11, together with the envelope made of the lowest and highest reconstructed values along the profile.
As a result, both PTR and MFI show similar capabilities to recover the general shape of each of the profiles. Nevertheless, MFI appears more sensitive to noise (larger envelope) than PTR but exhibits a higher capability to follow steep gradients at the edge of the profiles and in the core (e.g., for the peaked profile). This could explain why MFI is usually preferred for SXR tomography [19][20][21][22][23] and PTR for neutron tomography in the literature [18].

Detectors in situ cross-calibration correction
While statistical noise can only be smoothed out by a proper optimization of the regularization parameter, systematic errors originating from the imperfect knowledge of the detectors sensitivity persist through the reconstructions of different emissivity profiles. Therefore, they can in principle be detected and corrected, provided that the measurements database contains enough statistics.
In this section, we propose a methodology to perform a cross-calibration correction of detectors sensitivity and we test its capabilities and limitations, using the defined set of synthetic profiles in the 1D problem configuration. A calibration error coefficient is applied to each detector with a Gaussian probability distribution where σ 0.1. The method aims at using the redundancy of information between detectors to predict the expected measurements from a given detector. The benchmark between the expected and experimental measurements from this detector can then be used to estimate a cross-calibration correction factor. The iterative procedure is decomposed as follows: 1. Perform the tomographic reconstructions of the full measurements set while omitting the detector #1 out of the reconstruction process. 2. Use the reconstructed emissivity profiles to predict the measurements f rec,1 from the detector #1 using f rec,1 N p j 1 T 1 j ε rec, j . 3. Estimate the mean m f, 1 and variance σ f,1 of the measurements ratio f meas,1 / f rec,1 . It should be expected that m f,1 1 when the detector sensitivity is already well-known. 4. Repeat the steps 1 to 3 for all other detectors in order to estimate m f,1:N los and σ f,1:N los , as shown in Fig. 12. 5. Select the detectors for which the following condition is fulfilled: with K a constant determined empirically (here K 3/2). This condition ensures that the applied correction coefficient C i 1/m f,i is robust against the data statistical variability. 6. Apply the cross-calibration correction factor C k to the detector showing the highest calibration error, i.e., k arg max i m f,i − 1 and still respecting the condition 5. 7. Repeat the steps 1 to 6, until no detector is selected in the step 5.
In this calibration procedure, only one detector with the highest measurements deviation is recalibrated at each iteration step. This careful approach allows to mitigate the risk of improper detector calibration correction in the cases where the sensitivity of few neighboring channels is commonly overestimated (or commonly underestimated).
Obviously, the success of this method relies on the fact that the data variability σ f is smaller than the measurements deviation m f,i − 1 . Two critical elements for this purpose are the quality of the measurements database (high statistics and diversity of emissivity profiles) and the accuracy of the tomographic reconstructions. As it was demonstrated in the last sections, the reconstruction accuracy strongly depends on the proper choice of the regularization parameter λ thanks to the discrepancy principle χ 2 N (λ discr ) σ 2 noise . However, it should be noted that the optimal error tolerance σ noise will progressively decrease during the correction procedure, as the global detectors cross-calibration is expected to improve after each iteration. It is therefore proposed to overcome this issue by applying a relaxation of the σ noise value. Initially, the iterative process is run with a relatively high σ noise value. After the end of step 7 (when no detector can be corrected anymore), the whole process is repeated with a σ noise value decreased by 10%. The σ noise is therefore decreased progressively, until the tolerance error reaches a lower threshold (e.g., the level of statistical noise, here 1%).
As expected, it can be seen in Fig. 12 that the calibration of edge channels, where the signal intensity is the lowest, is the most challenging task. The evolution of the mean absolute calibration error after each iteration of the correction procedure is shown in Fig. 13, and the final set of correction coefficients is presented in Fig. 14. Although relatively successful, it appears that the cross-calibration procedure is less accurate for the edge channels (as discussed before), but also for the very core channels. The latter observation can be explained by the fact that the regularization procedure slightly oversmooths steep emissivity gradients in the core, affecting the tomographic reconstructions of the peaked and hollow profiles, as seen in Fig. 11.
Overall, the lack of diversity in the chosen set of emissivity profiles can partly explain the limitations observed for the correction of core and edge channels. In a real tokamak environment, this issue could be solved by including plasma discharges in which the plasma center is vertically oscillating, in order to increase the spatial coverage of the lines-of-sight. Upgrading the tomography algorithm, e.g., by including the magnetic equilibrium as a priori information in the reconstruction [23], should also lead to a significant improvement.

Detectors position correction
A similar procedure as described in the previous section is applied here in order to correct mispositioning of the detectors. A position error shift z LoS is applied to each channel with a Gaussian distribution where σ 0.1. In this case, a scan of the position of each detector is performed, in order to find the detector position that matches the best the measured signal. An example of position scan is presented in Fig. 15, where it is clearly visible that the optimal position shift is associated with the lowest data variability. The evolution of the mean absolute position error after each iteration of the correction procedure is shown in Fig. 16 and the final corrective set of position shifts is presented in Fig. 17. Interestingly, it can be noted that the global position error increases after the iteration #38, probably indicating an issue of measurements overfitting when the tolerance error becomes too low. Similarly as for the detector sensitivity cross-calibration, the method could benefit from an increase of the diversity of profiles in the database and from an upgrade of the tomography algorithm.
In practice, the sensitivity and position cross-calibration of the detectors should be performed simultaneously in a real tokamak environment. An example of such calibration procedure can be found for the SXR diagnostic of ASDEX Upgrade in [20].

Application to a more realistic 2D scenario
The goal of this section is to introduce the basic steps necessary to build a synthetic SXR tomography diagnostics in a tokamak environment, together with some tools to assess the capabilities of the 2D tomography algorithm.

Definition of the tokamak geometry
First of all, a magnetic equilibrium should be defined in order to develop a reference plasma scenario. The equilibrium can be taken from a reconstruction numerical code, experimentally or by modeling [25]. In the absence of these tools, we will use a simpler solution, namely the Soloviev equilibrium as an analytical solution of the Grad-Shafranov equation [26], where the magnetic flux ψ(R, Z ) is expressed as: with R 0 0.5m the major radius and c 0 B 0 / R 2 0 κ 0 q 0 , c 1 B 0 κ 2 0 + 1 / R 2 0 κ 0 q 0 and c 2 0 are constant, where we will choose B 0 2T for the magnetic field, κ 0 1.2 for the ellipticity and q 0 1.4 for the safety factor on the axis.
The plasma boundary (r b , z b ) is defined by the following equation [26]: where ψ b 0.049 here. The normalized magnetic flux ψ N is such that: with ψ 0 the magnetic flux on the axis. The resulting equilibrium is plotted in Fig. 18, where the plasma is additionally horizontally shifted by R shift 0.16 m. In the following, the normalized minor radius ρ √ ψ N will be used to label the magnetic flux surfaces.

Definition of the detection system
The most commonly used detection system for SXR tomography in tokamaks is based on silicon barrier diodes (SBD), though other technologies based on gaseous detectors such as Gas Electron Multipliers (GEM) and Low Voltage Ionization Chambers (LVIC) are under development [27][28][29]. The X-ray absorption efficiency of a SBD can be estimated by: with α the absorption coefficient of silicon in cm 2 /g, that can be found for instance in the NIST XCOM database [30] and where we will choose the depletion zone V 3µm, the diffusion length L p 200µm and the thickness of the substrate D 380µm, see [31]. A 50-µm-thick beryllium pinhole is added in the system, almost absolutely cutting the SXR spectrum below 1 keV. The resulting detector spectral response η S X R is displayed in Fig. 19 (left).
A set of five SXR cameras, each composed of 20 photodiodes, is introduced here as presented in Fig. 19 (right). Each line-of-sight is defined by the positions of the photodiode and the associated pinhole. The cameras are placed on the outer board of the vacuum chamber, to account for the lack of available space that usually occurs in tokamaks.  The SXR plasma emissivity ε(hν) in W.m −3 .eV −1 of the bulk plasma, usually composed of fully ionized species such as hydrogen, deuterium, tritium or helium, comes from the electron-ion (free-free) collisions and can be described by the Bremsstrahlung formula [32]: where n e and T e are the plasma electron density and temperature, respectively, Z eff i n i Z 2 i /n e denotes the plasma effective charge (Z eff 1.2 here) and G ff is a free-free Gaunt factor, usually close to unity in the SXR range hν [1 keV; 20 keV].
The presence of nonfully ionized impurities in the plasma, such as argon, iron or tungsten ions, can modify significantly the SXR emissivity, due to complex contributions from line radiation (bound-bound) and radiative recombination (free-bound). Besides, it requires solving the ionization equilibrium of each species. For simplicity, this aspect is kept out of the scope of this paper. Nevertheless, it can be noted that many valuable spectroscopic data are periodically released by the Open-ADAS project and can be used to implement impurity radiation and transport in numerical codes [33].
Equation (27) requires to define electron temperature and density radial profiles in order to determine the radiation power spectrum, for instance based on an integrated tokamak simulator [34]. Here, we will mimic the shape of usual tokamak profiles by using Eq. (20) of regular and flat profiles with T e (0) 3 keV and n e (0) 5 × 10 19 m −3 respectively, as presented in Fig. 20 (left). The associated profile of the radiation power spectrum is plotted in Fig. 21. By integrating over all photon energies and convoluting with the photodiodes spectral response, the total and SXR-filtered radiation profiles are obtained and plotted in Fig. 20 (right). It can be noted that the SXR detectors can roughly detect half of the total radiated power in this plasma scenario.
It is then possible to remap the radiation profile onto the magnetic equilibrium in order to obtain a 2D synthetic emissivity field on a fine mesh (100×100 resolution here), as depicted in Fig. 22. The associated line-integrated measurements are obtained by convoluting the emissivity field with the transfer matrix of the detection system, see Eq. (6), as shown in Fig. 23.
All the necessary tools have now been developed to investigate the tomographic capabilities of the defined detection system. Some tomographic tests are shown in the next sections, Page 21 of 30 706

First tomographic test
A first tomographic test (using MFI) of the emissivity distribution obtained in the previous section is performed and presented in Fig. 24, including 1% of Gaussian noise in the measurements. It can be seen that the geometrical coverage of the plasma allows a rather satisfying reconstruction of the emissivity field, where the position and intensity of the maximum are properly recovered, as well as the general shape of the distribution. The reconstruction error map displayed in Fig. 24 (right) allows to observe that the reconstruction error is close to zero in the plasma center, while it reaches values up to 15% at the plasma edge due to the difficulty to reconstruct accurately steep emissivity gradients.

Tomographic reconstructions versus the number of cameras
Since the available space for diagnostics port-plugs is usually quite limited (as well as the budget), studying the optimal number and position of cameras for the purpose of tomographic reconstructions should be a mandatory step in the diagnostic design. Here, we show how the Eur. Phys. J. Plus (2021) 136:706

Fig. 22
Simulated soft X-ray emissivity profile in the defined geometry number of cameras affects the quality of the reconstruction in our defined scenario. The obtained tomograms are presented in Fig. 25. It can be seen that the number of cameras impacts dramatically the quality of the reconstruction. At least two or three cameras with a different angle of view are necessary to recover the global shape of the emissivity distribution.

Reconstruction vs. noise
The robustness of the reconstruction against experimental uncertainties is one of the critical aspects to consider for a tomographic system, as already discussed in the previous sections.
In tokamaks, the level of noise in line-integrated measurements can largely vary from less  Tomographic reconstruction test for different numbers of camera used in the reconstruction than 1% to more than 10% depending on the energy detection range (soft x-ray, hard xray or gamma-ray), the photon statistics, the detector technology, the electronic acquisition system or the electromagnetic environment surrounding the diagnostic. Here, we perform tomographic tests for four different levels, i.e., 2%, 5%, 10% and 20%, of zero-mean Gaussian noise in the measurements as depicted in Fig. 26. The associated tomograms are presented in Fig. 27, showing a relatively good resilience against noise with a reconstruction error growing, for 5% of noise, from 15 to 20% at the plasma edge and from 1 to 5% in the plasma core. Even for 20% of noise, the emissivity error in the plasma core remains below 10-15%, despite the presence of significant artifacts around the last closed flux surface.

Anisotropic regularization: a priori information from the magnetic equilibrium
In tokamaks, the plasma emissivity is usually expected to exhibit steeper gradients in the direction perpendicular to the magnetic flux surfaces (i.e., radially), due to reduced particle and heat transport, than in the parallel ones (i.e., toroidal and poloidal directions). Using this a priori information on the emissivity distribution could therefore significantly help finding an appropriate solution for the tomographic problem. In practice, the knowledge about the where ∇ // , ∇ ⊥ denote the parallel and perpendicular components of the gradient, respectively, as depicted in Fig. 28. The anisotropy factor τ < 0.5 allows to impose a higher level of smoothness of the solution in the parallel direction.
To demonstrate the beneficial effect of anisotropic regularization in our simple scenario, a series of tomographic reconstruction tests, with 2% of Gaussian noise in the measurements, is presented in Fig. 29 for different τ values from τ 0.5 (equivalent to isotropic regularization) up to τ 0.01. It is clearly visible that adding the ME constraint reduces the presence of artifacts outside of last closed surface and diminishes the emissivity error in the plasma center. The ME constraint can also be useful to reconstruct more accurately poloidal asymmetries in the emissivity distribution, as discussed in the next section.
Nevertheless, including the ME constraint implies adding a free parameter τ as well as an additional source of error in the reconstruction, as any error in ME reconstruction (horizontal or vertical shift, distortions) could propagate in the tomographic inversion. To illustrate this aspect, a series of tomographic reconstruction tests with anisotropic regularization (τ 0.05) is presented in Fig. 30, where a horizontal shift of ME to the low-field side by R shift 2cm up to R shift 20 cm was applied. As a result, while a ME shift lower than 5 cm does not seem to impact dramatically the quality of the solution, errors of the level or larger than 10 cm appear to be detrimental for the tomographic reconstruction in the presented scenario.

Asymmetric emissivity profiles
It was assumed in the previous sections that the emissivity was homogeneous on each magnetic flux surface. This assumption relies on the fact that the parallel transport is few orders of magnitude faster than the perpendicular transport, such that electron temperature and density distribution is homogenized on each flux surface. This consideration is usually correct in plasmas with low-Z impurities; however, heavy impurities like tungsten can be subject to centrifugal [35], electrostatic [36] or friction forces [37] that induce an asymmetry of their poloidal distribution. For instance, Reinke et al. derived an analytical formula for the poloidal asymmetry of the impurity density n S induced by the centrifugal force and by ion cyclotron resonance heating (ICRH): T e Z e f f T e + T i η(r ) , (29) where n S n S +ñ S , ω 2 S is the rotation velocity and η(r ) a factor accounting for the temperature anisotropy ratio (T ⊥ /T ) of the minority species. From this expression, we will test the tomography algorithm by simply assuming that the centrifugal asymmetry take the following general form for the emissivity: with the coefficient K 1 here, leading to a low-field side (LFS) asymmetry. The asymmetric emissivity profile and its reconstruction, using the full set of five cameras, are presented in Fig. 31.

Fig. 30
Tomographic reconstruction tests including the ME constraint (τ 0.05) with a horizontal shift to the low-field side: R shift 2cm, R shift 5cm, R shift 10cm and R shift 20cm (from left to right).

Fig. 31
Tomographic reconstruction test in the case of a low-field side poloidal asymmetry and 1% of Gaussian noise, with the initial emissivity (top left), the reconstructed emissivity without ME constraint (top center), with ME constraint and τ 0.05 (top right) and the associated measurements (bottom) The high-field-side asymmetry that can be induced by a change of the electrostatic potential in the presence of off-axis ICRH is mimicked by: ε ε −η(r )r cos θ.
As a result, it can be observed that the reconstructed tomograms recover the general shape of the asymmetry as well as the radial shift of the emissivity maximum. Though the Fig. 32 Tomographic reconstruction test in the case of a high-field side poloidal asymmetry and 1% of Gaussian noise, with: the initial emissivity (left), the reconstructed emissivity without ME constraint (top center), with ME constraint and τ 0.05 (top right) and the associated measurements (bottom) fine details of the HFS asymmetry are smoothed out by the isotropic regularization (middle plots), the anisotropic regularization (right plots) seems to reconstruct more accurately the typical banana shape.

Conclusion
In this paper, we have introduced some methodology and tools to implement a tomography algorithm for fusion devices. Based on a simple 1D tomography problem, the Tikhonov regularization has been described in detail with a focus on the second-order Philips-Tikhonov regularization and the Minimum Fisher information method. The optimal reconstruction parameters have been studied, such as the proper choice of the emissivity spatial resolution and the regularization parameter, thanks to the L-curve corner method or the discrepancy principle. A methodology has been proposed to perform an in situ sensitivity and position cross-calibration of the detectors with an iterative approach, by using the information redundancy and data variability in a given set of reconstructed profiles. The method has been validated on a set of synthetic data and its limitations have been discussed. Finally, a more realistic 2D synthetic X-ray tomography diagnostic has been introduced and described, including the detector response and the plasma scenario. The possibility to include the magnetic equilibrium constraint with an anisotropic regularization operator has been investigated. The synthetic diagnostic has been tested in several situations, including the presence of poloidal asymmetries in the emissivity distribution, illustrating the capabilities and limitations of a tomography algorithm to recover specific emissivity patterns in a given geometry.