Nonlinear Problems of Equilibrium Charge State Transport in Hot Plasmas

Abstract

The general coupling between particle transport and ionization-recombination processes in hot plasma is considered on the key concept of equilibrium charge state (CS) transport. A theoretical interpretation of particle and CS transport is gained in terms of a two-dimensional (2D) Markovian stochastic (random) processes, a discrete 2D Fokker-Plank-Kolmogorov equation (in charge and space variables) and generalized 2D coronal equilibrium between atomic processes and particle transport. The basic tool for analysis of CS equilibrium and transport is the equilibrium cell (EC) (two states on charge and two on space), which presents simultaneously a unit phase volume, the characteristic scales (in space and time) of local equilibrium, and a comprehensive solution for the simplest nonlinear relations between transport and atomic processes. The space-time relationships between the equilibrium constant, transport rates, density distributions, and impurity confinement time are found. The subsequent direct calculation of the total and partial density profiles and the transport coefficients of argon impurity showed a strong dependence of the 2D CS equilibrium and transport on the atomic structure of ions. A model for recovering the recombination rate profiles of carbon impurity was developed basing on the CS equilibrium conditions, the derived relationships, the data about density profiles, plasma parameters and ionization rates.

1. Introduction

Understanding the behavior and control of impurities as an inevitable component of magnetically confined plasmas continue to be a topic of considerable importance in fusion research programs. Steady state mode of operation with impurity equilibrium is of vital to the success in modern experimental studies of improved confinement and stability of standard H-modes on tokamaks and stellarators [1,2].

In these devices, the radial transport of plasma particles across the magnetic field is accompanied by random change in their charge states (CS’s) due to ionization-recombination processes. The probabilistic imposition of these atomic processes on the transport of impurity particles occurs in an unusual phase space, which includes the discrete space of (Z + 1) CS’s. So, the resulting transport turns out to be at least two-dimensional (2D), while its structure involving an intersection of statistically independent processes is unexpectedly complicated. A consistent approach to this general 2D transport problem could suggest the use of multidimensional (mD) Fokker-Planck-Kolmogorov equation (FPKE) [3,4,5].

Meanwhile, for more than four decades, the description of impurity transport remained implicitly simplified and was reduced to a 1D diffusive-convective model (DCM), represented in 1D impurity transport codes (see, for example, the STRAHL code [6]). Here, the standard 1D continuity equation was supplemented in the right-hand side with also a standard 1D flux balance of ionization-recombination processes interpreted as random sources and sinks of charged particles. So, the required 2D description was just replaced by the sum of these 1D parts. The corresponding misinterpretation fell into an inevitable contradiction with the systematic approach that the theory of the generalized FPKE could provide. Indeed, the essential cross terms (mixed derivatives) representing the coupling of atomic processes with particle transport were also lost in the DCM equations. Moreover, the required discretization of these equations revealed and returned this 2D transport problem but only in the form of a 2D grid model (GM) in charge and radial variables (see discussion in [6] on page 6).

In the extensive practice of modeling medium and heavy impurities and with an increase in the rates of ionization-recombination it was found that the DCM sensitivity to the empirical transport remains extremely low [7,8,9,10,11,12,13] and just unacceptable for W impurity [14,15].

Nevertheless, the leading role in interpreting the behavior of impurities was assigned to various mechanisms of particle transport in the 6D phase space (of coordinates and velocities) in terms of the Boltzmann kinetic equation and its consequences carefully studied in the neoclassical transport theory [16]. But in the basic equations of neoclassical transport, the charge state as a mandatory variable was ignored along with impurity charge state space (see, e.g., Equation (2.2) in [16]).

A consistent critique of the standard 1D approach using DCM was recently given along with a completely new understanding of impurity transport in terms of the 2D CS transport [17]. The main unresolved physical problem outlined here was that 1D (radial) transport of particles would have to be clearly distinguished from the generalized 2D transport in the proper 2D phase space, where the particle transport is only part of CS transport. It is important, therefore, to make, first of all, a terminological distinction between particle and CS transports, that is, to introduce and develop the concept of CS transport. But the initial idea was to consider the role of ionization/recombination processes not as sources and sinks in the standard 1D transport equations, but on equal terms with particle transport, as a kind of transport phenomenon [18,19,20]. Due to mentioned imposition, the basic components of impurity transport turn out to be mixed according to probabilistic laws in a resulting generalized 2D transport. This understanding has been reformulated step by step in framework of traditional analysis via diffusive D and convective V coefficients in terms of CS transport [21,22] and finally allowed us to clarify the structure and probabilistic nature of impurity transport [17].

However, before any subsequent transition to equations describing the details of the distribution function with operators, the basic (first) principles for any subsequent analysis of the transport of charged particles are mandatory. Indeed, although, there is no particle without CS, as well as no CS of no particle, but there is a fundamental difference between the transport of particles and the transport of their CS’s due to ionization and recombination processes. In fact, let the probability of CS transport by particles be H(P), and by atomic processes be H(A), then the probability of events P + A is expressed, as is known, by the general Formula (1):

$$H(P+A)=H(P)+H(A)-H(P)\cdot H(A)$$

(1)

which shows that the generalized impurity transport includes three mandatory components: (i) CS transport by particles or p-transport, (ii) CS transport by atomic processes or a-transport, and, finally, (iii) CS transport represented by probabilistic intersection of a-transport and p-transport but irreducible to any of them. Note that the latter is the most represented case of events related to impurity transport in plasma. It is clear from Formula (1) that DCM is limited to analyzing only the first two components, implicitly (and erroneously) assuming that events P and A are incompatible.

Following the logic of Formula (1), the events caused by the intersection of $H(P)\cdot H(A)$ and ensuring the local occurrence of any ion, are initially and in principle non-local, mixed by products of their probabilities in subsequent expressions, and as a result are indistinguishable due to the probabilistic nature and logic of the sum of P + A.

The significant advantage of the CS concept of impurity transport is the self-consistency of the developed 2D approach. The 2D Markovian impurity equilibrium directly provides the possibility to calculate D and V using the 2D FPKE adopted to the discrete CS phase space [17]. This makes it possible to perform transport studies without any additional empirical coefficients, while all components of 2D transport are expressed in more general way. The additional possibilities of the local (core and edge) transport studies also follow from the proposed 2D modelling.

From a theoretical point of view, we are faced here with a complex 2D structure of the impurity equilibrium and with its inevitable symmetry along with the possibility of the self-consistent determination of D and V. Moreover, the standard 1D descriptions of impurity ionization equilibrium or the stationary distribution (radial profiles) of the total impurity density is just included in this 2D model of Markovian equilibrium. We assume a 2D Markov equilibrium of impurities, since it follows from the so-called ergodic principle [3]. In this case, the CS transport of charged particles should be considered as an equilibrium transport. The transport rates are locally determined by a 2D probabilistic balance between the motion of impurity particles and the change in their charges due to atomic processes. In addition, the particle conservation law implied in the DCM equations is just generalized to the conservation of impurity CS moving in the 2D phase space. It is that results in a symmetric representation of transport in charge and radial variables on equal terms as suggested by the generalized FPKE.

Thus, the concept of impurity CS transport allows one to draw an unexpected conclusion about the determining role of atomic processes in impurity transport [17].

Clear evidence of this interpretation can be revealed in experiments and related modelling practice. First of all, this is a coronal equilibrium (CE) of heavy impurities indeed found in the early studies of Cr, Ni, and Mo impurities on the TFR tokamak [23,24], and then confirmed by modelling practice for W on the ASDEX Upgrade tokamak [14,15]. This meant that the equilibrium assumed by the DCM between 1D particle transport and 1D balance of ionization-recombination processes could not be observed, since in fact we are dealing here with the 2D equilibrium of the plasma impurity.

We further note the observations of a noticeable drop in the diffusion coefficient D at the central impurity accumulation and a sharp difference in transport in the plasma core from transport in the rest periphery (see, for example, [25]). These results can be explained by the largest size of the central equilibrium cell (CEC), as clearly seen from the CS transport modelling [17].

The next important experimental evidence is the known coupling found between central accumulation of impurities and diffusion coefficient D: the greater the accumulation of the impurity in the plasma core, the lower D [17,26], and vice versa, that is, the “flat” distribution is systematically associated with an increase D to anomalous values (in the conventional comparison with neoclassical predictions) from the core to the plasma edge.

Finally, there are also two characteristic maxima on the total density profiles. The observed impurity accumulation takes place simultaneously at the edge and in the plasma core, as found, e.g., for W in JET [1] that consistently reproduced in the CS transport model [17].

