A Finite-Difference Approach for Plasma Microwave Imaging Profilometry

Abstract

Plasma diagnostics is a topic of great interest in the physics and engineering community because the monitoring of plasma parameters plays a fundamental role in the development and optimization of plasma reactors. Towards this aim, microwave diagnostics, such as reflectometric, interferometric, and polarimetric techniques, can represent effective means. Besides the above, microwave imaging profilometry (MIP) may allow the obtaining of tomographic, i.e., volumetric, information of plasma that could overcome some intrinsic limitations of the standard non-invasive diagnostic approaches. However, pursuing MIP is not an easy task due to plasma’s electromagnetic features, which strongly depend on the working frequency, angle of incidence, polarization, etc., as well as on the need for making diagnostics in both large (meter-sized) and small (centimeter-sized) reactors. Furthermore, these latter represent extremely harsh environments, wherein different systems and equipment need to coexist to guarantee their functionality. Specifically, MIP entails solution of an inverse scattering problem, which is non-linear and ill-posed, and, in addition, in the one-dimensional case, is also severely limited in terms of achievable reconstruction accuracy and resolution. In this contribution, we address microwave inverse profiling of plasma assuming a high-frequency probing regime when magnetically confined plasma can be approximated as both an isotropic and weak penetrable medium. To this aim, we adopt a finite-difference frequency-domain (FDFD) formulation which allows dealing with non-homogeneous backgrounds introduced by unavoidable presence of plasma reactors.

1. Introduction

Plasma diagnostics using active microwave techniques is one of the most promising tools to optimize the production of plasma in both large and compact reactors because of its non-invasive nature and relatively modest access requirements [1]. Among the others, microwave imaging reflectometry (MIR), interferometric and polarimetric techniques are currently under investigation to extract plasma proprieties useful to model and optimize the heating process. Such information is crucial to study and improve plasma generation based on Electron Cyclotron Resonance [2] as well as innovative heating schemes such as the ones employing Electrostatic Bernstein Waves in overdense plasma [3]. Unfortunately, the diagnostic approaches mentioned above show some limitations. As an example, MIR, a multifrequency radar-like technique, mainly used in large reactors, such as tokamak and stellarator, is commonly used to infer electronic density distribution and fluctuation, but suffers from the impossibility of performing diagnostics of non-increasing plasma profiles since inner reflecting cutoff layers are “hidden” from the outer ones [4]. On the other hand, line-integrated measurements, such as those actually provided by interferometric and polarimetric techniques, although employable on compact reactors such as the electron cyclotron ion source (ECRIS) [5], cannot give local information about electron density, but only a line-of-sight integrated value.

For these reasons, it is worth devoting interest towards microwave imaging tomographic approaches, such as those actually used in medical and subsurface microwave diagnostics [6], which are in principle able to “locally” monitor plasma uniformity and instabilities.

In a recent paper, we investigated the possibility of obtaining plasma imaging profilometry (MIP) by means of electromagnetic inverse scattering techniques requiring only measurements of the reflection coefficient in a quite large frequency band. In particular, adopting a high-frequency probing regime (i.e., probing frequencies much larger than the cyclotron frequency), i.e., when the plasma can be assumed to be both isotropic and weak scatterer, linear imaging recovery techniques are able to achieve quantitative characterization of plasma profiles [7]. The results have shown that if the recovery process is properly formulated as a sparse optimization problem, compressive sensing (CS)-based strategies are able to achieve nearly optimal results. However, free-space homogeneous background has been considered, which is an assumption suitable only if multiple reflections introduced by the metallic environment (reactor walls) are neglected or, in any case, filtered out.

With the aim to get a further step in the analysis of such a problem, in this paper we introduce a frequency-domain finite-difference (FDFD) formulation for the solution of the inverse scattering problem in order to take easily into account the effect of non-homogeneous backgrounds related to the metallic environment of the plasma chamber.

On the other hand, besides the modeling aspects, solution of the inverse scattering problem in presence of a metallic surface, or, in general, non-homogeneous background, can bring some advantages due to the possibility to achieve enhanced resolution in the reconstruction process [8,9]. Roughly speaking, this enhancement can be explained due to the multiple reflected waves which can bring additional information conveyed back by the recovery process. Interestingly, this may be of fundamental importance in the solution of the one-dimensional inverse scattering problem (also known as inverse profiling) wherein the number and kind of retrievable information is severely limited [10].

The paper is structured as follows. In Section 2, the adopted mathematical formulation for the solution of the forward and inverse scattering problem is described. In Section 3 the microwave imaging profilometry for inhomogeneous plasma slabs recovery is formulated. Finally, in Section 4, the developed approach is validated against simulated data dealing with homogeneous and non-homogeneous backgrounds as well as reflection-only and reflection-transmission measurement configuration. Conclusions follow.

2. Mathematical Formulation

Let us assume a bounded, simply connected, domain D located in an non-homogeneous background medium, wherein some electromagnetic primary sources $\underline{\mathbf{j}}{}_0$ are located in a region characterized by relative permittivity and permeability distribution ${\mu }_r\left(\underline{\mathbf{r}}{}^{\prime }\right)$ and ${ϵ}_r\left(\underline{\mathbf{r}}{}^{\prime }\right)$, $\underline{\mathbf{r}}{}^{\prime }$ spanning the space D. After normalizing the electric field, magnetic field and the primary sources according to $\widetilde{\underline{\mathbf{e}}}=\sqrt{\frac{{ϵ}_0}{{\mu }_0}}\underline{\mathbf{e}}$, $\widetilde{\underline{\mathbf{h}}}=\underline{\mathbf{h}}$ and ${\widetilde{\underline{\mathbf{j}}}}_0=c_0\underline{\mathbf{j}}{}_0$, with ${\displaystyle c_0=\frac{1}{\sqrt{{ϵ}_0{\mu }_0}}}$ the speed of light, Maxwell equations read:

