A new pseudolocal tomography algorithm is developed for soft X-ray(SXR) imaging measurements of the turbulent electron temperature fluctuations (δ T_{e}) in tokamaks and stellarators. The algorithm overcomes the constraints of limited viewing ports on the vessel wall (viewing angle) and limited number of lines of sight (LOS). This is accomplished by increasing the number of LOS locally in a region of interest. Numerical modeling demonstrates that the wavenumber spectrum of the turbulence can be reliably reconstructed, with an acceptable number of viewing angles and LOS and suitable low SNR detectors. We conclude that a SXR imaging diagnostic for measurements of turbulent δ T_{e} using a pseudolocal reconstruction algorithm is feasible.

## I. INTRODUCTION

Better understanding of turbulence-driven transport is vital to improving plasma confinement in tokamaks and requires measurements of turbulence-related quantities. While there are many electron density fluctuations *δn*_{e} diagnostics available today,^{1–3} there are fewer electron temperature fluctuations *δT*_{e} diagnostics. Correlation electron cyclotron emission (CECE) is currently the only diagnostic capable of measuring *δT*_{e} in the plasma core and is limited to low-k (*k*_{⊥}*ρ*_{s} ≲ 1), and cannot usually be used in spherical tokamaks because of frequency cutoffs.^{4} High-k (*k*_{⊥}*ρ*_{s} > 1) *δT*_{e} diagnostics are important to understanding electron transport and would be a new and important tool for cutting-edge magnetic fusion research.

In this paper, we propose a new temperature fluctuation, *δT*_{e}, diagnostic using soft x-ray imaging (SXR). A pseudolocal tomography algorithm, a concept that has not been applied in fusion research before, is used to reconstruct local *δT*_{e} from the measurements of line-integrated soft x-ray emissivity. We describe the configuration of a feasible diagnostic using the new algorithm and explore capabilities of a new diagnostic numerically. We create a model for the reconstructed wavenumber-frequency spectrum, which is what a diagnostic would measure. We show that the measurement can be optimized by varying the number of viewing angles of x-ray detectors and the number of lines of sight (LOSs) at each viewing angle.

Soft x-ray imaging is a well-known diagnostic in magnetic fusion research, commonly used for two applications: (1) measuring electron temperature^{5,6} and (2) tomographically reconstructing the dynamic 2D emissivity profiles to study MHD phenomena.^{7} A number of tomographic reconstruction algorithms have been developed in fusion, including Cormack method,^{8} maximum entropy,^{6,9} linear regularization,^{6,10} minimum Fisher information,^{6,11} Gaussian process tomography,^{12} and deep learning.^{13} These algorithms successfully overcome the limitation of sparse measurements in the projection space,^{8} remove the high-frequency noise artifacts, and provide a smooth profile of emissivity or electron temperature. However, turbulent fluctuations are much smaller in amplitude compared to the noise profile and cannot be recovered by these algorithms. New tomographic reconstruction algorithms are required to enable measurement of fluctuations.

The major difference between measurement requirements for the *T*_{e} profile and turbulent *δT*_{e} is the different spatial and temporal scales. *δT*_{e} is at a much smaller scale both in space and time. Theoretically, if there were an unlimited number of viewing angles and LOSs and high enough sampling rate of x-ray detectors, then it should be possible to measure *δT*_{e}. The temporal resolution achieved in today’s tokamaks is about several 100 kHz,^{8,14} which is sufficient for the measurement of low-k *δT*_{e}, and is close to the range for high-k *δT*_{e}. With the advancement of instruments for radiation detection, the temporal resolution will eventually meet the requirement for a high-k *δT*_{e} diagnostic. Therefore, in this work, we assume high temporal resolution and focus on how to improve the spatial resolution subject to the constraints of the limited number of viewing angles and LOSs.

## II. PSEUDOLOCAL TOMOGRAPHY