Meanwhile, the 2D imposition of ionization/recombination on the particle transport results in a fundamental nonlinearity of the 2D transport analysis, which aim is to find and study the particle transport. The setting of the problem cannot be reduced only to the standard 1D linear equations in principle. A more general 2D problem consists in directly modeling the impurity density profiles (both partial and complete) observed in experiments to find the corresponding D and V, using only the atomic database of ionization/recombination rate coefficients, as shown in [17].

However, this 2D analysis gives the nonlinear equations for even a simplest 2D system of four charge states (two in space and two in charge called the equilibrium cell (EC))—leads to a transcendental equation [17]. The general solution was developed basing on a number of reducing 2D schemes and the use of the pseudo-state technique. Thus, the solution of the nonlinear transport problem is represented by two main stages: first, by reducing the general discrete GM of impurity equilibrium to a set of local equilibrium cells and, second, by finding their transport coefficients that provide the calculation of density profiles along with the transport coefficients.

Nevertheless, such direct modelling of the total density profiles suffers from the problems associated with the usual and significant uncertainties of atomic rates. As well known, the most uncertainties appear in the data of cx- and dielectronic recombination rates [8,14,15].

To solve this uncertainty problem, in this paper we propose a new formulation of the transport nonlinear problem along with the developing CS transport model. It is based on the possibility of analyzing the equilibrium transport by separating its absolute and relative scales, which represent the equilibrium constant λ and the corresponding rate ratios [17]. Then, the input data in this model are the typical radial profiles of the total impurity density observed in experiments and ionization rate coefficients. So, we continue here the previous analysis of the CS equilibrium with its absolute scale, λ, which turns out to be common to all emerging nonlinear problems [17].

We further show that an important relationship can be established between λ, the standard transport coefficients and the impurity confinement time ${\tau }_p$. A generalization of the EC concept is developed in the form of the reduced equilibrium cell (REC). The REC concept allows us to proceed directly to impurity equilibrium analysis based on a wide variety of the density profiles as mentioned above.

The structure of the paper is as follows: In the second part, we consider the symmetry features of the basic equations of the CS equilibrium and transport, analyzing from this point of view the same basic formulas obtained earlier in [17]. An extended set of possible schemes for reducing GM to easily analyzed simplifications is presented, also taking into account possible symmetries of the proposed approach. The matrix description of the impurity density profiles is developed and its important relationship with the GM equilibrium is shown through the important scale factor λ. Its dependence on the atomic structure of the impurity ions, the most abundant in plasma, is analyzed in detail in Section 3. The novel CS model and the basic results of recovery of the relative density profile of hydrogen neutrals that determine the impurity charge exchange recombination rate are given in Section 4. Section 5 presents the conclusions.

2. Charge State Equilibrium and Transport

In the classical review of Braginsky [27], an equation representing the impurity transport was given as follows Equation (2):

$$\frac{\partial n_k\left(r,t\right)}{\partial t}+{\nabla }_r{\Gamma }_k\left(r,t\right)=Q_k\left(r,t\right)$$

(2)

where $n_k\left(r,t\right)$ is the density of impurity particles in k-th CS (k = 0, 1,..., Z), ${\Gamma }_k\left(r,t\right)$ is the particle flux of the k-th CS and ${\displaystyle \sum _k{Q}_k\left(r,t\right)}=0$. All changes of $n_k\left(r,t\right)$ due to ionization-recombination processes are described in Equation (2) by the term $Q_k\left(r,t\right)$, but the expression for it was not considered in the review. Nevertheless, it is clearly seen, that the proper phase space for Equation (2) has two independent variables, k and r (and time). Consequently, it would have to be considered in the 2D phase space of these k and r on equal terms using the system approach suggested by FPKE. In this case, the 2D transport represented in Set (2) must also be consistently analyzed in a general series of discrete mD Markovian processes.

2.1. The Discrete Transport and Equilibrium Conditions

In the general mD case, the differential probability distribution function $g\left(x_1,x_2,\dots ,x_m\right)$ with independent variables $x_1,x_2,\dots ,x_m$ can be represented by the FPKE [3,4,5], shown as Equation (3):

$$\frac{\partial g}{\partial t}=-{\displaystyle \sum _i\frac{\partial }{\partial x_i}}\left(A_{i}g\right)+{\displaystyle \sum _{i,j}\frac{{\partial }^2}{\partial x_i\partial x_j}}\left(B_{ij}g\right),$$

(3)

where $A_i$ is a vector of convective fluxes, $B_{ij}$ is a tensor of generalized coefficients of diffusion in the mD phase space, ${\displaystyle \int g{\displaystyle \prod d{x}_i}}=1$. It is also assumed here that the variables x_i and x_j are changed almost continuously. The function $g\left(x_1,x_2,\dots ,x_m\right)$ must satisfy the Smoluchowski Equation with integration over mD phase volume ${\displaystyle \prod d{x}_i}$. The derivation of Equation (3) can be found elsewhere [4,5].

$Q_k\left(r,t\right)$ in Equation (2) can be generally found from comparing Equations (2) and (3). The mixed derivatives included in the resulting expression in the 2D case of interest are related to the volume $\partial r\cdot \partial k$, which is directly given in these derivatives and just not present in the standard DCM. (see, Equation (2.2) in [16]). It should be noted that the mandatory integration of the general (mD) form of the kinetic equation on velocities to derive Equations (2) or (3) cannot, in principle, eliminate the differentials of variables with their product, $\partial r\cdot \partial k$, that determines the unit phase volume. Their spatial scale represents the non-locality of the 2D coupling between a- and p-transports in the EC. Therefore, if the conventional 1D equation does not contain them, then it follows that it contradicts and does not agree with both Formula (1) and the mD FPKE (3).

To represent atomic processes in convenient discrete form, Equations (2) and (3) must obviously be modified. Consider the discrete 2D GM, which is usually used to discretize Equation (2) in transport codes [6]. We introduce a grid probability distribution function, $g_{kn}$, where k = 0, 1, 2, ..., Z and n is the radial index instead of r; n = 1, 2, ..., N, N = m + 1, where m is the number of cells. Since $g_{kn}$ is normalized as follows ${\displaystyle \sum _{k,n}^{}{g}_{kn}}=1$, we get ${\displaystyle \sum _{k=0}^Z{g}_{kn}}=p_n$, that is, the total local (on n) normalized particle density and ${\displaystyle \sum _{n=1}^N{g}_{kn}}=q_k$, that is, the total normalized density of all ions with the same charge k. $p_n$ and $q_k$ are related to $g_{kn}$ by the distribution functions $f_{kn}$ and ${\phi }_{kn}$ according to Equation (4):

$$g_{kn}=p_n{f}_{kn}=q_k{\phi }_{kn}$$

(4)

and normalizations: ${\displaystyle \sum _{k=0}^Z{f}_{kn}}=1$, ${\displaystyle \sum _{n=1}^N{\phi }_{kn}=1}$.

The CS dynamics on the grid are considered as the random jumps to neighboring positions (grid nodes): by k due to a-transport and/or by n due to p-transport. The a-transport rates are $S_{kn}$ and $R_{kn}$, that is, the total rates of ionization and recombination. The p-transport rates (also in s⁻¹) could be denoted as $w_{kn}^p$ if the CS motion by p-transport is directed from n to n + 1 and as $u_{kn+1}^p$ if the motion is directed from n + 1 to n. Note that these introduced transport coefficients are just a modified form of the conventional D and V and present the same, albeit unidirectional, particle fluxes. So, for rough numerical estimates, D and V could be obtained from these 1D matrix coefficients simply as $D\propto \left(w_{kn}^p+u_{kn+1}^p\right)/2$ and $V\propto \left(u_{kn+1}^p-w_{kn}^p\right)$ for the most abundant k states (see in detail in [28]).

The discrete 2D form of Equation (3) assumes the choice of the necessary system of pairwise incompatible events (in accordance with the general probability summation formula) associated with states without matching indices. Therefore, we consider pairs of diagonally conjugate states: (k, n) and (k + 1, n + 1); (k + 1, n) and (k, n + 1) and the corresponding diagonal rates $w_{kn},u_{k+1n+1}$ and ${w^{\prime }}_{k+1n},{u^{\prime }}_{kn+1}$. Then the 2D equation (for example, for the first of the symmetric pairs) is obtained as Equation (5):