$$\begin{array}{c}{c}_0\underline{\nabla }\times \widetilde{\underline{\mathbf{e}}}=-\mathrm{j}\omega {\mu }_r\left(\underline{\mathbf{r}}{}^{\prime }\right)\widetilde{\underline{\mathbf{h}}}\\ c_0\underline{\nabla }\times \widetilde{\underline{\mathbf{h}}}=\mathrm{j}\omega {ϵ}_r\left(\underline{\mathbf{r}}{}^{\prime }\right)\widetilde{\underline{\mathbf{e}}}+{\widetilde{\underline{\mathbf{j}}}}_0\end{array}$$

(1)

where $\omega$ is the angular frequency and the time factor $exp\left\{\mathrm{j}\omega t\right\}$ is assumed and dropped for the time-harmonic fields. Normalized electric and magnetic fields, $\widetilde{\underline{\mathbf{e}}}$ and $\widetilde{\underline{\mathbf{h}}}$, are of the same order of magnitude and this is an advantage in formulating the perfectly matched layer (PML) [11]. If, for the sake of a simpler notation, we drop the tilde $\tilde{}$ symbol and the $c_0$ is included in the $\underline{\nabla }\times$ definition, we can write:

$$\begin{array}{c}\underline{\tilde{\nabla }}\times \underline{\mathbf{e}}=-\mathrm{j}\omega {\mu }_r\left(\underline{\mathbf{r}}{}^{\prime }\right)\underline{\mathbf{h}}\\ \underline{\tilde{\nabla }}\times \underline{\mathbf{h}}=\mathrm{j}\omega {ϵ}_r\left(\underline{\mathbf{r}}{}^{\prime }\right)\underline{\mathbf{e}}+\underline{\mathbf{j}}{}_0\end{array}$$

(2)

In the above formulation, to satisfy different boundary conditions, a lossy PML [12] can be introduced to truncate the computational domain, as commonly done in numerical solver. By choosing proper PML parameters, various boundary conditions can seamlessly introducing ranging from a perfect absorbing boundary condition to a PEC (or PMC) boundary condition, the latter case when no PML is used [11].

The background field (i.e., the field without the unknown profile) defined as a plane wave propagating in a non-homogeneous, eventually lossy, background with permittivity and permeability, respectively, given by:

$${ϵ}_b\left(\underline{\mathbf{r}}{}^{\prime }\right)={ϵ}_b^{\prime }\left(\underline{\mathbf{r}}{}^{\prime }\right)-j{ϵ}_b^{″}\left(\underline{\mathbf{r}}{}^{\prime }\right)={ϵ}_b^{\prime }\left(\underline{\mathbf{r}}{}^{\prime }\right)-j\frac{\sigma \left(\underline{\mathbf{r}}{}^{\prime }\right)}{{ϵ}_0\omega }$$

$${\mu }_b\left(\underline{\mathbf{r}}{}^{\prime }\right)={\mu }_b^{\prime }\left(\underline{\mathbf{r}}{}^{\prime }\right)-j{\mu }_b^{″}\left(\underline{\mathbf{r}}{}^{\prime }\right)$$

satisfies the following equations:

$$\begin{array}{c}\underline{\tilde{\nabla }}\times \underline{\mathbf{e}}{}_{bck}=-\mathrm{j}\omega {\mu }_b\left(\underline{\mathbf{r}}{}^{\prime }\right)\underline{\mathbf{h}}{}_{bck}\\ \underline{\tilde{\nabla }}\times \underline{\mathbf{h}}{}_{bck}=\mathrm{j}\omega {ϵ}_b\left(\underline{\mathbf{r}}{}^{\prime }\right)\underline{\mathbf{e}}{}_{bck}+\underline{\mathbf{j}}{}_0\end{array}$$

(3)

To formulate the scattering problem, we can add and subtract the quantities $\mathrm{j}\omega {ϵ}_b\left(\underline{\mathbf{r}}{}^{\prime }\right)\underline{\mathbf{e}}$ and $\mathrm{j}\omega {\mu }_b\left(\underline{\mathbf{r}}{}^{\prime }\right)\underline{\mathbf{h}}$ at the right hand side member of Equation (2) to obtain:

$$\begin{array}{c}\underline{\tilde{\nabla }}\times \underline{\mathbf{e}}=-\mathrm{j}\omega {\mu }_b\left(\underline{\mathbf{r}}{}^{\prime }\right)\underline{\mathbf{h}}-\mathrm{j}\omega \left[{\mu }_r\left(\underline{\mathbf{r}}{}^{\prime }\right)-{\mu }_b\left(\underline{\mathbf{r}}{}^{\prime }\right)\right]\underline{\mathbf{h}}\\ \underline{\tilde{\nabla }}\times \underline{\mathbf{h}}=\mathrm{j}\omega {ϵ}_b\underline{\mathbf{e}}+\mathrm{j}\omega \left[{ϵ}_r\left(\underline{\mathbf{r}}{}^{\prime }\right)-{ϵ}_b\left(\underline{\mathbf{r}}{}^{\prime }\right)\right]\underline{\mathbf{e}}+\underline{\mathbf{j}}{}_0\end{array}$$

(4)

and, finally, subtracting (3) to (4), the equations governing the scattered fields, defined as $\underline{\mathbf{e}}{}_{sct}=\underline{\mathbf{e}}-\underline{\mathbf{e}}{}_{bck}$ and $\underline{\mathbf{h}}{}_{sct}=\underline{\mathbf{h}}-\underline{\mathbf{h}}{}_{bck}$ , are given by:

$$\begin{array}{c}\underline{\tilde{\nabla }}\times \underline{\mathbf{e}}{}_{sct}=-\mathrm{j}\omega {\mu }_b\left(\underline{\mathbf{r}}{}^{\prime }\right)\underline{\mathbf{h}}{}_{sct}-\mathrm{j}\omega {\mu }_b\left(\underline{\mathbf{r}}{}^{\prime }\right)\left[\frac{{\mu }_r\left(\underline{\mathbf{r}}{}^{\prime }\right)-{\mu }_b\left(\underline{\mathbf{r}}{}^{\prime }\right)}{{\mu }_b\left(\underline{\mathbf{r}}{}^{\prime }\right)}\right]\underline{\mathbf{h}}\\ \underline{\tilde{\nabla }}\times \underline{\mathbf{h}}{}_{sct}=\mathrm{j}\omega {ϵ}_b\left(\underline{\mathbf{r}}{}^{\prime }\right)\underline{\mathbf{e}}{}_{sct}+\mathrm{j}\omega {ϵ}_b\left(\underline{\mathbf{r}}{}^{\prime }\right)\left[\frac{{ϵ}_r\left(\underline{\mathbf{r}}{}^{\prime }\right)-{ϵ}_b\left(\underline{\mathbf{r}}{}^{\prime }\right)}{{ϵ}_b\left(\underline{\mathbf{r}}{}^{\prime }\right)}\right]\underline{\mathbf{e}}\end{array}$$

(5)

wherein the quantities under the square brackets can be easily recognized as the contrast functions:

$${\chi }_{\mu }(\omega ,\underline{\mathbf{r}}{}^{\prime })=\frac{{\mu }_r(\omega ,\underline{\mathbf{r}}{}^{\prime })-{\mu }_b(\omega ,\underline{\mathbf{r}}{}^{\prime })}{{\mu }_b(\omega ,\underline{\mathbf{r}}{}^{\prime })},$$

$${\chi }_{ϵ}(\omega ,\underline{\mathbf{r}}{}^{\prime })=\frac{{ϵ}_r(\omega ,\underline{\mathbf{r}}{}^{\prime })-{ϵ}_b(\omega ,\underline{\mathbf{r}}{}^{\prime })}{{ϵ}_b(\omega ,\underline{\mathbf{r}}{}^{\prime })},$$

which relate the dielectric and magnetic properties of the anomalies, i.e., scattering objects, to those of the embedding background medium at each frequency.

Adopting a dual grid [13,14] and a proper discretization, in operator notation, the above equations can be rewritten as:

$$\left[\begin{array}{cc}\mathcal{C}& \mathrm{j}\omega {\mathcal{M}}_b\\ -\mathrm{j}\omega {\mathcal{E}}_b& \mathcal{C}\end{array}\right]\left[\begin{array}{c}\underline{\mathbf{e}}{}_{sct}^{\prime }\\ \underline{\mathbf{h}}{}_{sct}^{\prime }\end{array}\right]=\left[\begin{array}{cc}0& -\mathrm{j}\omega {\mu }_b{\chi }_{\mu }\\ \mathrm{j}\omega {ϵ}_b{\chi }_{ϵ}& 0\end{array}\right]\left[\begin{array}{c}\underline{\mathbf{e}}{}^{\prime }\\ \underline{\mathbf{h}}{}^{\prime }\end{array}\right],$$

(6)

wherein ${\left(\right)}^{\prime }$ is used to indicate the field integrated over the appropriate line element of the primal or dual mesh ($\underline{\mathbf{e}}{}^{\prime }\approx \underline{\mathbf{e}}\overline{\Delta L}$, $\underline{\mathbf{h}}{}^{\prime }\approx \underline{\mathbf{h}}\widetilde{\Delta L}$). Moreover, the matrix operator $\mathcal{C}$ known in the literature as incidence matrix (see for example section §7.3.3 and Equation (7.1) of Tonti’s book [13] or [15]) denotes the proper derivative operations and is a sparse matrix. Moreover, ${\mathcal{M}}_b$ and ${\mathcal{E}}_b$ are matrices accounting for the magnetic and dielectric non-homogeneous background properties, respectively, wherein ${\mathcal{E}}_b\approx {ϵ}_b\frac{\widetilde{\Delta S}}{\overline{\Delta L}}$ and ${\mathcal{M}}_b\approx {\mu }_b\frac{\overline{\Delta S}}{\widetilde{\Delta L}},$ with $\overline{\Delta L}$ and $\overline{\Delta S}$ refer to lines and surfaces of the primal mesh, and $\widetilde{\Delta L}$ and $\widetilde{\Delta S}$ refer to the dual mesh [13]. In the remaining of the paper for the sake of a simpler notation we drop the ${}^{\prime }$ symbol.

The formal solution of the linear system (6) can be rewritten as:

$$\left[\begin{array}{c}\underline{\mathbf{e}}{}_{sct}\\ \underline{\mathbf{h}}{}_{sct}\end{array}\right]={\mathcal{H}}_b^{-1}\left(\left[\begin{array}{cc}0& -\mathrm{j}\omega {\mu }_b{\chi }_{\mu }\\ \mathrm{j}\omega {ϵ}_b{\chi }_{ϵ}& 0\end{array}\right]\left[\begin{array}{c}\underline{\mathbf{e}}\\ \underline{\mathbf{h}}\end{array}\right]\right).$$

(7)

The above equation is the key equation for the inverse scattering problem adopting a finite-difference formulation. As a matter of fact, the measured scattered fields, i.e., the scattered field recorded on the data domain S (for the case at hand, the same positions $z=z_S$ at different frequencies), can be represented by the following operator notations:

$$\left[\begin{array}{c}\underline{\mathbf{e}}{}_{sct}(z_S,{\omega }_n)\\ \underline{\mathbf{h}}{}_{sct}(z_S,{\omega }_n)\end{array}\right]={\mathcal{M}}_S\left\{{\mathcal{H}}_b^{-1}\left(\left[\begin{array}{cc}0& -\mathrm{j}\omega {\mu }_b{\chi }_{\mu }\\ \mathrm{j}\omega {ϵ}_b{\chi }_{ϵ}& 0\end{array}\right]\left[\begin{array}{c}\underline{\mathbf{e}}\\ \underline{\mathbf{h}}\end{array}\right]\right)\right\},$$