This paper presents a new tomographic algorithm called pseudolocal tomography customized for SXR diagnostics in tokamaks. Pseudolocal tomography has been used in other fields but not fusion.^{15} Unlike global tomography, pseudolocal tomography uses both global measurements and measurements localized around a particular region of interest (ROI). Figure 1(a) shows the configuration of the diagnostic at any viewing angle *θ*. Instead of a fan-beam configuration, we adopt a parallel-beam configuration to illustrate our new algorithm. In Fig. 1(a), the black circle is the poloidal section of a tokamak, the X-ray detectors are assembled in the camera (orange strip), the yellow circle is the ROI for measuring *δT*_{e}, the blue dense parallel segments are the localized LOS and the red sparse parallel segments are the global LOS. The dense LOS is what distinguishes pseudolocal tomography from global tomography (where all LOS are evenly spaced). The addition of the dense LOS (blue) is needed to obtain the information about *δT*_{e} in the ROE. Outside the ROI, the spacing of LOS is large (red) because they are only used to reconstruct the global *T*_{e} profile. By reducing LOS, the diagnostic cost is reduced.

A qualitative explanation for how to reconstruct the turbulence is helpful. In conventional tomography with evenly spaced LOSs, the reconstruction can be done via the back projection formula^{16}

where $Rf\u0302(\theta ,k)=\u222bRf(\theta ,s)e\u2212i2\pi skds$ is the Fourier transform of *Rf* and *Rf* is the line-integration of *f* whose expression is $Rf(\theta ,s)=\u222bx\u20d7\u22c5\theta \u0302=sf(s\theta \u0302+t\theta \u0302\u22a5)dt$. *Rf*(*θ*, *s*) is fully determined by the viewing angle *θ* and LOS’s coordinate *s*. The two unit vectors are defined as $\theta \u0302=(cos\u2061\theta ,sin\u2061\theta )$ and $\theta \u0302\u22a5=(sin\u2061\theta ,\u2212cos\u2061\theta )$. When applying this algorithm in practice, a ramp filter is often added to filter out the high-frequency noise, giving a more stabilized reconstruction [thus, this algorithm is usually called filtered back projection (FBP)]. FBP is the primary method in tomographic reconstruction in many fields, like medicine, because of ease of implementation. However, in fusion research, FBP is seldom used due to irregular sampling and insufficient sampling density in (*θ*, *s*) space. Limited efforts have been made to apply FBP in fusion plasmas,^{17} showing that it is possible to get a satisfactory reconstruction when an appropriate interpolation scheme is implemented in projection space.

The new pseudolocal algorithm introduced here is a variation on FBP, with the addition of a new localized sampling option. The modified reconstruction formula for pseudolocal tomography becomes

where *f* is divided into two parts: the profile $f\u0304$ and the fluctuations $f\u0303$. The two parts are reconstructed separately with different k-resolutions. Figure 1(b) shows a typical wavenumber spectrum of the line-integrated signal $Rf\u0302(k)$ at a certain angle *θ*. The red curve corresponds to the global profile $f\u0304$, living in very low-k region $(0,kcutofflow)$ of the spectrum. To reconstruct this part, we need to cover the whole low-*k* region, and the spacing of LOSs in real space should be no bigger than $\Delta xglb=2\pi /kcutofflow$, matching the red LOS in Fig. 1(a). The requirement on the wavenumber spacing Δ*k*_{glb} is 2*π*/Δ*k*_{glb} > 2*a*, that is, covering the whole poloidal cross section. The blue curve corresponds to the fluctuating part $f\u0303$. It spans a wide range of wavenumbers $(kcutofflow,kcutoffhigh)$. To reconstruct this part, the spacing of LOS in real space should be no larger than $\Delta xloc=2\pi /(kcutoffhigh\u2212kcutofflow)\u22482\pi /kcutoffhigh\u226a\Delta xglb$, matching the blue LOSs in Fig. 1(a). The requirement for the wavenumber spacing Δ*k*_{loc} is that Δ*k*_{loc} should be small enough to resolve the spectrum of the fluctuations. However, Δ*k*_{loc} should also not be too small; otherwise, it requires too many local LOSs.

The separate reconstruction with different wavenumbers or spatial resolutions can be realized as long as the above requirements for Δ*x*_{glb}, Δ*k*_{glb}, Δ*x*_{loc}, Δ*k*_{loc} are fulfilled. It is the natural separation of the global profile and the local fluctuations in wavenumber space that allows for pseudolocal tomography to provide such separate reconstructions.

## III. NUMERICAL MODEL

The pseudolocal algorithm proposed needs to be verified computationally before being implemented. A simplified model is developed to test this algorithm. The setup of the model is illustrated still in Fig. 1(a). The poloidal cross section is a unit circle with radius *a* = 0.6 m, and the plasma surfaces are concentric circles. A camera carrying a number of multi-color x-ray detectors is installed at each of several viewing angles in the same poloidal cross section. The purple dots in the ROI are called measuring points because they are the points where we want to measure *δT*_{e} using the pseudolocal tomography. The arrangement of measuring points determines which direction the total wavenumber vector projects onto. In our study, we let measuring points align with the radius of the poloidal cross section so as to measure the radial wavenumber *k*_{r}.