$$\frac{\partial g_{kn}}{\partial t}=u_{k+1n+1}{g}_{k+1n+1}-\left(w_{kn}+u_{kn}\right)\cdot g_{kn}+w_{k-1n-1}{g}_{k-1n-1},$$

(5)

where the rate matrix is tridiagonal, Jacobian, singular (it has zero value in the spectrum of its real eigenvalues) and similar to some symmetric real matrix. At any time, t, of the temporal evolution of Equation (5), there is an ergodic limit (see [3]), which is determined by a complete set of these matrices. The stationary solution of Equation (5) is simultaneously the local equilibrium condition, shown as Equation (6):

$$\begin{array}{l}{w}_{kn}{g}_{kn}=u_{k+1n+1}{g}_{k+1n+1},\\ {w^{\prime }}_{k+1n}{g}_{k+1n}={u^{\prime }}_{kn+1}{g}_{kn+1}.\end{array}$$

(6)

The 2D equilibrium is symmetric, assuming the equalities of symmetric pairs of counter CS fluxes. Shown as Equation (7):

$$\begin{array}{l}{w}_{kn}^p{g}_{kn}=u_{kn+1}^p{g}_{kn+1},\\ S_{kn}{g}_{kn}=R_{k+1n}{g}_{k+1n},\end{array}$$

(7)

which can be proved strictly in general for any 2D structures (see in [17]). Therefore, note, that Equation (6) directly follow from these Equation (7) and independent of Equation (5).

The nonlinearity of the problems that arise in the analysis of the CS transport is directly follows from Equations (6) and (7), so we get Equation (8):

$$\frac{g_{k+1n+1}}{g_{kn}}=\frac{w_{k+1n}^p{S}_{kn}}{u_{k+1n+1}^p{R}_{k+1n}}=\frac{w_{kn}^p{S}_{kn+1}}{u_{kn+1}^p{R}_{k+1n+1}}=\frac{c_1{w}_{k+1n}^p{S}_{kn}+c_2{w}_{kn}^p{S}_{kn+1}}{c_1{u}_{k+1n+1}^p{R}_{k+1n}+c_2{u}_{kn+1}^p{R}_{k+1n+1}},$$

(8)

where $c_1$ and $c_2$ are arbitrary constants to be found. By defining Equation (9):

$$\begin{array}{l}{w}_{kn}\equiv c_1{w}_{k+1n}^p{S}_{kn}+c_2{w}_{kn}^p{S}_{kn+1},\\ u_{k+1n+1}\equiv c_1{u}_{k+1n+1}^p{R}_{k+1n}+c_2{u}_{kn+1}^p{R}_{k+1n+1},\end{array}$$

(9)

we obtain the general case of conditions (6) expressed in terms of products of conditional probabilities of particle and charge transport processes. Formally, using the normalization of $g_{kn}$ with Equations (8) and (9), we obtain the set of rate equations. But even simplest expressions for the diagonal CS transport in the single EC give a strong nonlinear relationship between the probability distribution and transport rates (see, e.g., Equation (13) below or Equation (21) in [17]).

Thus, the rates $\left\{w_{kn}^p\right\},\left\{u_{kn}^p\right\}$ and $\left\{S_{kn}\right\}$, $\left\{R_{kn}\right\}$ are related in a nonlinear way to $\left\{w_{kn}\right\},\left\{u_{kn}\right\}$ in accordance with Equations (1) and (5)–(9). Using Equation (4), Equality (7) can be rewritten as Equation (10):

$$\begin{array}{l}{w}_{kn}^p{\phi }_{kn}=u_{kn+1}^p{\phi }_{kn+1},\\ S_{kn}{f}_{kn}=R_{k+1n}{f}_{k+1n}.\end{array}$$

(10)

where both distribution functions are normalized (see after Equation (4)), hence Set (10) can be interpreted as a generalized equilibrium or double coronal equilibrium (DCE). For each cell, there are symmetric pairs of such equalities that define the local DCE conditions. A local unit of phase volume has a pair of necessary transport ratios, shown in Figure 1.

Each impurity particle within a local unit of the discrete phase space is characterized by four states. Consider the relationships: ${\alpha }_{kn}=S_{kn}/R_{k+1n}=g_{k+1n}/g_{kn}$ and ${\beta }_{kn}=w_{kn}^p/u_{kn+1}^p=g_{kn+1}/g_{kn}$. Using the entered relations, the Equations (7) and (10) can be rewritten again in the form of Equation (11):

$$\frac{{\alpha }_{kn}}{{\alpha }_{kn+1}}=\frac{{\beta }_{kn}}{{\beta }_{k+1n}},$$

(11)

where both values in the left part obviously depend only on the rates of ionization and recombination and relate only to a-transport, while in the right part—only to p-transport. Equation (11) defines the impurity equilibrium as the necessary local coupling of a- and p-transports, but within the size of each cell. A transport that satisfies the condition (11) could be called an equilibrium transport.

Whatever the mechanisms of p-transport in the plasma (collisional or turbulent), local equilibrium of impurities can only be provided by the equilibrium p-transport rates. Consequently, the local equilibrium transport is determined by statistically independent rates of ionization and recombination. Since the p-transport corresponds to the a-transport, that is, it turns out to be equilibrium, a generalized DCE of impurity (10–11) occurs.

Thus, local EC is simultaneously: (i) a unit phase space volume and the characteristic scale of the local equilibrium, (ii) self-consistent equilibrium system, which provides a comprehensive analysis of the rate ratios of a- and p-transports, (iii) a crucial tool for analyzing a general equilibrium of impurity CS transport.

Moreover, we argue that there is a fundamental difference in the basic estimates of the impurity equilibrium in the 1D and 2D cases. In contrast to the 1D case, 2D equilibrium relations are established between the fluxes of the same type of processes in each 2D cell (see Equations (7) and (10)) and, in addition, the relations of the relations themselves, as follows from (11). Mixed equilibrium fluxes (see, Equations (8) and (10)), which are the superposition of a-and p-transports, occur only in symmetric diagonal directions with inevitable nonlocal coupling. A similar mixture of atomic processes and particle transport was first derived from the 2D solution of standard equations [20,21].

The radial scale of the nonlocal coupling between particle transport and ionization/recombination is determined by the m value, that is by the local scale of the EC. The largest one is related to the CEC.

Equation (11) couples only the rate ratios, which can be reproduced for impurities with different Z and with significantly different rates. But only they determine the equilibrium transport along with the resulting stationary density distributions. The corresponding similarities of the equilibrium rate ratios can be suspected in experiments, e.g., when there is an apparently similar stationary density profiles of the central accumulation of C [29,30], Ar [31,32] and W [1,33].

Contrary, the absolute values of the atomic rates determine the equilibrium scale, λ, that is, an equilibrium constant of steady state impurity [17]. It turns out to be the same to all cells and reducing GM schemes, as will be shown below. Therefore, the difference in the equilibrium of impurities appears primarily as a function λ(Z), which strongly depends on the atomic structure and the rates of the most abundant ions. For light and mid-Z impurities, these are K- and L-shells of H-, He- and Li-like ions. Normalizing the atomic rates by this λ allows us to develop a generalized elementary 2D model of equilibrium that could be called a reduced equilibrium cell (REC). Formally the REC is just a special case of the EC with λ = 1.

2.2. Reduced Equilibrium Cell

The CS equilibrium and transport in the EC introduced elsewhere [17] is completely determined by the rates of ionization and recombination processes, $S_1$, $S_2$ and $R_1$, $R_2$ respectively. The relations ${\alpha }_1=S_1/R_1$, ${\alpha }_2=S_2/R_2$, $s=S_2/S_1$, $\tilde{r}=R_2/R_1$, ${\beta }_1=w_1/u_1$, ${\beta }_2=w_2/u_2$ are either directly obtained from the given rates of atomic processes, or can be obtained from the equilibrium analysis. In particular, in the case of the EC the Equation (11) is presented in Equation (12):

$$\gamma \equiv \frac{{\alpha }_2}{{\alpha }_1}=\frac{{\beta }_2}{{\beta }_1}=\frac{s}{\tilde{r}}.$$

(12)