(8)

where ${\mathcal{M}}_S$ is an operator that extracts the field values at the observation point $z=z_S$ at different frequencies. Substituting ${[\underline{\mathbf{e}}{}_{sct},\underline{\mathbf{h}}{}_{sct}]}^T={[\underline{\mathbf{e}},\underline{\mathbf{h}}]}^T-{[\underline{\mathbf{e}}{}_{bck},\underline{\mathbf{h}}{}_{bck}]}^T$ into (7) , we obtain the domain (or object) equation:

$$\left[\begin{array}{c}\underline{\mathbf{e}}{}_{bck}\\ \underline{\mathbf{h}}{}_{bck}\end{array}\right]=\left[\begin{array}{c}\underline{\mathbf{e}}\\ \underline{\mathbf{h}}\end{array}\right]-\left\{{\mathcal{H}}_b^{-1}\left(\left[\begin{array}{cc}0& -\mathrm{j}{\mu }_b\omega {\chi }_{\mu }\\ \mathrm{j}\omega {ϵ}_b{\chi }_{ϵ}& 0\end{array}\right]\left[\begin{array}{c}\underline{\mathbf{e}}\\ \underline{\mathbf{h}}\end{array}\right]\right)\right\},$$

(9)

that can be sampled with an operator, say it ${\mathcal{M}}_D$, that extracts fields only in the imaging domain D at the discretization grid points $z_D$ [16], i.e., :

$$\left[\begin{array}{c}\underline{\mathbf{e}}{}_{bck}(z_D,{\omega }_n)\\ \underline{\mathbf{h}}{}_{bck}(z_D,{\omega }_n)\end{array}\right]-\left[\begin{array}{c}\underline{\mathbf{e}}(z_D,{\omega }_n)\\ \underline{\mathbf{h}}(z_D,{\omega }_n)\end{array}\right]={\mathcal{M}}_D\left\{{\mathcal{H}}_b^{-1}\left(\left[\begin{array}{cc}0& -\mathrm{j}\omega {\mu }_b{\chi }_{\mu }\\ \mathrm{j}\omega {ϵ}_b{\chi }_{ϵ}& 0\end{array}\right]\left[\begin{array}{c}\underline{\mathbf{e}}\\ \underline{\mathbf{h}}\end{array}\right]\right)\right\}.$$

(10)

Equations (8)–(10) are the basic equations to solve the inverse scattering problem which is non-linear since both the contrast functions ${\chi }_{ϵ}$ and ${\chi }_{\mu }$ and the field $\underline{\mathbf{e}}$ and $\underline{\mathbf{h}}$ are unknowns of the problem [17]. Moreover, the problem is ill-posed, too, and effective regularization strategies must be adopted for stable solutions against the presence of noise on data [18].

3. Microwave Imaging Profilometry

3.1. Linearized Approach for Plasma Slabs

Let us consider a one-dimensional non-homogeneous slab extending in the range $[z_1,z_2]$ embedded in a non-homogeneous background, see Figure 1a. Normal incidence plane waves are used to probe the slab at difference frequencies with phase reference equal to zero at $z=z^{\left(\mathrm{View}\right)}(z\ge z_1)$. Accordingly, the scattered field can be gathered at a given abscissa $z_1\le z_S=z^{\left(\mathrm{Meas}\right)}\le z_{min}$ (reflection measurement) or $z_{max}\le z_S=z^{\left(\mathrm{Meas}\right)}\le z_2$ (transmission measurement).

To retrieve nonmagnetic (${\mu }_r\left(\underline{\mathbf{r}}{}^{\prime }\right)={\mu }_b$) dielectric profile of inhomogeneous isotropic plasma slabs, we consider a first order linearized approach under the Rytov approximation, instead of solving the exact optimization problems stated by Equations (8)–(10). Linearized approaches, such as the Born (BA) and Rytov (RA) approximations, can be used to recover small and extended (with respect to the wavelength) weak slabs, respectively, i.e., slabs showing low contrast ($\left|\right|{\chi }_{ϵ}\left|\right|<<1$). In particular, the Rytov approximation is useful for scatterers far extended in terms of the probing wavelengths (anyway, RA reduces to BA when the scatterers become smaller and smaller with respect to the probing wavelengths [19]). In this respect, to fulfill model hypothesis, it is worth noticing that assuming a relatively high-frequency probing band (with respect to the cyclotron frequency) entails several positive fallouts. First, at probing frequencies much higher than the cyclotron frequency, the plasma can be assumed as an isotropic medium since the permittivity tensor becomes diagonal with the terms all equals each other [2,20]. Second, at probing frequencies such that ${\omega }^2>>{\nu }^2$ the plasma contrast becomes:

$${\chi }_{ϵ}=-\frac{{\omega }_p^2}{{\omega }^2+{\nu }^2}-\mathrm{j}\frac{{\omega }_p^2\nu }{\omega ({\omega }^2+{\nu }^2)}\approx -\frac{{\omega }_p^2}{{\omega }^2}-\mathrm{j}\frac{{\omega }_p^2\nu }{{\omega }^3},$$

(11)

wherein:

$${\omega }_p=\sqrt{\frac{n_e{e}^2}{m_e{ϵ}_0}},$$

(12)

$n_e$, e and $m_e$ being the electron density, electron charge, and electron mass, respectively. Under the above approximations, considering the frequency difference (FD) formulation described in Section 2, the unknown total electric field can be approximated by the incident field according to the RA, to obtain the following linear equation:

$$\underline{\mathbf{f}}{}_{sct}={\mathcal{M}}_S\left\{{\tilde{\mathcal{R}}}_b^{-1}{\chi }_{ϵ}\underline{\mathbf{e}}{}_{bck}\right\},$$

(13)

wherein $\underline{\mathbf{f}}{}_{sct}=\underline{\mathbf{e}}{}_{sct}/\underline{\mathbf{e}}{}_{bck}$ and ${\tilde{\mathcal{R}}}_b^{-1}$, replacing ${\tilde{\mathcal{H}}}_b^{-1}$ unless unessential constant, takes into account the normalization by the incident field values at the measurement point $z^{\left(\mathrm{Meas}\right)}$ introduced by the Rytov approximation [19] and the frequency dependence of the (multifrequency) scattering operator. As a result, in order to solve for the plasma constitutive parameters ${\omega }_p^2$ and $\nu$, the frequency dependent complex valued problem (13), can be recast in terms of frequency independent real valued problem as already exploited in [7], i.e.,:

$$\left[\begin{array}{c}Re\left(\underline{\mathbf{f}}{}_{sct}\right)\\ Im\left(\underline{\mathbf{f}}{}_{sct}\right)\end{array}\right]=\left[\begin{array}{cc}-\frac{1}{{\omega }^2}{\mathcal{M}}_{S}Re({\tilde{\mathcal{R}}}_b^{-1}\underline{\mathbf{e}}{}_{bck})& \frac{1}{{\omega }^3}{\mathcal{M}}_{S}Im({\tilde{\mathcal{R}}}_b^{-1}\underline{\mathbf{e}}{}_{bck})\\ -\frac{1}{{\omega }^3}{\mathcal{M}}_{S}Im({\tilde{\mathcal{R}}}_b^{-1}\underline{\mathbf{e}}{}_{bck})& -\frac{1}{{\omega }^2}{\mathcal{M}}_{S}Re({\tilde{\mathcal{R}}}_b^{-1}\underline{\mathbf{e}}{}_{bck})\end{array}\right]\left[\begin{array}{c}{\omega }_p^2\\ {\omega }_p^2\nu \end{array}\right].$$

(14)

Problem (14) can be solved as a whole or by splitting it into two subsequent subproblems wherein the first step copes with the retrieval of the plasma frequency (i.e., the real part of the plasma contrast) and the second one dealing with the retrieval of the collision rate which is related to the imaginary part of the plasma contrast. This approach has been proposed in [7] wherein, considering the very low effect of the plasma losses on scattered field, the overall problem can be “decoupled”. This notwithstanding, in the following, we concentrate our attention only on the retrieval of the plasma angular frequency. As a matter of fact, experimental studies state that the collision rate is almost one order of magnitude smaller that the cyclotron frequency, which in turn, is much lower that the probing frequency. This entails to cope with a lossless plasma and the whole entire problem (14) will be reduced to the following one:

$$\left[\begin{array}{c}Re\left(\underline{\mathbf{f}}{}_{sct}\right)\\ Im\left(\underline{\mathbf{f}}{}_{sct}\right)\end{array}\right]=\left[\begin{array}{c}-{\omega }_{max}^2\frac{\omega }{{\omega }^2}{\mathcal{M}}_{S}Re\left({\tilde{\mathcal{R}}}_b^{-1}\underline{\mathbf{e}}{}_{bck}\right)\\ -{\omega }_{max}^2\frac{\omega }{{\omega }^3}{\mathcal{M}}_{S}Im\left({\tilde{\mathcal{R}}}_b^{-1}\underline{\mathbf{e}}{}_{bck}\right)\end{array}\right]\left[\begin{array}{c}\frac{{\omega }_p^2}{{\omega }_{max}^2}\end{array}\right],$$

(15)

with ${\omega }_{max}^2$ the maximum probing angular frequency used as normalization constant introduced for numerical stability of the inversion process [7].

3.2. Sparsity-Promoting Recovery Approaches

Equation (15) is a linear but still ill-posed one which must be solved in a regularized fashion. However, as stated in [10] and investigated in the case of plasma [7], minimum energy solutions approaches, such as Tikonov regularization [21,22] are completely unable to recover neither the shape nor the constitutive parameters of arbitrarily shaped plasma slabs. For this reason, CS-based approaches are worth to be considered, since they are in principle suitable to overcome the intrinsic limitations of the minimum energy solution related to the kind and number of actual parameters which can be conveyed back from multifrequency single view inverse scattering problem.

Let us consider a generic linear problem of the following kind:

$$\mathsf{\Phi }x=f,$$

(16)

x is the $N\times 1$ unknown vector and $\mathsf{\Phi }$ is the $M\times N$ matrix, the so-called sensing matrix, and f the $M\times 1$ data vector of the multifrequency scattered field ordered according to (15). Let us suppose that a convenient representation matrix $\mathsf{\Psi }$ exists such that the unknown can be expanded as $x=\mathsf{\Psi }s$ with only a small number of the coefficients s different from zero. According to CS theory, it is possible to solve the inverse problem (16) even if $M<<N$ but is sufficiently larger than the number S of coefficients different from zero (with $S<M<N$). However, it is worth underlining that the number of the measurements M should be anyway in the order of the degrees of freedom of the scattered field [7], which in the case of the slab can be calculated as [10]:

$$I=2L_{slab}\sqrt{{ϵ}_b}\left(\frac{{\lambda }_M}{c_0}-\frac{{\lambda }_m}{c_0}\right),$$

(17)