This numerical modeling includes two steps. The first step is synthetic data generation, and the second step is the reconstruction using pseudolocal tomography. For the first step, we assume the profiles of *n*_{e} and *T*_{e} are in the form of $(1\u2212r2)\alpha $ and are a combination of the equilibrium profile and fluctuation profile. The default values are $\alpha ne=0.5$ and $\alpha Te=0.8$. The fluctuations are generated from the power spectrum,

where $kr\u0302=kr/krc,k\theta \u0302=k\theta /k\theta c,f\u0302=f/fc$ are normalized to the center radial wavenumber *k*_{rc}, center poloidal wavenumber *k*_{θc}, and center frequency *f*_{c} of the spectrum. This expression is an abstraction based on the results from a full-physics nonlinear gyrokinetic simulation.^{18} The fluctuations are added on top of the electron density and temperature profile, and the combined electron density and temperature are used to calculate the emissivity via the formula $\epsilon \u221dne2Tee\u2212EcTe$, where *E*_{c} is the cutoff energy determined by the thickness of the foil. At least two *E*_{c}s are required to infer *T*_{e} from the fit of *ɛ* against *E*_{c}. In our model, we have three *E*_{c} values in an arithmetic sequence. This emissivity function is simplified. It does not contain the contribution of discontinuous lines or the effects of impurities. The local emissivity is used to calculate the synthetic measurements, which are the brightness or the line-integration of emissivity along the given LOS specified by the coordinate (*θ*, *s*).

For the second step, generated synthetic measurements are used to reconstruct *δT*_{e}. Pseudolocal tomography is applied to reconstruct the local emissivity from the brightness. The local emissivity, with corresponding *E*_{c}s, is used to infer *T*_{e} and *δT*_{e} by eliminating the time-averaged component from *T*_{e}. Note that due to insufficient sampling in (*θ*, *s*) space, the reconstructed *δT*_{e} can be quite different from the synthetic *δT*_{e}, but the up-and-down pattern among the measuring points may be retained since the reconstruction error due to numerical dicretization is somewhat systematic for all measuring points. Instead of comparing *δT*_{e} in real space, the Fourier transform of the reconstructed *δT*_{e} is used to obtain the “measured” spectrum and compare with synthetic *δT*_{e} spectrum.

## IV. OPTIMIZING DIAGNOSTIC DESIGN