Solutions of the equations describing the temporal evolution of densities $p_{1/2}=g_{1/3}+g_{2/4}$, $q_{1/2}=g_{1/2}+g_{4/3}$, $g_{1/3}$, $g_{2/4}$(for short, the subscripts for pairs are changed) include non-zero eigenvalues of the EC matrix equations. As was shown in the EC analysis [17], these turn out to be equal, i.e., $\lambda ={\lambda }^{\prime }={\lambda }_p={\lambda }_q$, where $\lambda =w+u$, ${\lambda }^{\prime }=w^{\prime }+u^{\prime }$ are the constants of the diagonal balance as follows from the FPKE adopted for the EC, ${\lambda }_p=\overline{w}+\overline{u}$, where $\overline{w}=\left(w_1+w_2{\alpha }_1\right)/\left(1+{\alpha }_1\right)$, $\overline{u}=\left(u_1+u_2{\alpha }_2\right)/\left(1+{\alpha }_2\right)$, and ${\lambda }_q=\overline{S}+\overline{R}$, where $\overline{S}=\left(S_1+S_2{\beta }_1\right)/\left(1+{\beta }_1\right)$ and $\overline{R}=\left(R_1+R_2{\beta }_2\right)/\left(1+{\beta }_2\right)$.

Now consider the REC, in which, as noted above, all rates are normalized to λ, that is, ${\tilde{S}}_1=S_1/\lambda$, ${\tilde{S}}_2=S_2/\lambda$, ${\tilde{R}}_1=R_1/\lambda$ and ${\tilde{R}}_2=R_2/\lambda$. We also get that ${\tilde{w}}_1=w_1/\lambda$, ${\tilde{w}}_2=w_2/\lambda$ and ${\tilde{u}}_1=u_1/\lambda$, ${\tilde{u}}_2=u_2/\lambda$. However, ${\alpha }_1$, ${\alpha }_2$, ${\beta }_1$ and ${\beta }_2$ remain the same. In this case, the definition of the REC is a set of equalities $\lambda ={\lambda }^{\prime }={\lambda }_p={\lambda }_q=1$. The REC scheme is shown in Figure 2.

The basic equations derived for the EC are the same for the REC also, as shown in Equation (13):

$$\begin{array}{l}1+{\kappa }^{\prime }xv=\frac{x}{x_e}+\frac{v}{v_e},\\ x+v=\frac{\left(1+{\kappa }^{\prime }\right)}{s}\frac{\left(1+{\alpha }_1\right)\left(1+{\alpha }_2\right)}{\left(1+{\beta }_1\right)\left({\alpha }_1+{\alpha }_2{\beta }_1\right)},\\ \frac{x}{1+{\kappa }^{\prime }x}+\frac{v}{1+v}=\frac{{\alpha }_2}{\left(1+{\alpha }_2{\beta }_1\right)}\cdot \frac{\left(1+{\kappa }^{\prime }\right)}{s}\cdot \frac{\left(1+{\alpha }_1+{\beta }_1+{\alpha }_2{\beta }_1\right)}{\left(1+{\beta }_1\right)\left({\alpha }_1+{\alpha }_2{\beta }_1\right)},\end{array}$$

(13)

where $v\left(x=0\right)\equiv v_e=\left(\kappa -1\right)/\left({\kappa }^{\prime }-\kappa \right)$, $x\left(v=0\right)\equiv x_e=\left(\kappa -1\right)/\left(1-{\kappa }^{\prime }\kappa \right)$ in the notations $x=u_2/R_2={\tilde{u}}_2/{\tilde{R}}_2$, $v=u_1/S_2={\tilde{u}}_1/{\tilde{S}}_2$, ${\kappa }^{\prime }={\beta }_{1}s$, $\kappa ={\alpha }_2\left(1+{\beta }_1/{\alpha }_1\right)/\left(1+{\alpha }_2{\beta }_1\right)$ and ${\tilde{w}}_1={\beta }_1{\tilde{u}}_1$ ($w_1={\beta }_1{u}_1$),${\tilde{w}}_2={\beta }_2{\tilde{u}}_2$ ($w_2={\beta }_2{u}_2$).

For the REC analysis we add to Equation (13) two expressions for ${\alpha }_1\left(p,{\beta }_1\right)$ and ${\alpha }_2\left(p,{\beta }_1\right)$(where $p\equiv p_2/p_1$) taken from the equation ${\lambda }_q=\overline{S}+\overline{R}=1$ and the relation $p={\beta }_1\left(1+{\alpha }_2\right)/\left(1+{\alpha }_1\right)$ in the form of Equation (14):

$$\begin{array}{l}{\alpha }_1=\frac{{\left(1+{\beta }_1\right)}^2}{\left(1+p\right)\left(1+{\beta }_1-{\tilde{S}}_1-{\tilde{S}}_2{\beta }_1\right)}-1,\\ {\alpha }_2=\frac{p}{{\beta }_1}\left(1+{\alpha }_1\right)-1.\end{array}$$

(14)

The final equation, which is obtained from Set (13–14), has the form of Equation (15):

$${\Phi }_P\left({\tilde{S}}_1,{\tilde{S}}_2,p,{\beta }_1\right)=0,$$

(15)

where the values ${\tilde{S}}_1,{\tilde{S}}_2,p$ are assumed to be given. Equation (15) differs from that for the EC [17], shown in Equation (16):

$${\Phi }_S\left({\alpha }_1,{\alpha }_2,s,{\beta }_1\right)=0,$$

(16)

which can be resolved with respect to ${\beta }_1$ for a given ${\alpha }_1,{\alpha }_2$, $s=S_2/S_1={\tilde{S}}_2/{\tilde{S}}_1$. An equilibrium with a given ionization base of the cell, s, could be called an S-equilibrium.

The symmetric version of the EC with subsequent R-equilibrium can be based on recombination processes (for given ${\alpha }_1$ and ${\alpha }_2$). This R-equilibrium is given by a triple ${\alpha }_1^{-1},{\alpha }_2^{-1},\tilde{r}$ with other version of final Equation (17):

$${\Phi }_R\left({\alpha }_1^{-1},{\alpha }_2^{-1},\tilde{r},{\beta }_1\right)=0.$$

(17)

It is important to note that the solution ${\beta }_1$ of Equations (15)–(17) is not a continuous function, since in general there are two equilibrium regions defined by inequalities (18):

$${\kappa }^{\prime }<\kappa <1,$$

(18)

and Equation (19):

$${\kappa }^{\prime }<1<\kappa .$$

(19)

The transition between regions (18) and (19) occurs in jumps in p, x and $v$ [17] at the bifurcation points of Equations (15)–(17). These points are the boundaries of the regions (18) and (19).

For other elementary equilibrium schemes involving one or more cells, the solutions of the final Equations (16) and (17) coincide in principle. However, real impurity equilibrium systems in a plasma with different equilibrium bases (s or $\tilde{r}$ in S- or R-equilibrium respectively) can differ significantly by their equilibrium rate profiles and by the conditions that implement the equilibrium at the variable boundaries. The point is that any of the three options is possible: (i) when the equilibrium can have both equilibrium bases and different profiles, (ii) when only one of the two equilibrium options is realized, and, finally, (iii) when there is no equilibrium for given distributions of parameters. Thus, the modelling shows that the complex cases turn out to have different equilibrium bases and respective boundary conditions.

Thus, the solution of Equation (15) can be used to model the CS equilibrium, which exactly corresponds to the given (or known from experiments) profiles of the total impurity density with known ionization rates. The modelling using this approach are considered below. However, unlike the settings of the direct Equations (16) and (17), which provide the calculated scale λ, the problem (15) assumes that λ is known in advance, e.g., from approximate estimates (which are discussed below) or from a preliminary approximate or direct calculation of the equilibrium.

Direct calculations of the equilibrium of the Ar impurity in the JET tokamak showed [17] that the last cell in the CCs differ from the others in a number of important properties. So, for m-th cell we get ${\alpha }_2\left(EC\right)={\alpha }_{jm+1}\left(CCs\right)$ in the region (19). We also obtain that $\gamma \equiv {\alpha }_2/{\alpha }_1$>>1, s >> 1, ${\beta }_1$<< 1. Typical values are γ = 20–100, s = 4–20, ${\beta }_1=w_1/u_1$ = $2\times {10}^{-4}-2\times {10}^{-3}$. In this case, the analysis can be simplified and, in fact, gives the position of the impurity equilibrium boundary. In the region (19), it turns out that $\kappa \approx {\alpha }_2$, $v<<x\approx 1/{\tilde{R}}_2$, ${\tilde{S}}_1\approx 1/\left(1+p\right)$.

Further, the analysis of the equilibrium in the last EC allows us to calculate the limit values of the boundary densities that satisfy the simple criterion $v\approx 0$ for several values of γ as a function $p\left(s,v=0\right)$. Figure 3 shows the dependence of the limit ratio of densities in the last EC. The equilibrium for a denoted value γ is possible only at densities higher than the calculated limit values.