$L_{slab}$ being the extent of the inhomogeneous slab, and ${\lambda }_m$ and ${\lambda }_M$ the minimum and maximum wavelength used to probe the scenario. Accordingly, it is possible to solve the inverse linear problem $\mathsf{\Phi }\mathsf{\Psi }s=\mathcal{P}s=f$ by means of different optimization constrained problems [23,24] based on the minimization of the ${\ell }_1$-norm of the unknown when represented in a proper sparse basis. The commonly adopted one in imaging problems, is the rectangular pulse basis function, commonly known as pixel and voxel in 2D and 3D, respectively. However, adopting such a basis, sparse recovery approaches are successful only if the scatterer resembles a point-like scatterer, whereas, in the case at hand, the plasma embeds almost the whole chamber which, in turn, is several probing wavelengths extended.

On the other hand, if a step-wise constant representation of the unknown profile is adopted, which is a reasonable assumption for the problem at hand, ${\ell }_1$ optimization can be still exploited. One possibility is offered CS-based approaches involving as objective function the first derivative of the unknown [7]. For these reasons, we could consider the following recovery problems:

$$\begin{array}{ccc}\hfill \mathrm{argmin}{\parallel \frac{d\left({\tilde{\omega }}_p^2\right)}{dz}\parallel }_{{\ell }_1}& \hfill \mathrm{subject}\mathrm{to}& \hfill {\parallel \mathcal{P}{\tilde{\omega }}_p^2-f\parallel }_2\le \delta ,{\tilde{\omega }}_p^2\ge 0,\end{array}$$

(18)

wherein ${\tilde{\omega }}_p^2=\frac{{\omega }_p^2}{{\omega }_{max}^2}$. In (18), ${\parallel ·\parallel }_{{\ell }_1}$ and ${\parallel ·\parallel }_2$ denote the ${\ell }_1$- and $L^2$-norms, respectively, and the parameter $\delta$, which depends on the measurement error (noise on data) and model error (field approximation), is a positive defined parameter. In (18) the minimization of the ${\ell }_1$-norm promotes the search of solutions with a sparse derivative among all the solutions consistent with the measured data (within the given error threshold $\delta$). In addition, taking into account the real positive value of the unknowns, a further constraint (${\tilde{\omega }}_p^2\ge 0$) can be enforced in (18) in order to reduce the search space by convex optimization and improve the reconstruction results. The approach (18) is known as basis pursuit denoising (BPDN) or least absolute shrinkage and selection operator (LASSO) problem [24].

In [7], the sparsity-promoting approaches (18) have been proved to be capable in retrieving synthetic step-wise constant axial plasma profile which can be represented with few coefficients in terms of their first derivative. However, as the 1D single view multifrequency inverse scattering problem is severe limited, by its very nature, in terms of number of degrees of freedom (I), the number of non-null coefficient of any suitable representation basis should stay very few in order to achieve reliable reconstructions. This means that when smoother ad smoother profiles are considered, the adopted sparsity-promoting approach (18) is expected to be unsatisfactory.

For this reason, a “relaxed” version of (18) can be considered, i.e.,:

$$\begin{array}{ccc}\hfill \mathrm{argmin}{\parallel \mathcal{P}{\tilde{\omega }}_p^2-f\parallel }_2& \hfill \mathrm{subject}\mathrm{to}& \hfill {\parallel \frac{d\left({\tilde{\omega }}_p^2\right)}{dz}\parallel }_{{\ell }_1}\le \gamma ,{\tilde{\omega }}_p^2\ge 0,\end{array}$$

(19)

wherein the role of the objective function with that of the constraints are exchanged. Approach (19) is expected to work better than (18) with profiles which are not exactly step-wise constant. Also, in this case, the assessment of the parameter $\gamma$ is not a simple task. In this respect, prior information can be exploited to assess the optimal choice of the regularization parameter that mainly depends on the sparsity of the unknown at hand. However, the optimal choice of this threshold cannot be a priori established, pointing out one could a priori establish the range of the threshold and compare the results obtained for different values of $\gamma$.

4. Numerical Assessment Towards Benchmark Examples

The two benchmarks we refer to as “flat-top” and “hollow-core”, respectively, represent expected electron density axial distribution in ECRIS-like devices [5]. The synthetic data have been generated considering the convolution between a rectangular window and a Hann window large 10 cm and 0.7 cm, respectively with a discretization grid of 256 cells. As a result, smoother profiles than those tackled in [22] can be considered. The total and the background fields have been simulated in the frequency range [17.5–34.5] GHz considering 87 evenly spaced frequency values in the solution of the forward problem stated by (3), thus the scattered field data has been obtained as difference between them. The solution of the forward problem can be performed by means of a standard conjugate gradient scheme as that proposed in [16]. The choice of the adopted frequency band is related to several requirements. The first one is the need to consider magnetically confined plasma as isotropic and weak scattering medium, and such assumption is fulfilled when the probing frequency is much higher than the cyclotron frequency, hereafter considered at $2.45$ GHz. The second one, is related to the need for working with a large bandwidth in order to achieve the number of degrees of freedom I as large as possible, see Equation (17). Last, but not least, at this frequency band, also in centimeter-sized compact reactors, the propagation can be reasonably approximated as one-dimensional one along the chamber’s axis. According to the above reasonings, we have set the frequency band of the fundamental mode in a rectangular waveguide WR34 that can be used to feed a horn antennas opening in the front and back walls of the plasma reactor [3].

We have considered three kind of measurement/background configurations which resemble possible experimental setup:

• A single transmitting and receiving antenna measuring the reflection coefficient in a free-space homogeneous background (reflection-only measurement);
• Single transmitting and receiving antenna measuring the reflection coefficient in presence of a PEC surface;
• A transmitting and receiving antenna measuring the reflection coefficient and a receiving antenna measuring the transmission coefficient in a homogeneous free-space background (reflection and transmission measurement).

The three cases above can be seen a possible operational scenario where the reflection from metallic environments introduced by the plasma reactor can be neglected or not. This is the case, for example, when large and small reactors are in order, respectively.