First we examine the number of viewing angles and LOS per viewing angle required for a reconstruction of synthetic spectrum. We then assess if this is feasible in practice. We scan the parameter space (*N*_{θ}, *N*_{glb}, *N*_{loc}, *N*_{mp}), where *N*_{θ} is the number of viewing angles, *N*_{glb} is the number of global LOSs, *N*_{loc} is the number of local LOSs and *N*_{mp} is the number of measuring points. We compare the reconstructed spectrum with the synthetic spectrum identify the best choices for diagnostic design. To do a practical scan of the 4D parameter space, we change one parameter at a time while keeping the other three fixed. Key parameters are the electron temperature and density at *r* = 0 are *T*_{e0} = 0.5 keV and *n*_{e0} = 5.34*e*19 m^{−3} respectively; local ion sound gyro-radius in the ROI is *ρ*_{s} = 0.007 m; the maximum fluctuation levels are *δn*_{e}/*n*_{e} = *δT*_{e}/*T*_{e} = 5%; central radial and poloidal wavenumber in the ROI for electron temperature and electron density $krc,Te=krc,ne=k\theta c,Te=k\theta c,ne=\rho s\u22121$; spacing of local LOSs is Δ*x*_{loc} = *λ*_{rc,Te}/8 = (2*π*/*k*_{rc,Te})/8; spacing of global LOSs is 2*a*/*N*_{glb}; central frequencies are *f*_{c,Te} = *f*_{c,ne} = 1 MHz; sampling rate is 6 MHz; radial center of the ROI *r*_{0}/*a* = 0.6; energy threshold of foils *E*_{c}/*T*_{e0} = (0.5, 1.0, 1.5). We assume that SNR of the measurement does not impact geometric configuration of the diagnostic so we set the SNRs for all LOSs to be infinite. The effect of noise is studied independently in the Sec. V.

Figure 2(a) shows an example of the synthetic spectrum (a) and reconstructed spectrum (b) for the diagnostic parameter number of viewing angles *N*_{θ} = 10, number of global LOSs *N*_{glb} = 5, number of measuring points *N*_{mp} = 21, and number of local LOSs *N*_{loc} = 21. They are obtained by smoothing the raw spectrum data to remove the random noise. The smoothing is done by Fourier transforming the raw spectrum. The Fourier transform of the spectrum of the fluctuations is localized in the low frequency region, while the Fourier transform of the spectrum of the noise is an evenly distributed background. Removing the evenly distributed background from the total Fourier transform and then inverse Fourier transforming it give us the smoothed spectrum. The red stars indicate the center (or peak) location (*k*_{rc}, *f*_{c}) of the spectra. The red circles represent the contour lines whose level is 90% of the peak value. While the shape of the spectrum defined by the contour line is not sensitively affected by the diagnostic parameters, the center coordinate *k*_{rc} is rather sensitive to the diagnostic parameters. Therefore, we can look at the relative error of *k*_{rc}s between the spectra to quantify the quality of the reconstruction. The definition of relative error is

Figure 3 shows *ϵ*_{k} plotted against four diagnostic parameters (*N*_{θ}, *N*_{glb}, *N*_{loc}, *N*_{mp}). Each data point is obtained by running the same set of parameters for 20 times and taking the average *ϵ*_{k} over these 20 runs. This reduces statistical uncertainties. We first study the relationship between *ϵ*_{k} and *N*_{θ}; other parameters are empirically chosen to be sufficiently large (*N*_{glb}, *N*_{mp}, *N*_{loc}) = (19, 21, 39). From the curve of *ϵ*_{k} ∼ *N*_{θ}, we can find that *N*_{θ} = 10 is optimal since our goal is to minimize *ϵ*_{k} with as few viewing angles as possible. One interesting observation is that *ϵ*_{k} does not decrease monotonously with *N*_{θ} as we might expect. This has to do with the numerical implementation of the algorithm. One possible explanation is the convergence of signed *ϵ*_{k} is oscillatory. There are many zeros of the signed *ϵ*_{k} ∼ *N*_{θ} curve, and we just take the first zero crossing.

Next, we examine the relationship between *ϵ*_{k} and *N*_{glb}. Since we already found the optimal *N*_{θ} = 10, we keep this value fixed in the remainder of the analysis and fix (*N*_{loc}, *N*_{mp}) = (39, 21). The minimum *ϵ*_{k} is achieved at *N*_{glb} = 13. However, if *ϵ*_{k} < 5% is acceptable, we can choose *N*_{glb} to be 5. This is important from a practical standpoint since the required number of global LOSs is at least halved.

Next, we examine the impact of *N*_{mp}, the number of measuring points. Because the local spacing is always set to be 1/8 of *λ*_{c,Te}, the central wavelength of *δT*_{e}, changing *N*_{mp} is equivalent to changing the number of wavelengths measured. Using more *N*_{mp} means more information about the fluctuations, aiding in reconstruction. For this scan, we do not fix *N*_{loc} but instead fix the ratio *N*_{loc}/*N*_{mp} = 2, since all measuring points are surrounded by some number of local LOSs. The purple curve in Fig. 3 shows how *ϵ*_{k} changes with *N*_{mp} and 25 is the optimal value of *N*_{mp}, but *N*_{mp} = 21 can also be chosen based on the *ϵ*_{k} < 5% standard.

Finally, *N*_{loc}, the number of local LOSs, is changed, while other parameters are fixed at their optimal values (*N*_{θ}, *N*_{glb}, *N*_{mp}) = (10, 5, 21). We can see *ϵ*_{k} does not change much with *N*_{loc} and is always below 5%. This allows choosing the minimum *N*_{loc} = *N*_{mp}. Now, we find the optimal setting of all the four parameters, which is (*N*_{θ}, *N*_{glb}, *N*_{mp}, *N*_{loc}) = (10, 5, 21, 21), and the total number of LOSs is *N*_{total} = *N*_{θ} · (*N*_{glb} + *N*_{loc}) = 260. This diagnostic layout is a large number of views but is consistent with many existing systems, for example, the tokamak à configuration variable (TCV) is equipped with ten pinhole cameras and a total of 200 LOSs.^{6} The stellarator W7-X has 20 pinhole cameras with a total of 360 LOSs.^{14} We conclude that the proposed diagnostic setup is achievable.

## V. MODELING THE EFFECTS OF SIGNAL-TO-NOISE RATIO (SNR)

In this section, the effect of signal-to-noise ratio (SNR) on the reconstructed spectrum is studied. The diagnostic parameters are set to the optimal values identified in Sec. IV: (*N*_{θ}, *N*_{glb}, *N*_{mp}, *N*_{loc}) = (10, 5, 21, 21), and we now vary SNR. The SNR depends on the intensity of both signal and noise. The noise can be modeled as the sum of two types of noises, a Poisson-type noise due to the statistical nature of x-ray emission and a Gaussian-type of additive electronic noise.^{19} We model these two types of noises separately to see their individual effects. For the Poisson-type noise, the SNR is proportional to the square root of signal (brightness),^{20} and the additive Gaussian has SNR ∝ brightness. Note that the brightness in our model depends on the LOS and also the energy threshold *E*_{c}, so each LOS has a different SNR for both types of noise. The following two effects lead to low SNR: low emissivity, such as in the edge LOS, and large *E*_{c}. In our model, the SNR for each LOS is calculated via

where *α* = 1/2 for Poisson noise and *α* = 1 for Gaussian noise, brightness_{m} is the maximum brightness received in all LOSs, SNR_{m} is the maximum SNR corresponding to brightness_{m}, which is the parameter we vary. For each LOS, the SNR is given by Eq. (5).

Figure 2(a) is the synthetic spectrum for the case we study, (b) is the reconstructed spectrum of *δT*_{e} without adding noises; (c) and (d) are both reconstructed for Poisson noise with SNR_{m} = 1000, but (c) does nothing to the noise, while (d) filters out the noise before the reconstruction. The reconstructed spectrum (b) basically recovers the synthetic spectrum (a), implying that the diagnostic works well for the noiseless case. The reconstructed spectrum in (c) is dominated by the noise in the background even if it is still possible to recognize the shape of the spectrum of *δT*_{e} similar to the ones shown in (a) and (b). This means our diagnostic fails if SNR_{m} is 1000 or less. However, if a low pass filter is added to remove the noise in the measured brightness, the reconstructed spectrum still recovers the synthetic spectrum to a great degree, shown in (d). The filter is designed based on the fact that the spectrum of *δT*_{e} is smooth (low frequency) and the spectrum of the noise is fluctuating (high frequency). Compared with (a) and (b), the spectrum in (d) is a little more elongated along the f-axis in shape and still has non-negligible components in the background. These effects are minor if the precision requirement is not very high. Thus, the addition of the filter makes the diagnostic work again for SNR_{m} = 1000. Without the filter, the least requirement for SNR_{m} is ∼1500, and this limit can be pushed to ∼700 using this filter.

The same study is carried out with Gaussian noise (not shown here). For Gaussian noise, the diagnostic without the denoising filter fails for SNR_{m} ≲ 2000, which is a higher requirement than the one for Poisson noise. This is reasonable since the SNR of Gaussian noise decreases more rapidly with the brightness, compared to Poisson noise. Adding the denoising filter, the least requirement for SNR_{m} is also reduced to ∼1500. When both Poisson noise and Gaussian noise are present, the SNR_{m} requirement should also be ≳1000. This SNR_{m} is completely acceptable, since a much higher SNR_{m} such as 5000 has been reported in previous work.^{14}

## VI. DISCUSSION

We conclude that a SXR imaging diagnostic for measurements of electron temperature fluctuations in a tokamak or stellarator is feasible, using a new algorithm we developed. One next step is to use gyrokinetic simulations to generate test data to provide more realistic turbulence characteristics. In addition, future work will consider a more realistic emissivity model that includes the effects of impurities.

## ACKNOWLEDGMENTS

This work is supported by US DOE contract DE-SC0019089. Computer simulations were carried out at the NERSC Center, supported by the Office of Science of the U.S. DOE under Contract No. DE-AC02- 05CH11231 and at the MIT-PSFC partition of the Engaging cluster at the MGHPCC facility (www.mghpcc.org) funded by DoE grant number DE-FG02-91-ER54109.

## DATA AVAILABILITY

Data for this study are available upon reasonable request.

## REFERENCES

*et al.*,

*et al.*,