2.3. Reduction Schemes

The Markovian equilibrium, assumed in GM for stationary plasma conditions, allows us to find the equilibrium transport by reducing GM by the pseudo-state technique. The main reducing schemes are shown in Figure 4.

The most obvious reduction options are local summations of the initial grid distribution function $g_{kn}$ for each n and for each k. It is easy to find several reducing schemes where the nodes are convenient simple sums of CS’s that preserve the original density distribution $\left\{p_n\right\}$ with the normalization. Thus, the formal algebraic procedure allows to reduce the analysis of nonlinear rate transport equations to the solution of the transcendental Equations (15)–(17). So, the GM equilibrium can be successively reduced to a local equilibrium of just four pseudo-states of EC, which are the sums of the corresponding GM quadrants around the single nodes of local EC, that is $G_i={\displaystyle {\sum }_i{g}_{kn}}$, where i = 1, 2, 3, 4 with normalization preserved ${\displaystyle {\sum }_i{G}_i}=1$.

Invariance of the resulting distributions can be illustrated by comparing densities, for example, in simple schemes HL, EC, and SP, where they are connected by obvious equalities shown in Equation (20):

$$\begin{array}{l}{\displaystyle \sum _{i=1}^n{p}_i\left(\mathrm{HL}\right)=}{P}_1(\mathrm{SP})=P_1(\mathrm{EC})=G_1(\mathrm{EC})+G_2(\mathrm{EC}),\\ {\displaystyle \sum _{i=n+1}^{m+1}{p}_i\left(\mathrm{HL}\right)=}{P}_2(\mathrm{SP})=P_2(EC)=G_3(\mathrm{EC})+G_4(\mathrm{EC}).\end{array}$$

(20)

So, for the CCs and HL, the density distribution $\left\{p_n\right\}$ is common, i.e., a reduction invariant. In the same way, the horizontal (HL) and vertical (VL) lines of pseudo-states preserve $\left\{p_n\right\}$ and $\left\{q_k\right\}$ distributions, which are shared with the original GM (see Equations (3) and (4)). Summation of Equation (5) results in the following set of the HL shown in Equation (21):

$$\frac{\partial p_n}{\partial t}=u_{n+1}{p}_{n+1}-\left(u_n+w_n\right)p_n+w_{n-1}{p}_{n-1},$$

(21)

where $u_n={\displaystyle \sum _{k=0}^Z{u}_{kn}^p{f}_{kn}}$ and $w_n={\displaystyle \sum _{k=0}^Z{w}_{kn}^p{f}_{kn}}$, but for the VL to a similar set, shown in Equation (22):

$$\frac{\partial q_k}{\partial t}=R_{k+1}{q}_{k+1}-\left(R_k+S_k\right)q_k+S_{k-1}{q}_{k-1}$$

(22)

where $R_k={\displaystyle \sum _{n=1}^N{R}_{kn}{\phi }_{kn}}$ and $S_k={\displaystyle \sum _{n=1}^N{S}_{kn}{\phi }_{kn}}$.

Another important invariant of the reduction is the CS fluxes between the original neighboring states and their reducing pseudo-states. In this case, the changed transition rates between all new nodes must correspond to the pseudo-state densities. The equalities of probabilistic fluxes determine these rates between pseudo-states. The SP and EC rates are related to each other as follows Equation (23):

$$\begin{array}{l}w\cdot P_1=w_1{G}_1+w_2{G}_2,\\ u\cdot P_2=u_1{G}_4+u_2{G}_3.\end{array}$$

(23)

Obviously, in the case of the CCs scheme, the proper equilibrium is generally provided by the correct choice of the j and j + 1 levels of ionic charge k. It is seen that the choice is determined by the most abundant impurity ions in the plasma, i.e., their corresponding charge and external atomic levels. These j and j + 1 could be called the equilibrium levels of steady state impurity.

The first stage of the GM reduction to the CCs is the summation along vertical lines below and above the levels $j/j+1$ respectively. The resulting CCs pseudo-states are the following sums Equation (24):

$$G_{jn}={\displaystyle \sum _{k=0}^j{g}_{kn}},G_{j+1n}={\displaystyle \sum _{k=j+1}^Z{g}_{kn}}.$$

(24)

The fluxes between levels $j/j+1$ in the original and reducing schemes are assumed to remain unchanged, which allows us to find the rates, shown in Equation (25):

$$S_{jn}=S_{kn}{g}_{kn}/G_{jn},R_{j+1n}=R_{k+1n}{g}_{k+1n}/G_{j+1n}.$$

(25)

The sequence of the ratios derived from the Equations (24) and (25) can be called the equilibrium function, shown in Equation (26):

$${\alpha }_{jn}=G_{j+1n}/G_{jn}=S_{jn}/R_{j+1n},$$

(26)

which determines $f_{jn}\left(\left\{{\alpha }_{jn}\right\}\right)$ and similarly the function ${\phi }_{jn}\left(\left\{{\beta }_{jn}\right\}\right)$, since Equation (27):

$${\beta }_{jn}=G_{jn+1}/G_{jn}.$$

(27)

Besides the density distribution and CS fluxes between neighboring states of original scheme there is another important reduction invariant. It is the equilibrium constant λ discussed in detail in the next section.

3. Equilibrium Constant

The reduction invariants allow us to obtain a general solution to the original transport nonlinear problem. But there is a common for all reducing schemes and, therefore, the most important equilibrium invariant given by the value λ, which should be considered in detail.

Averaging the atomic rates and its ratios over the corresponding φ functions for the n-th reducing EC and for levels $j/j+1$ of CCs, we get Equation (28):

$$\begin{array}{l}\overline{S}\left(n\right)=\frac{S_1\left(n\right)+{\beta }_1\left(n\right)S_2\left(n\right)}{1+{\beta }_1\left(n\right)}={\displaystyle \sum _{n=1}^{m+1}{S}_{jn}{\phi }_{jn}},\\ \overline{R}\left(n\right)=\frac{S_1\left(n\right)+{\beta }_1\left(n\right)S_2\left(n\right)}{{\alpha }_1\left(n\right)+{\alpha }_2\left(n\right){\beta }_1\left(n\right)}={\displaystyle \sum _{n=1}^{m+1}{R}_{j+1n}{\phi }_{j+1n},}\\ \overline{\alpha }=\overline{S}\left(n\right)/\overline{R}\left(n\right)=\frac{{\alpha }_1\left(n\right)+{\alpha }_2\left(n\right){\beta }_1\left(n\right)}{1+{\beta }_1\left(n\right)}.\end{array}$$

(28)

It is seen that $\overline{S}\left(n\right)$ and $\overline{R}\left(n\right)$ remain the same for all reducing cells of the CCs. Then the sum in Equation (29):

$$\lambda \left(n\right)=\overline{S}\left(n\right)+\overline{R}\left(n\right)$$

(29)

is also a common constant. The same conclusions can be drawn in the symmetric case, shown in Equation (30):

$$\lambda \left(k\right)=\overline{w}\left(k\right)+\overline{u}\left(k\right).$$

(30)

The intersection of vertical and horizontal reducing lines and imposed of the EC’s (with the corresponding pairs of n and k) reveals that for any reducing local EC there is only a single common constant $\lambda \left(n\right)=\lambda \left(k\right)=\lambda$. Note that the non-zero eigenvalues of all reducing SP are also equal to λ. In other words, λ is a constant, that determines all possible cases of impurity equilibrium in GM.

The important relationship between λ and the equilibrium transport rates, density distributions, $\left\{p_n\right\}$, and, finally, the characteristic times of the impurity confinement, ${\tau }_p$, can be found as follows.

3.1. Time Scale of Transport Rates

In the case of steady state equilibrium ($\partial p_n/\partial t=0$), the solution of Equation (21) is determined by the equalities of the same particle fluxes from neighboring pseudo-states and between them. Comparing these fluxes for HL and SP and using equation $\lambda =w+u$, we get the equalities corresponding to the conservation of fluxes for each pair of n and n + 1 states of HL and the sequence of SP pseudo-states, shown in Equation (31):

$$\begin{array}{l}w\left(n\right)+u\left(n+1\right)=\lambda ,\\ p_n{w}_n={\displaystyle \sum _{i=1}^n{p}_i}\cdot w\left(n\right)=p_{n+1}{u}_{n+1}={\displaystyle \sum _{i=n+1}^{m+1}{p}_i}\cdot u\left(n+1\right),\end{array}$$

(31)

that gives a complete set of equations. The solution is expressed by Formula (32):