To simulate the above measurement/background configurations, an absorbing, or free-space condition, can be very easily introduced, in the case at hand, by considering losses in PML layers ($z_1-d\le z\le z_1$ and $z_2\le z\le z_2+d$) where the medium is gradually modified introducing losses according appropriate permittivity and permeability ${ϵ}_r[1-\mathrm{j}\sigma \left(z\right)]$, ${\mu }_r[1-\mathrm{j}\sigma \left(z\right)]$ to have ideally no reflections from the PML [11,12,14].

For the sake of simplicity, we have chosen a simple polynomial grading for the PML layers conductivity profiles: $\sigma \left(z\right)={\sigma }_{max}{\left(\frac{z_1-z}{d}\right)}^m$ for $z_1-d\le z\le z_1$ and $\sigma \left(z\right)={\sigma }_{max}{\left(\frac{z-z_2}{d}\right)}^m$ for $z_2\le z\le z_2+d$; but more complex profiles can be used [25]. On the other hand, when PEC layer is considered, the PML is completely removed and perfect reflecting surface boundary conditions are enforced at $z=z_2$.

It is worth pointing out that for the solution of the forward problem, the plasma profiles are assumed to be lossy with the collision rate profile one order of magnitude smaller than the plasma angular frequency. As a result, this introduces a model error in the scattered field data which avoids the so-called “inverse crime”. Furthermore, the scattered field data have been corrupted by random gaussian white noise of $5\%$. Finally, the number of the cells in the computational grid are set to $N=256$.

Problem (19) has been solved by means of the environment CVXPY 1.0 for the solution of convex optimization problems [26]. In doing this, the choice of the threshold value $\gamma$ has been set equal to the exact ${\ell }_1$-norm of the actual unknown. Although this exact value cannot be chosen in any actual experiment, it allows an understanding of how the value of $\gamma$ can be set on the base of expected a priori information on the degree of sparsity of the profile. Furthermore, for our investigation, such a choice allows the making of a fair comparison among the achieved results in the three different cases. For the two benchmarks profiles at hand it is $\gamma =6.22$ (flat-top) and $\gamma =7.14$ (hollow-core)

Finally, the required time to solve the inverse scattering problems takes very negligible time (below few seconds on a standard PC) due to the one-dimensional nature of the problem and then the size of the relevant matrix operator.

4.1. Reflection-Only Measurement with Single Antenna

We have considered free-space background with PML layers with thickness of 2 cm, $\sigma =6$, $m=4$, for a sketch the reader can refer to Figure 1a without considering the rightmost antenna. On the other hand, we have considered a non-homogeneous background with a PEC boundary condition on the right hand side of the computational domain ($z_2=25$ cm) to simulate the plasma chamber with only one antenna, see Figure 1b.

The imaging domain, $[z_{min},z_{max}]$, is long 20 cm with the source point $z^{\left(\mathrm{View}\right)}$ and the measurement point $z^{\left(\mathrm{Meas}\right)}$ placed at $2.5$ cm far from the imaging domain both in the case of single antenna and two antennas.

The imaging results are reported in Figure 2 and Figure 3 for both the homogeneous and non-homogeneous background, respectively. As it can be see, in the case of PEC layer, the reconstructions brings better result with a reconstruction error (defined as $err=\frac{{∥x_{true}-x_{rec}∥}_2^2}{{∥x_{true}∥}_2^2}$ ) of $9.39\%$ for the flat-top and $13.03\%$ for the hollow core, with the plasma profile not heavily underestimated as in the case of the homogeneous background, wherein the reconstruction errors achieve $47.85\%$ for the flat-top and $27.02\%$ for the hollow core. This means that the presence of a reflecting surface allows to the achievement of better results in terms of accuracy of the reconstructed profiles. On the other hand, it is worth underlining that while the reconstruction of the flat-top profile is almost satisfactory, Figure 3a, in the case of the hollow-core profile, neither the presence of the PEC allows imaging of the main feature, although differences about the two different recovered profiles are evident, see Figure 3a,b. This may be due to a lower degree of sparsity of the hollow-core profile than that of the flat-top one and to the requirement by (19) to find solutions (profiles) with bounded ${\ell }_1$ norm compatible with the objective function.

4.2. Reflection-Transmission Measurement with Two Antennas

In this case, two antennas are considered, see Figure 1a, with one (active) antenna radiating the probing field and measuring the reflection coefficient, and a second (passive) antenna measuring only the transmitting coefficient. Opposite to the case with a single antenna, in this case the non-homogeneous background accounting for two PEC surfaces has not been considered. First, from a modeling point of view, the use of PEC surfaces entails an electromagnetic cavity which enforces solutions only at the resonance frequencies given by $f_n=\frac{n{c}_0}{2L_{\mathrm{cavity}}}$, i.e., when the cavity length ($L_{\mathrm{cavity}}$) is an integer multiple of the probing half-wavelength, for a given finite set of integer $n^{*}$ indexes depending on the adopted frequency band. Secondly, when the plasma partially embeds the cavity, the resonance frequencies shift in frequency, so that the expected large error on the background field approximation prevent using a linear model such as the Rytov approximation. Finally, from an applicative point of view, since in general, the front and back walls of the plasma reactors host two aperture antennas [3], just a (possible small) fraction of the transmitted microwave power will be reflected, while the most part will be gathered by the antennas. For these reasons, we have considered only homogeneous background.

The imaging results are reported in Figure 4 and the reconstruction errors achieve $err=1.01\%$ and $err=3.56\%$ for the flat-top and hollow core, respectively. Accordingly, it can be shown as the presence of a second passive antenna achieve better results in terms of accuracy of the retrieved profiles when compared to the single antenna in free-space homogeneous background. On the other hand, while the reconstruction of the flat-top profile is nearly optimal, also in this case, the reconstruction of the hollow-core profile is not satisfactory to apprise the main feature of the profile, although some small differences between the two retrieved profile are evident (as in the case of single antenna and PEC surface) and the same possible explanations hold true.