$$\begin{array}{l}{w}_n=\lambda \cdot {\tilde{w}}_n\\ u_{n+1}=\lambda \cdot {\tilde{u}}_{n+1}\end{array}$$

(32)

where ${\tilde{w}}_n$ and ${\tilde{u}}_{n+1}$ are the reduced particle transport rates shown by Formula (33):

$$\begin{array}{l}{\tilde{w}}_n={\displaystyle \sum _{i=1}^n{p}_i}\left(1-{\displaystyle \sum _{i=1}^n{p}_i}\right)/p_n,\\ {\tilde{u}}_{n+1}={\displaystyle \sum _{i=1}^n{p}_i}\cdot \left(1-{\displaystyle \sum _{i=1}^n{p}_i}\right)/p_{n+1}.\end{array}$$

(33)

Using Formula (32), Equation (25) can be rewritten as Equation (34)

$$\frac{\partial \overrightarrow{p}}{\partial t}=N\cdot \overrightarrow{p}=\lambda \cdot \tilde{N}\cdot \overrightarrow{p},$$

(34)

where N is the matrix of Equation (25), $\tilde{N}$ is the reduced matrix combined from the values (33), shown in Formula (35)

$$\tilde{N}\left(\left\{p_n\right\}\right)=\left|\begin{array}{cccc}-{\tilde{w}}_1& {\tilde{u}}_2& 0& 0\\ {\tilde{w}}_1& -\left({\tilde{w}}_2+{\tilde{u}}_2\right)& {\tilde{u}}_3& \dots \\ \dots & \dots & \dots & {\tilde{u}}_{m+1}\\ 0& \dots & {\tilde{w}}_m& -{\tilde{u}}_{m+1}\end{array}\right|=B^{-1}T\cdot \left|\begin{array}{cccc}0& 0& \dots & 0\\ 0& -{\lambda }_1& \dots & 0\\ \dots & \dots & \dots & 0\\ 0& \dots & 0& -{\lambda }_m\end{array}\right|\cdot T^{t}B$$

(35)

where $\left\{0,{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_m\right\}$ is the spectrum of its absolute eigenvalues (let $0<{\lambda }_1<{\lambda }_2<\dots <{\lambda }_m$), determined by the density profile according to Equation (33), the matrix B is diagonal, and T is orthogonal, and T^t is the transposed orthogonal matrix (see, e.g., in [20]).

Thus, the value λ determines the scale $\tau =\lambda \cdot t$ of temporal evolution $\left\{p_n\left(\tau \right)\right\}$. The eigenvalues of $\tilde{N}$ is related only to the stationary profile $\left\{p_n\right\}$, as follows from Equations (33)–(35).

3.2. Impurity Equilibrium and Transport Coefficients

The discrete structure of CS transport, finite equilibrium regions, and bifurcations during transitions between them strongly limit the possibility of describing impurity transport by differential equations that require continuity of the corresponding functions along with their derivatives in the entire domain of definition. On the contrary, the proposed matrix approach uses a minimum of input information to analysis and most closely corresponds to the physical nature of probabilistic (random) CS transport. Indeed, Equations (21) and (34) represent a matrix analogue of the standard 1D continuity equation for the total impurity density shown in Equation (36):

$$\frac{\partial n_Z\left(r,t\right)}{\partial t}+{\nabla }_r{\Gamma }_Z\left(r,t\right)=0$$

(36)

where ${\Gamma }_Z\left(r,t\right)={\displaystyle \sum {\Gamma }_k\left(r,t\right)}$. Note that for a stationary impurity profile $n_Z\left(\rho \right)$ obtained in experiment, the corresponding discrete profile $\left\{p_n\right\}$ is also well known, since $\partial \ln p_n/\partial \rho =\partial \ln n_Z/\partial \rho$, where $\rho =r/a$, $a$ is the minor radius of the plasma column. Next, we note that in the accepted cylindrical geometry, it is convenient to set ${\rho }^2=n/N$ in accordance with Equation (21), and then $N=2\left(m+1\right)$. Comparing the coefficients before the second derivatives in Equation (21) (after expanding by n as a continuous parameter, see, e.g., in [17]) and (36) and using the Equation (32), we obtain the particle diffusion coefficient, D, in Formula (37):

$$D=\frac{a^2}{8{\left(m+1\right)}^3}\lambda \tilde{D}$$

(37)

where $\tilde{D}\left(\left\{p_n\right\}\right)=\left({\tilde{w}}_n+{\tilde{u}}_n\right)/2$ is a reduced diffusivity, which is easily calculated using Equation (33). The reduced convective velocity, $\tilde{v}$, is obtained by a standard way $\tilde{v}=\tilde{D}/a\cdot \partial \ln p_n/\partial \rho$.

The strong dependence $D\propto {\left(m+1\right)}^{-3}$ occurs in Equation (37) for several reasons. First, using the expansion in series on n of Equation (21) we get the second derivative depending on n with corresponding coefficient (see, e.g., in [17,18,19,20]). Second, the comparison can only be made using the initially accepted cylindrical coordinate system with argument ${\rho }^2=n/\left(m+1\right)$, which provides a correct description of the radial states (layers) with equal plasma volumes. Third, unlike Equation (36), the Equations (21) and (34) are linear on n. But the argument ${\rho }^2$ appears before the second derivative in the steady state solution of Equation (36), that is, in the confluent hypergeometric equation to be solved (see, for example, in [22,34,35,36]). Then, to convert the original coefficient before second derivative of linear Equation (21) into that of the standard Equation (36) with that required ${\rho }^2$, it is necessary to include the correcting scale factor ${\left(m+1\right)}^{-1}$ into the expression of D, since the factor $n$ is just included in the diffusivity profile $\tilde{D}\left(n\right)$.

Now we can directly calculate the reduced transport profiles $\tilde{D}$ and $\tilde{v}$ together with the profiles $\left\{{\tilde{w}}_n\right\}$ and $\left\{{\tilde{u}}_n\right\}$ corresponding to the given profile $n_Z\left(\rho \right)$. We can directly derive $\left\{w_n\right\}$ and $\left\{u_n\right\}$, $\left\{{\tilde{w}}_n\right\}$ and $\left\{{\tilde{u}}_n\right\}$ using the invariance of pseudo-state fluxes of HL and EC systems shown in Formulas (38):

$$\begin{array}{l}{p}_n{w}_n=G_1\left(n\right)w_1\left(n\right)+G_2\left(n\right)w_2\left(n\right),\\ p_{n+1}{u}_{n+1}=G_3\left(n\right)u_2\left(n\right)+G_4\left(n\right)u_1\left(n\right)\end{array}$$

(38)

where for each n we can use the following obvious expressions, shown in Formulas (39) and (40):

$$\begin{array}{l}{G}_1=1/\left(1+{\alpha }_1+{\beta }_1+{\alpha }_2{\beta }_1\right),G_2={\alpha }_1/\left(1+{\alpha }_1+{\beta }_1+{\alpha }_2{\beta }_1\right),\\ G_3={\alpha }_2{\beta }_1/\left(1+{\alpha }_1+{\beta }_1+{\alpha }_2{\beta }_1\right),G_4={\beta }_1/\left(1+{\alpha }_1+{\beta }_1+{\alpha }_2{\beta }_1\right),\end{array}$$

(39)

$$\begin{array}{l}{w}_1=S_2{\beta }_1\cdot v,w_2=\frac{S_2{\beta }_1}{{\alpha }_1}\cdot x,\\ u_1=S_2\cdot v,u_2=\frac{S_2}{{\alpha }_1}\cdot x,\end{array}$$

(40)

then from Equations (13) and (38)–(40) we get the relations of the HL and EC values shown by Expressions (41):

$$\begin{array}{l}{w}_n=\frac{S_2{\beta }_1\left(x+v\right)}{p_n\left(1+{\alpha }_1+{\beta }_1+{\alpha }_2{\beta }_1\right)},\\ u_{n+1}=\frac{S_2{\beta }_1\left(x+v\right)}{p_{n+1}\left(1+{\alpha }_1+{\beta }_1+{\alpha }_2{\beta }_1\right)}.\end{array}$$

(41)

Typical profiles $n_Z\left(\rho \right)$, for example, of C and Ar are well known from experiments (see, e.g., in [29,30,31,32,33,34]). The profiles $n_Z\left(\rho \right)$ can be (i) peaked in the plasma core, (ii) almost flat, or (iii) slightly hollow. Figure 5a shows these typical profiles observed in experiments, and Figure 5b shows the corresponding profiles $\left\{{\tilde{w}}_n\right\}$ and $\left\{{\tilde{u}}_n\right\}$ calculated on Equation (33).

Numerically ${\tilde{w}}_n$ and ${\tilde{u}}_n$ close to each other in almost the entire profile. So, we conclude that $\tilde{D}\approx {\tilde{w}}_n\approx {\tilde{u}}_n$ (see Equation (37)). The closer to the center of the plasma, the better the approximate ratio is reproduced. The profiles A₁ and A₂ with central accumulation correspond to a decrease in D in the plasma core, which is in qualitative agreement with the experiments [26] and modelling [17].

3.3. Impurity Equilibrium and Confinement Time

The confinement time of impurity particles, ${\tau }_p$, is conventionally determined by the smallest of the eigenvalues of the divergence operator [34,35], which represents the particle transport in the standard Equation (36). In the proposed matrix case, the analogue of this value, as follows from Equations (34) and (35), is the product $\lambda \cdot {\lambda }_1$. We denote ${\lambda }_E\equiv {\lambda }_1$ the smallest (non-zero) absolute eigenvalue of the spectrum of the matrix $\tilde{N}$ (see explanation to Formula (35)). Then ${\lambda }_E$ can be found from the spectral decomposition $\left\{0,{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_m\right\}$ of the matrix $\tilde{N}\left(\left\{p_n\right\}\right)$ constructed by Formulas (32) and (33). Using this matrix approach, we obtain a relationship between ${\lambda }_E$ and m. Figure 6 shows the calculations of these dependencies ${\lambda }_E\left(m\right)$ for the profiles $n_Z\left(\rho \right)$ shown in Figure 5a.

From the definition of ${\tau }_p$, Equations (34) and (35) and calculations shown in Figure 6, we get Formula (42):

$${\tau }_p=1/\left(\lambda \cdot {\lambda }_1\right)\approx m/\lambda ,$$

(42)

where the last approximate ratio is performed within a good accuracy, especially for flat (F) density profiles. From the expressions (37) and (42) we conclude that $D\propto \lambda \propto 1/{\tau }_p$. The resulting estimates of the value ${\tau }_p$ for the various cases of equilibrium of C and Ar impurities are presented below.

3.4. Types of Impurity CS Equilibrium

The atomic structure of impurity ions is directly manifested in their CS equilibrium. It follows from Formulas (28) and (29) that changes in the equilibrium levels $j/j+1$ can affect both λ and ${\alpha }_{jn}$. An additional factor here is, as mentioned above, the type of the equilibrium base—ionization (S) or recombination (R). This can be shown by modelling of impurity equilibrium.

Figure 7 shows calculations (by the TICS code [17]) of the argon impurity equilibrium for the temperature and density profiles typical for a variety of discharge conditions in JET [31,32,34], in particular, for a plasma with $T_e\left(0\right)=4$ keV, $T_i\left(0\right)=4$ keV, and $n_e\left(0\right)=$4.5 × 10¹⁹ m⁻³. The selected profiles of temperatures, electron density, and relative density of hydrogen neutrals ${\xi }_n=n_n/n_e$ remain the same as previously used for modeling Ar transport in JET.

Figure 7a shows the profiles of the total argon density for three typical cases of the equilibrium observed in experiments. Here, discrete data points are connected by splines. The first profile, denoted by S, is a strongly peaked and obtained in the case of S-equilibrium with the base $j/j+1$ for He/H—like ions. The next profile denoted by R has two maxima (at the center and at the edge), calculated for R-equilibrium with the base $j+1/j$ for He/Li—like ions. Finally, that denoted by R_f, calculated for almost flat in the core and peaked at the plasma edge, also related to R-equilibrium (with He/Li—like ions). The calculations shown in Figure 7b show that a change of the levels He/Li (in R-equilibrium) to He/H (in S-equilibrium) leads to a significant decrease in D. This is obviously due to the low ionization and recombination rates of He-like and H-like ions.

Figure 7c shows the profiles ${\xi }_n\left(\rho \right)$ that determine the rate profiles of recombination through charge-exchange of impurity ions on hydrogen neutrals and, hence, their corresponding equilibrium functions (26). Two regions of these equilibrium profiles for S-equilibrium (He/H) and for R-equilibrium (He/Li) are separated from each other by a narrow transition region depending on the plasma parameters that provide equilibrium. The gradual transition of impurity from R-equilibrium to S-equilibrium related to an increase in the average charge $\langle Z\rangle ={\displaystyle \sum _{k=0}^{Z}k\cdot f_k}$.

This transition requires a noticeable decrease in the charge-exchange rates across the entire plasma cross-section, and especially at the periphery. Even small fluctuations in the plasma parameters (e.g., ${\xi }_n$) at the equilibrium boundary can cause a transition between the S- and R-equilibrium regions along with a significant change in the values of D and ${\tau }_p$. But since this is a transition from the regions (18) and (19), it can only occur in a jump. Figure 7d shows how the corresponding changes in λ occurs when trying to gradually vary the parameters of CS equilibrium.

The modeling calculations of CS equilibrium of Ar were performed for a variety of $T_e\left(0\right)$, $T_i\left(0\right)$ and profiles ${\xi }_n\left(\rho \right)$ close to the modelled [17] conditions at JET [31,32,34]. Thus, the studied CS equilibrium of the types discussed above can be represented in a very wide range of plasma parameters, revealing the corresponding range of λ and ${\tau }_p$.

Figure 8 summarizes these variations of equilibrium values λ (a) and ${\tau }_p$ (b) caused by the change in the temperatures $T_e\left(0\right)$ and $T_i\left(0\right)$ from 1.5 to 8 keV and by the corresponding (to equilibrium) variations of ${\xi }_n\left(\rho \right)$ (from the core to the equilibrium boundary) at an almost constant plasma density profile with $n_e\left(0\right)=4.5\times {10}^{19}{\mathrm{m}}^{-3}$. Full symbols correspond to S-, and open symbols to R-equilibrium type. The growth of $T_e\left(0\right)$ leads to a noticeable increase in λ along with the increase $\overline{\alpha }\left(T_e\right)$. Increase ${\xi }_n\left(0\right)$ changes a lot $\overline{\alpha }\propto 1/\overline{R}\left[{\xi }_n\left(0\right)\right]$, but it does not change much $\lambda =\overline{S}\left(1+1/\overline{\alpha }\right)$, since $\overline{\alpha }\approx 4$. An increase in $T_e\left(0\right)$ results in an increase $\overline{\alpha }\left(T_e,T_i\right)=\overline{S}\left(T_e\right)/\overline{R}\left(T_i\right)$, that is, to a shift of the obtained points on Figure 8 to the right. Figure 8b shows the corresponding calculations of ${\tau }_p$ on Formula (42).

These estimations of ${\tau }_p$ generally correspond to the measurements and known ${\tau }_p$ scaling for L-, I- and H-modes [38,39]. In particular, the estimates ${\tau }_p\approx 5-20$ ms for $j/j+1$ = 15/16 (Li- and He- like Ar) (see in Figure 8b) turn out to be in good agreement with scaling for L-mode plasmas [38] and for I-mode in Alcator C-Mode [39], but taking into account its higher plasma density $n_e\left(0\right)\simeq 2\times {10}^{20}{\mathrm{m}}^{-3}$, since $\lambda \propto n_e$. Moreover, the calculated dependence ${\tau }_p\left(T_e\right)$ clearly shows the same tendency of slow decreasing with strongly increase in $T_e\left(0\right)$ from 1.5 up to 4 keV (see, e.g., Figure 4 in [39]). The transition from L- to H-mode plasmas results in a stepwise increase in ${\tau }_p\approx 50-1000$ ms [38,39,40,41] and corresponds to impurity equilibrium (both S- and R-type) for higher levels $j/j+1$= 16/17 (He- and H- like Ar) and $j/j+1$= 17/18 (H-like and full stripped Ar ions).

3.5. The Radial Scale of the CS Equilibrium

As follows from Formulas (37) and (42), the radial number of equilibrium cells m plays a significant role in the scale factors of the impurity CS transport. The number of radial cells m was initially determined by the number of times on average the particle changes its CS moving from the center to the grid boundary. However, there is an inevitable contribution in m of the other transport processes. Therefore, it is not entirely clear how m is defined. Indeed, it could be assumed that collisional and/or turbulent transport in the plasma can contribute to m along with atomic processes, e.g., as a simple sum $m=m_{at}+m_{coll}+m_{turb}$. The simulations of the CS equilibrium transport show that it is necessary $m\approx 10$ for C, $m\approx 12-14$ for Ar and $m\approx 22-24$ for W. But for C, the largest excess over the reasonable estimates $m_{at}\approx 2-3$ ($\le 6$ at least) is obtained, whereas for medium and heavy impurities, it can be assumed that $m_{coll}+m_{turb}\le 2-5$. It is also worth noting that $\lambda \propto n_e$. But in accordance with Formula (42) $m\approx \lambda {\tau }_p$. So, if we assume that $m_{coll}\propto n_e$, then from $D\propto \lambda /{\left(m+1\right)}^3$, we can finally get a slowly decreasing dependence $D\left(n_e\right)$ as observed in experiments with light impurities (see, e.g., in [12,42]).

Meanwhile, the simulation results are very weakly dependent on the value of m. For example, Figure 9a shows calculations of the dependence $\lambda \left(m\right)$ for two types of equilibrium that differ in the total density profile (see Figure 7a). So, λ is weakly sensitive to m, that is, to the radial dimensions of EC’s over a wide range. Similarly, the value λ is weakly sensitive to $T_e\left(0\right)$, while maintaining the shape of the profile $T_e\left(\rho \right)$. It should be remembered that the ionization rates usually vary slowly near their maxima, depending on the temperature.

4. The Recovery of CS Equilibrium on the Density Profile

As noted above, it is possible to change the formulation of the problem of equilibrium modelling to use the total density profile, $n_Z\left(\rho \right)$, usually known from experiments (or assumed to be given) with the data about ionization rates. Then, the recombination rates could be recovered along with the equilibrium transport rates. Such a formulation can be based on solving Equation (15) after reducing GM to a set of separate REC’s. However, to solve it, it is necessary to have at least approximate data on λ as an input parameter, while the exact value is calculated as Equation (43):

$$\lambda ={\displaystyle \sum _{n=1}^{m+1}\left(S_{jn}{\phi }_{jn}+R_{j+1n}{\phi }_{j+1n}\right)},$$

(43)

where, initially, knowledge of the functions ${\phi }_{jn}$ and ${\phi }_{j+1n}$ is required. Nevertheless, the estimates show that using the total density profile $\left\{p_n\right\}$ (i.e., the profile $n_Z\left(\rho \right)$ taken form the experiment), the approximation ${\lambda }_{est}\approx \lambda$ can be found by the Equation (44)

$${\lambda }_{est}\approx {\displaystyle \sum _{n=1}^{m+1}\left(S_{jn}+R_{j+1n}\right)}\cdot p_n$$

(44)

Calculations of λ using Equations (43) and (44) for C and Ar are compared in Figure 10. It is seen that all these data can be represented within a good accuracy by two approximate relations. For impurity equilibrium within the K-shell, choosing $j/j+1$ of He- and H-like ions, we get ${\lambda }_{est}\left(K\right)=0.68\cdot \lambda$ for both C and Ar, whereas within $j/j+1$ of He- and Li-like ions we obtain ${\lambda }_{est}\left(KL\right)=0.87\cdot \lambda$. So, a reasonable approximation ${\lambda }_{est}\approx \lambda$ to the exact value depends on knowledges of (i) the basis of the equilibrium and (ii) the atomic structure of the most abundant impurity ions. Even rough estimates of the rates of ionization-recombination processes, normalized by an approximate value ${\lambda }_{est}$, provide the necessary input data for modelling using Formula (15) to reduce GM to REC.

Experiments with carbon pellet injection on the LHD stellarator [43] were chosen for modeling, since a large variety of carbon density profiles was observed there. The direct calculation of carbon equilibrium and its quasi-stationary temporal evolution was performed in [17]. These data turn out to be helpful when setting up the proposed recovery model of the profile ${\xi }_n\left(\rho \right)$. So, the recovery model directly uses as input data and reproduces exactly the profiles $n_C\left(\rho \right)$ observed in the LHD. The calculated profiles ${\xi }_n\left(\rho \right)$ (and the corresponding charge-exchange recombination rates) are found by an iterative procedure from the calculations of equilibrium. The convergence of iteration procedure in the recovery case is achieved in 8–10 iterations, which is noticeably faster compared to TICS, where it is achieved in only 30–40 iterations.

Figure 11a shows the experimental profiles $n_C\left(\rho \right)$ adopted to discrete case, that is, $\left\{p_n\right\}$ at $t_1=1.84$ s, $t_2=2.04$ s, $t_3=2.24$ s, $t_4=2.54$ s, which are the input data of the model under consideration, except for the profile with carbon accumulation in the center, denoted by the letter A.

As noted in [17], this profile A was modelled using data from the L-mode at $t_1=1.84$ s and by adjusting the profile ${\xi }_n\left(\rho \right)$, which determines the cx-recombination rates providing 2D equilibrium. Figure 11b shows the diffusion coefficient profiles $D\left(\rho \right)$ calculated by Formula (37). Here we obtain the main feature of equilibrium CS transport in the plasma core, also found in experiments [26], in modelling [17] and revealed by the calculations of the profiles $\left\{{\tilde{w}}_n\right\}$ and $\left\{{\tilde{u}}_n\right\}$ presented above in Figure 5b: the lower the relative impurity density in the core, the higher $D\left(0\right)$. Conversely, the more impurity accumulates in the center, the smaller $D\left(0\right)$ is required for such an equilibrium.

Figure 11c shows the calculation of the quasi-stationary temporal evolution of the equilibrium function, its successive noticeable flattening and displacement of the sharp edge to the periphery of the plasma column along with the maximum of the density profile $n_C\left(\rho \right)$. The resulting profiles ${\xi }_n\left(\rho \right)$ also change sequentially showing a flattening in the plasma core and a decrease in the density of neutrals at the plasma periphery, where the maximum of $n_C\left(\rho \right)$ shifts. It should be noted that, in general, the obtained novel results are in reasonable agreement with the data of direct modeling of the equilibrium CS transport by the TICS code [17].

5. Conclusions

The impurity CS transport is a typical transport process in plasma, which finds its place in a range of other transport processes. This fundamentally 2D transport is consistently represented by symmetric differential and discrete forms of FPKE. The concept of impurity CS equilibrium transport provides an understanding of the crucial role of atomic processes in the behavior of plasma impurities. The impurity equilibrium observed in experiments together with a large variety of impurity density profiles, both total and partial, in a stationary plasma is essentially a double coronal equilibrium (DCE). DCE is a direct consequence of the symmetric properties of the discrete FPKE (7) and of the conservation law of CS’s moving in the proper phase space.

Nonlinear problems of analysis of 2D CS equilibrium and transport are significantly simplified by reducing the original GM to a whole set of simple schemes (see Figure 4) using the pseudo-state technique together with reduction invariants. The most important of them is the equilibrium constant, λ, which determines the general time scale of the impurity equilibrium in a stationary plasma. It provides crucial information to understand the scale of equilibrium transport rates and confinement time of impurities in hot plasmas.

In particular, the difference of the behavior and transport of impurities in plasma is determined by their ionization-recombination rates values that is, by the dependence $\lambda \left(Z\right)$, analysis of which is the immediate task. Naturally, the impurity equilibrium and transport in hot plasma depend on the atomic structure of impurity ions that define the equilibrium levels for the most abundant ions as discussed above.

A significant difference of 3–10 times both in the values of D (see Figure 7b) and in the values ${\tau }_p$ (see Figure 7 and Figure 8) in various types of CS equilibrium suggests their direct coupling with the observations of impurities in the L-, I -, and H-modes in plasma [25,39]. The estimates ${\tau }_p$ made on the basis of CS transport modelling of C and Ar are in good agreement with the known experimental scaling [38] and measurements [39].

The coupling of the profile $\left\{p_n\right\}$ with the rates of particle transport makes it possible to directly analyze typical profiles $\tilde{D}\left(\rho \right)$ and $\tilde{v}\left(\rho \right)$ observed in experiments, in particular, for the central accumulation or any other distribution of impurities in the plasma cross-section. It is shown that using REC analysis, data on $\left\{p_n\right\}$, ionization-recombination rates and the proposed approximation of λ (see Equation (44)), it is also possible to solve the recovery problem for rate profiles along with equilibrium CS transport analysis. In particular, it is possible to calculate the relative profile of neutral hydrogen, which determines the charge-exchange recombination rate of light and medium impurities.

Thus, the concept of impurity CS transport opens up significant prospects for direct analysis of the rates of ionization-recombination processes in a hot quasi-stationary plasma.