5. Discussion of the Results through Singular Value Decomposition Analysis

To understand the different results achieved under the two different considered measurement configurations and related backgrounds, we have considered the SVD analysis of the relevant scattering operator $\mathcal{P}$ for the different analyzed cases.

First of all, it is worth underlining that even if we are not considering minimum energy solution approaches ($L^2$-norm), the SVD analysis can be retained still valid to understand the main features among the different measurement configurations and backgrounds, since we are considering sparse promoting approaches rather than exact compressive sensing-based approaches. Indeed, (19) entails a standard minimum energy solution objective function equipped with a regularization (penalty) ${\ell }_1$-norm enforcing step-wise sparsity. For this reason, we numerically perform the SVD of the relevant scattering operator with the aim to give a qualitative interpretation of the different results shown in the previous Section.

Let us denote with $\{u_n$, ${\sigma }_n$, $v_n\}$ the singular value decomposition of the matrix operator $\mathcal{P}$ mapping the data-to-unknown relationship in (19), wherein $u_n$ represents the left singular vectors (spanning the object space), $v_n$ the right singular vectors (spanning the data space) and ${\sigma }_n$ the singular values, in such way $\mathcal{P}{u}_n={\sigma }_n{v}_n$ and ${\mathcal{P}}^{+}{v}_n={\sigma }_n{u}_n$, wherein ${}^{+}$ stand for the adjoint of the operator [18].

For all the measurement configurations adopted in the numerical analysis, we have analyzed the following metrics based on the SVD analysis, i.e.,:

• behavior of the singular values —the logarithmic plot of the singular values as ordered in non-increasing fashion. Indeed, as the scattering operator is a compact one [18], its singular values exhibits an exponential decay after a given threshold index I (analytically expressed by (17)), which indicates the maximum number of the degrees of freedom, and hence of the parameters which can be conveyed back by the recovery procedure;
• spectral coverage (SC) defined as:

  $$SC=\sum _{n=0}^{N_T}{\left|{\tilde{U}}_n\right|}^2$$

  (20)

  wherein $\tilde{U}$ is the Fourier Transform of the left singular vectors and $N_T$ the truncation index used as regularization parameter in TSVD approach [10]. It is a measure of the class of profiles which can be actually retrieved in the object domain by the inversion process [10].
• point spread function (PSF) defined as [27]:

  $$PSF(z-z_0)=\sum _{n=0}^{N_T}{u}_n\left(z\right)u_n{\left(z_0\right)}^{*}$$

  (21)

  wherein ${\left(\right)}^{*}$ stands for conjugation, and $z_0$ the abscissa with respect the PSF is considered. For the case at hand, we set $z_0=12.5$ cm, which is the center of the imaging domain. The point spread function is a direct measure of the ultimate attainable spatial resolution [27].

The SVD analysis is reported in Figure 5. It can be seen that the singular values are different for the considered configurations. Indeed, for the PEC case the lower order singular values show a greater magnitude while vanishing to zero (like in the case of the free-space background configuration). On the other hand, in the case of two antennas in homogeneous background, the magnitude of the singular values keeps almost always greater than the single antenna (without and with PEC).

The SVD analysis give us the fundamental answer to understand the differences in the reconstruction process for the three adopted measurement configurations. Indeed, two main comments are in order. First, the presence of a perfect reflecting plane (PEC) introduces an increase in the singular values magnitude as compared to the homogeneous case. This is actually shown by the blue and red dashed line in Figure 5a. Interestingly, such an increase entails the advantage to deal with a more stable inverse scattering problem and with an increase of the information content [9]. Second, when we move to the case of two antennas, it can be easily understood as this configuration is almost equivalent to the case of single antenna operating in presence of a PEC surface. Indeed, this can be explained considering that the PEC surface acts as a “secondary” antenna which radiates (reflects) the field in such a way a virtual transmission coefficient can be measured by the actual antenna, as superimposed to the reflection coefficients that it would be measured without the PEC surface. As a result, when two antennas are considered, the reconstruction result achieves better results as in presence of a PEC surface. This analysis shows that the use of a second passive antenna is almost equivalent to have a reflecting surface behind the plasma slab. Anyhow, in all those cases wherein free-space approximation can be assumed (f.i., large meter-sized machine), two antennas can be used to achieve an enhancement in the reconstruction process. This has an immediate fallout in the design of the measurement setup depending on the constraints enforced by the experimental facilities.

On the other hand, the SVD analysis shows also that SC and PSF attained in three measurement configurations are identical, unless some very negligible differences due to numeric evaluation of the singular value system. These are confirmed by the same reconstruction performance achieved by the single antenna with PEC surface and two antennas in homogeneous background, for which the harmonic content is almost exactly the same. As a result, the improvement in the reconstruction in presence of a PEC surface can be attributed only to the difference in the singular values spectrum.

6. Conclusions

In this paper, we have addressed microwave imaging profilometry through the solution of an inverse scattering problem with a finite-difference formulation which takes into account, under the adopted assumptions, the effects of the metallic case represented by the plasma reactor. It has been shown that the presence of a metallic surface profoundly benefits the inversion recovery process when faced with sparsity-promoting approaches. In addition, the developed analysis shows the equivalence among two possible measurement configurations which are a single transmitting-receiving antenna in presence of a PEC surface and two antennas (one active and the second passive) able to measure also the transmission coefficients. Although the imaging capabilities are enhanced in presence of totally reflecting surface and with two antennas, microwave imaging profilometry still shows severe limitations in terms of accuracy in the recovery of arbitrary-shaped electron density profiles, i.e., not exactly step-wise constant profiles. Further activities are currently being devoted to exploit more efficient unknown representations and possible a priori information in order to achieve better reconstructions accuracy and to deal with a class of profiles wider than that considered in this paper.