Explanation of Experimentally Observed Phenomena in Hot Tokamak Plasmas from the Nonequilibrium Thermodynamics Position

In studying the hot plasma behavior in tokamak devices, the classical approach for collisional processes is traditionally used. This approach leaves unexplained a number of phenomena observed in experiments related to plasma energy confinement. Further, it is well known that tokamak plasma is always turbulent and self-organized. In the present paper, we show that the nonequilibrium thermodynamics approach allows us to explain many observed dependences and paradoxes; for example, puffing of impurities results in confinement improvement if zones of plasma cooling by impurities and additional plasma heating are not overlapped. The analysis of the experimental results shows the important role of radiation losses at the plasma edge in the processes determining its total energy confinement. It is shown that the generally accepted dependence of energy confinement on plasma density is not quite adequate because it is a consequence of dependence on radiation losses. The phenomenon of the appearance of internal transport barriers and magnetic islands can also be explained by plasma self-organization. The obtained results may be taken into account when calculating the operation of a future tokamak reactor.


Introduction
For more than 60 years, the behavior of hot plasma in tokamak devices (toroidal chamber with a strong longitudinal magnetic field and longitudinal current) has been studied. The aim of such research is to create a controlled thermonuclear fusion reactor that does not produce radioactive waste and does not pollute the atmosphere. The relevance of the problem has grown over the years.
Many years of experiments with tokamak devices show that plasma is always turbulent and cannot be described by classical equations for collisional processes, even if we increase the transport coefficients up to some anomalous value. Many interesting phenomena and regimes have been observed. The library of so-called "modes of operation" includes L-mode, H-mode, RI-mode, Super Shot, Advanced Tokamak, Internal Transport Barrier (ITB), and so on.
Many authors have tried to build some fitting models that can describe plasma behavior in some area of experimental conditions. In these models, the flux depends on the gradient of local plasma parameters. However, it appears that for each tokamak and each experimental condition, different models are needed. 2 of 15 In some models, a new term, the so-called "pinch flux", appeared, which is the heat flux to the plasma center enhancing the classical value. This helped in some cases but did not solve all the problems.
It was shown (see, for example, [1]) that many nonlocal phenomena, for example, the energy confinement improvement associated with edge radiation cooling, cannot be described by means of the local flux-gradient expressions.
Powerful modern computers gave hope for calculating processes in plasma from the first-principles approach. The mighty gyrokinetic codes that had to calculate the behavior of all instabilities were created. However, plasma has many degrees of freedom and numerous instabilities and waves, which also interact nonlinearly with each other, and it is impossible to take into account all these badly known processes. Even the best modern computers do not allow us to perform the calculations to describe the magnetohydrodynamic (MHD) processes.
Scientists had to provide the information for the experimental reactor project. They collected much experimental information, but without the possibility of explaining processes by physics, they turned to building a scaling-empirical dependence of the energy confinement time τ E = W d /P in on the discharge parameters. Here, W d is stored plasma energy and P in is injected heating power. The international thermonuclear reactor ITER was designed and built on the base of this scaling [2], but any scaling is only valid in the parameter range where it is found. In our case, the scaling requires extensive extrapolation-more than one order of value. This may be not so dangerous for independent external parameters, but the scaling includes an internal parameter, the plasma density, which depends on τ E itself. Besides that, we need scalings for transition from one mode of operation to another. Difficulties would strongly decrease if we could understand the physics of the processes that are determined by turbulence behavior.
We offer another approach to explain the abovementioned experimental data, which is based on the idea of plasma self-organization. It is well known that tokamak plasma is turbulent and self-organized [3][4][5][6][7]. It was shown experimentally that the normalized pressure profile p N (ρ) is independent of any parameters. Here, ρ = r/(IR/kB) 1/2 , where I is the plasma current, B is the longitudinal magnetic field, k is the plasma cross section elongation, and r is the minor tokamak radius. In contrast with collisional systems, in turbulent plasmas, the particles and their energy cannot pass to the wall with different speeds. Plasma escapes by clusters, so there is no diffusion, only thermal conductivity. An energy balance equation has to hold constant the pressure profile. So, we need to find such an energy balance equation that can support the normalized pressure profile p N (ρ). In [8,9], Dyabilin and Razumova suggested using a Smoluchowski-like equation, which can maintain a given self-organized pressure profile.
In this work, the explanation of the observed experimental phenomena from the plasma self-organization position is offered.

Self-Organization
The self-consistent profile phenomena can be interpreted in the frame of a thermodynamic approach, when the self-consistent solutions correspond to the minimum of the free energy: where S and E are entropy and energy, respectively, and θ is some effective temperature. Following the terminology from [5], we call parameter θ the "magnetic temperature". For the energy balance, we have [8] ∂p ∂t = ∇· θp ξ ∇ ln(p/p c ) + Q s (2) Entropy 2020, 22 where p(r) is the plasma pressure, p c (r) is the self-organized plasma pressure, Q s is a source of heating and radiation cooling, ξ in the Smoluchowski equation describes dissipation in the turbulent system, and the transport coefficient (thermal conductivity) κ = θ/ξ. The heat flow density is given by In the general case, we need to add the neoclassical flux into Equations (2) and (3). For electrons, it is small compared with turbulent flux, but for ions, sometimes it can be significant.
In [8,9], Dyabilin and Razumova considered the magnetic confinement system in general and obtained pressure profiles that depended on some parameter γ. Comparison of these profiles with those measured experimentally in tokamak experiments shows with good accuracy that for tokamaks γ ∝ q L , θ depends on external discharge parameters and the plasma stored energy: θ ∝ p 0 β 0 /q L , where p 0 is the plasma pressure in the center, β 0 = 8πp 0 /B 2 , and q L is the safety factor at the edge.
The first term in brackets in Expression (3) corresponds to the outward gradient flow, and the second term characterizes the thermal pinch, which depends only on the shape of the self-organized profile. The external source of heat and loss terms Qs, distorting the pressure profile, leads to an increase in the level of turbulence and, in turn, increases the transport coefficients. Changes induced by turbulence take place on a short time scale. Experimentally, this is demonstrated, for example, by the speed of the heat wave that propagates along the minor radius after the electron cyclotron resonance (ECR) heating is switched on/off. The time of the heat pulse propagation is two orders shorter than the plasma energy confinement time [10]. The stronger the external influences on the plasma, leading to a distortion in the self-consistent pressure profile, the higher the free energy level (i.e., the turbulent activity), the greater the radial heat flux associated with the distortion of p(r), and the greater the transport coefficient κ = θ/ξ, as it depends on the level of turbulence. Figure 1 shows that the normalized pressure profile p N within the limits of experimental errors is the same for tokamaks with different dimensions, heating methods, profiles, and intensity of heating. We see that the real p N (ρ) is always not far from the self-organized one p c (ρ). This means that the brackets in Equation (3) change slightly and the profile is not totally regulated by the gradient increase but by the change in ξ. The same follows from the calculations of specific experiments. However, even when there is no radial distorting heat flux (i.e., when the plasma free energy has a minimum), there still is a significant level of turbulence and, therefore, a significant value of energy transport.
In the general case, we need to add the neoclassical flux into Equations (2) and (3). For electrons, it is small compared with turbulent flux, but for ions, sometimes it can be significant.
In [8,9], Dyabilin and Razumova considered the magnetic confinement system in general and obtained pressure profiles that depended on some parameter γ. Comparison of these profiles with those measured experimentally in tokamak experiments shows with good accuracy that for tokamaks L q γ ∝ , θ depends on external discharge parameters and the plasma stored energy: where p0 is the plasma pressure in the center, , and qL is the safety factor at the edge.
The first term in brackets in Expression (3) corresponds to the outward gradient flow, and the second term characterizes the thermal pinch, which depends only on the shape of the self-organized profile. The external source of heat and loss terms Qs, distorting the pressure profile, leads to an increase in the level of turbulence and, in turn, increases the transport coefficients. Changes induced by turbulence take place on a short time scale. Experimentally, this is demonstrated, for example, by the speed of the heat wave that propagates along the minor radius after the electron cyclotron resonance (ECR) heating is switched on/off. The time of the heat pulse propagation is two orders shorter than the plasma energy confinement time [10]. The stronger the external influences on the plasma, leading to a distortion in the self-consistent pressure profile, the higher the free energy level (i.e., the turbulent activity), the greater the radial heat flux associated with the distortion of p(r), and the greater the transport coefficient / κ θ ξ = , as it depends on the level of turbulence. Figure 1 shows that the normalized pressure profile pN within the limits of experimental errors is the same for tokamaks with different dimensions, heating methods, profiles, and intensity of heating. We see that the real pN(ρ) is always not far from the self-organized one pc(ρ). This means that the brackets in Equation (3) change slightly and the profile is not totally regulated by the gradient increase but by the change in ξ . The same follows from the calculations of specific experiments.
However, even when there is no radial distorting heat flux (i.e., when the plasma free energy has a minimum), there still is a significant level of turbulence and, therefore, a significant value of energy transport.  [7] and self-consistent profile pc [8] calculated from thermodynamic approach vs normalized radius ρ = r/(IR/B) 1/2 . The total heat flux Γ determined by the heating power deposited into the plasma can be represented as the sum of two fluxes:  [7] and self-consistent profile p c [8] calculated from thermodynamic approach vs normalized radius ρ = r/(IR/B) 1/2 . The total heat flux Γ determined by the heating power deposited into the plasma can be represented as the sum of two fluxes: where Γ 0 is the heat flux associated with the self-consistent pressure profile, where the turbulence level is minimal, and Γ 1 is the part of flux Γ associated with the pressure profile distortion. In order to accumulate the entire energy of additional heating in the plasma with the best confinement, we must realize a deposition profile that coincides with the profile of the self-consistent pressure. Then, Γ 1 = 0, and it does not spoil the confinement. This option can be implemented more or less in ohmic mode, if the radiation loss does not significantly distort the profile of the deposited power. In such a case, the current density profile j(r) will be proportional to T 3/2 e , if the aspect ratio is not too low. Being self-organized, plasma has to form such T e (r) and n e (r) profiles, where P OH (r) coincides with p c (r). In this case, This experimental result is presented in Figure 2. Being generated in the region of the local heating, Γ 1 flows along the radius, distorting p c (r) and increasing transport. We can express the thermal conductivity coefficient as where χ = 1/ξ = χ 0 + χ 1 and θ χ has a dimension of diffusion.
Entropy 2019, 21, x FOR PEER REVIEW 4 of 15 where Γ0 is the heat flux associated with the self-consistent pressure profile, where the turbulence level is minimal, and Γ1 is the part of flux Γ associated with the pressure profile distortion. In order to accumulate the entire energy of additional heating in the plasma with the best confinement, we must realize a deposition profile that coincides with the profile of the self-consistent pressure. Then, Γ1 = 0, and it does not spoil the confinement. This option can be implemented more or less in ohmic mode, if the radiation loss does not significantly distort the profile of the deposited power. In such a case, the current density profile j(r) will be proportional to 3 2 / e T , if the aspect ratio is not too low. Being self-organized, plasma has to form such Te(r) and ne(r) profiles, where POH(r) coincides with pc(r). In this case, This experimental result is presented in Figure 2. Being generated in the region of the local heating, Γ1 flows along the radius, distorting pc(r) and increasing transport. We can express the thermal conductivity coefficient as

Most Unclear Experimental Fact: RI mode
Among the many unexplained phenomena is the dependence of the energy confinement time , then it saturates. It is included in the ITER scaling. It seems to be very important, since the physical cause of this dependence remained unclear.
To understand this phenomenon, special experiments were performed in the T-10 tokamak (major radius of toroidal plasma R = 1.5 m, minor radius a = 0.3 m), where, in addition to ohmic heating (OH), plasma was heated by ECR power [12]. It has been shown that in conditions where the plasma density is maintained constant by a feedback system (Figure 4), the stored energy Wd is independent of the plasma density ( Figure 4a); Wd may be increased by 40% due to the Ne injection with a constant density (within 10%). In Figure 4b, the stored energy Wd is shown as a function of radiation losses Prad for the same shots. The form of this curve is typical for the ( ) E e n τ dependence shown in Figure 3.
After the initial increase of ( ) E rad P τ in the left part of the curve (while the plasma density is the same),

Most Unclear Experimental Fact: RI mode
Among the many unexplained phenomena is the dependence of the energy confinement time τ E = W d /P in on the plasma density observed experimentally ( Figure 3) [11]. At low densities, τ E ∝ n e , then it saturates. It is included in the ITER scaling. It seems to be very important, since the physical cause of this dependence remained unclear.
To understand this phenomenon, special experiments were performed in the T-10 tokamak (major radius of toroidal plasma R = 1.5 m, minor radius a = 0.3 m), where, in addition to ohmic heating (OH), plasma was heated by ECR power [12]. It has been shown that in conditions where the plasma density is maintained constant by a feedback system (Figure 4), the stored energy W d is independent of the plasma density ( Figure 4a); W d may be increased by 40% due to the Ne injection with a constant density (within 10%). In Figure 4b, the stored energy W d is shown as a function of radiation losses P rad for the same shots. The form of this curve is typical for the τ E (n e ) dependence shown in Figure 3. After the initial increase of τ E (P rad ) in the left part of the curve (while the plasma density is the same), it becomes independent of P rad . Thus, we observed that the plasma stored energy is strongly dependent on radiation losses and nearly independent of the plasma density. it becomes independent of Prad. Thus, we observed that the plasma stored energy is strongly dependent on radiation losses and nearly independent of the plasma density. The effect of confinement increases with the increase of radiation losses after different gas injections into the plasma was discovered in the 1990s by Ongena et al. and Messiaen and Ongena [13,14]. This phenomenon has been named "RI-mode".  In Figure 5, the dependence of Prad on the density for the OH case is presented. In the main part of the curve, Prad is proportional to ne. Why? The injection of the main working gas H2 or D2 leads to a strong charge-exchange flux to the wall and an increase of the light impurity influx into the plasma. For H2, it has to be less intensive than for D2. Hence, when injecting working gas, we also injected light impurities. The intensity of the impurity flux depends on the charge-exchange flux, the kind of hydrogen isotope, and the wall conditions. it becomes independent of Prad. Thus, we observed that the plasma stored energy is strongly dependent on radiation losses and nearly independent of the plasma density. The effect of confinement increases with the increase of radiation losses after different gas injections into the plasma was discovered in the 1990s by Ongena et al. and Messiaen and Ongena [13,14]. This phenomenon has been named "RI-mode".  In Figure 5, the dependence of Prad on the density for the OH case is presented. In the main part of the curve, Prad is proportional to ne. Why? The injection of the main working gas H2 or D2 leads to a strong charge-exchange flux to the wall and an increase of the light impurity influx into the plasma. For H2, it has to be less intensive than for D2. Hence, when injecting working gas, we also injected light impurities. The intensity of the impurity flux depends on the charge-exchange flux, the kind of hydrogen isotope, and the wall conditions. The effect of confinement increases with the increase of radiation losses after different gas injections into the plasma was discovered in the 1990s by Ongena et al. and Messiaen and Ongena [13,14]. This phenomenon has been named "RI-mode".
In Figure 5, the dependence of P rad on the density for the OH case is presented. In the main part of the curve, P rad is proportional to n e . Why? The injection of the main working gas H 2 or D 2 leads to a strong charge-exchange flux to the wall and an increase of the light impurity influx into the plasma. For H 2 , it has to be less intensive than for D 2 . Hence, when injecting working gas, we also injected obviously depends on plasma energy confinement as energy and density in a turbulent plasma move to the wall with the same speed. How we can explain τE dependence on Prad?
If the radiation loss Prad reduces the flux Γ1 (Equation (3)), it should lead to a decrease in 1 1 ( ) χ Γ in the region where cooling takes place and outside it. Hence, the confinement improves because the plasma is self-organized, and an increase in the edge gradient ∇ p means its growth over the whole plasma cross section. When the increase of the radiation losses leads to 1 / (2 2 ) rad P R r π π ≈ Γ , then 1 1 ( ) χ Γ will be small, and the improvement in confinement with increasing Prad saturates. A further increase in Prad will lead to a distortion of the pressure profile and the excitation of instabilities. So, χ0 and χ1 in Equation (6) have somewhat different physical senses: χ0 is a value determined by the character of the turbulent flux, which plasma uses in its self-organization when the free energy is minimal. It cannot be decreased in any case and is independent of any parameters. χ1 depends on the distortion factor-the value and profile of the input power-which cannot be deposited into the self-consistent profile, but it is also independent of the other discharge parameters.
The coefficient / κ θ ξ = should be independent of local parameters ne(r), Te(r), and so on; otherwise, there can be no self-organization. All initial instabilities are important for turbulence development, but after the self-organization-mixing and averaging-the direct link between these instabilities and confinement disappears. It exists only for a time scale shorter than the turbulent relaxation time, which may be characterized by the heat/cold wave propagation time and is two orders of value less than τE. The coefficient ( ) can be determined from experiments with various Pin. Details of χ0 and χ1 determination are discussed in Section 5.

Dependence of Energy Confinement on Radiation Losses
It was shown in [14] that Ne, Ar, or Si injection or combinations yields the same improvement in energy confinement. It is important to check that this is valid also for He, which has Z = 2 and radiates at the very edge of the plasma. In T-10 experiments [12], it was shown that Wd saturated at the same level for Ne and He injection for other similar conditions. Plasma in this tokamak was bounded with a circular metal limiter. It was possible, besides the usual working gas D2, to inject the noble gases Ne and He. Increasing the D2 injection, we also increased the fluxes of C and O from the wall. It was also possible to decrease the flux of these light impurities by lithiation of the wall and limiter surfaces, which strongly decreased the radiation losses. Figure 6 shows the dependence of the energy confinement time τE on the intensity of radiation losses for ohmic discharges with two different magnetic fields. The curves in Figure 6 may be divided We see that τ E is independent of plasma density. In contrast, however, the plasma density obviously depends on plasma energy confinement as energy and density in a turbulent plasma move to the wall with the same speed. How we can explain τ E dependence on P rad ?
If the radiation loss P rad reduces the flux Γ 1 (Equation (3)), it should lead to a decrease in χ 1 (Γ 1 ) in the region where cooling takes place and outside it. Hence, the confinement improves because the plasma is self-organized, and an increase in the edge gradient ∇p means its growth over the whole plasma cross section. When the increase of the radiation losses leads to P rad /(2πR 2πr) ≈ Γ 1 , then χ 1 (Γ 1 ) will be small, and the improvement in confinement with increasing P rad saturates. A further increase in P rad will lead to a distortion of the pressure profile and the excitation of instabilities.
So, χ 0 and χ 1 in Equation (6) have somewhat different physical senses: χ 0 is a value determined by the character of the turbulent flux, which plasma uses in its self-organization when the free energy is minimal. It cannot be decreased in any case and is independent of any parameters. χ 1 depends on the distortion factor-the value and profile of the input power-which cannot be deposited into the self-consistent profile, but it is also independent of the other discharge parameters.
The coefficient κ = θ/ξ should be independent of local parameters n e (r), T e (r), and so on; otherwise, there can be no self-organization. All initial instabilities are important for turbulence development, but after the self-organization-mixing and averaging-the direct link between these instabilities and confinement disappears. It exists only for a time scale shorter than the turbulent relaxation time, which may be characterized by the heat/cold wave propagation time and is two orders of value less than τ E . The coefficient (χ 0 + χ 1 (Γ 1 )) can be determined from experiments with various P in . Details of χ 0 and χ 1 determination are discussed in Section 5.

Dependence of Energy Confinement on Radiation Losses
It was shown in [14] that Ne, Ar, or Si injection or combinations yields the same improvement in energy confinement. It is important to check that this is valid also for He, which has Z = 2 and radiates at the very edge of the plasma. In T-10 experiments [12], it was shown that W d saturated at the same level for Ne and He injection for other similar conditions. Plasma in this tokamak was bounded with a circular metal limiter. It was possible, besides the usual working gas D 2 , to inject the noble gases Ne and He. Increasing the D 2 injection, we also increased the fluxes of C and O from the wall. It was also possible to decrease the flux of these light impurities by lithiation of the wall and limiter surfaces, which strongly decreased the radiation losses.  Figure 6 shows the dependence of the energy confinement time τ E on the intensity of radiation losses for ohmic discharges with two different magnetic fields. The curves in Figure 6 may be divided into three parts. The first part for low P rad is characterized by a linear confinement increase with P rad ; in the second part, we see independence of τ E from P rad , while too high P rad intensity leads to strong cooling of the plasma edge, pressure profile distortion, excitation of MHD instabilities, and τ E degradation. This effect is stronger for lower q L (shots with B = 1.9 T) because p c (r) distortion takes place near more important low q = m/n number rational surfaces.
injected. The dependence of plasma confinement on the amount of injected impurities, that is, on the radiation value at the plasma edge, was studied.
In Figure 7, we show the results obtained with strong ECR heating (PEC = 0.85 MW). Without Ne (curve 1), the plasma stored energy Wd was a little bit lower than for ohmic discharge. Applying Ne injection, as expected, we increased (restored) the stored energy. The stronger Ne injection led to a further improvement in the energy confinement (curves (2) and (3)).
An examination of these results with Equation (2) confirms the previous conclusion that the plasma energy confinement is independent of the plasma density (the initial densities for all shots were nearly the same) but depends on the radiation losses. An increase of Γ1 led to an increase in χ1, which, in turn, depended on Γ1. Ne injection reduced Γ1 in the outer part of the plasma and created an improved confinement zone. Although a small amount of Ne did not lead to a significant increase in the number of electrons in the plasma, the plasma density increased due to a decrease in χ1 at the edge.  According to the explanation given above, in the first part of the curve with linear increase of τ E (P rad ), the flux Γ 1 distorting the plasma pressure profile is partly compensated by P rad . In the second part of the curve, the flux Γ 1 is compensated totally in the radiation zone and outside it. In this case, Γ 1 = 0 and χ is minimal; χ = χ 0 . This is the best possible confinement (without the transport barrier, ITB). The difference in confinement for two series of shots is in accordance with θ(q L ) dependence.
Once the increase in radiation losses leads to an improvement in confinement, a decrease in radiation losses should lead to degradation of confinement and loss of a significant fraction of the plasma. Using lithiation of the walls and limiter in the discharge chamber, it was possible to strongly decrease radiation losses. In the following experiments, the dependence of plasma confinement on radiation losses at the edge under strong additional ECR heating (by a factor of 8 greater than ohmic heating) and a wide change of P rad was studied. The scenario of the experiment was as follows: A series of shots with the same parameters was carried. At the ohmic stage of the discharge at time t = 500 ms, additional on-axis ECR heating with a power of 0.85 MW was used. At t = 600 ms, Ne was injected. The dependence of plasma confinement on the amount of injected impurities, that is, on the radiation value at the plasma edge, was studied.
In Figure 7, we show the results obtained with strong ECR heating (P EC = 0.85 MW). Without Ne (curve 1), the plasma stored energy W d was a little bit lower than for ohmic discharge. Applying Ne injection, as expected, we increased (restored) the stored energy. The stronger Ne injection led to a further improvement in the energy confinement (curves (2) and (3)).
An examination of these results with Equation (2) confirms the previous conclusion that the plasma energy confinement is independent of the plasma density (the initial densities for all shots were nearly the same) but depends on the radiation losses. An increase of Γ 1 led to an increase in χ 1 , which, in turn, depended on Γ 1 . Ne injection reduced Γ 1 in the outer part of the plasma and created an improved confinement zone. Although a small amount of Ne did not lead to a significant increase in the number of electrons in the plasma, the plasma density increased due to a decrease in χ 1 at the edge.  (1) is for a discharge without Ne injection. Despite the powerful heating, the stored energy is lower than for the ohmic discharge. Curves (2) and (3) show the improvement of the energy confinement for shots with stronger Ne injection.
As we need to use some gas injection for confinement increase, we have to know, should these ions accumulate in the hot core of plasma or it is possible to avoid this effect?
In ohmic discharges with low radiative losses in the plasma center, the profile of the input power has almost the same shape as the pressure profile, and the flux Γ1 is close to zero. Thus, for those conditions, the best energy confinement is obtained (in this case, Wd does not increase with Prad). However, this also means that the pinch flux (second term in Equation (3)) is everywhere nearly comparable to the gradient flux, and impurities are able to penetrate into the hot plasma core. In such plasmas, the addition of a neoclassical flux may be important for ions. If their collision frequency is high, they may be not involved in turbulent movement and their confinement will be determined by usual collisional equations. In such a neoclassical case, ions will be distributed in plasma with a law: In the "clean" plasma at the beginning of the ( ) E rad P τ dependence, Zeff = 1, but the energy confinement is poor. With increasing radiation losses at the edge, the pinch flow in Equation (3) will increase, the impurity ions will accumulate in the plasma core, and Zeff will increase. Its value depends on the intensity of the radiating gas injection. The highest Zeff will be reached when the plasma stored energy ( ) d r a d W P reaches saturation due to the insertion into the plasma of more radiating ions; then, the energy confinement reaches its maximum. In this case, the radial distribution of the impurity ions and the working gas will be the same due to the specifics of turbulent transport (plasma participates in the motion "as a whole"). So, the effective plasma charge ( ) eff Z r const = . This is supported by experimental results [15]. This may be a problem for reactors, since we will be forced to inject radiating impurities into the plasma volume, but we do not want to have strong bremsstrahlung radiation.

Calculations of χ0 and χ1
During the impurity gas puffing, the increase of radiation losses at the plasma edge results in a decrease of the radial heat flux Γ1 (up to Γ1 = 0); hence, the coefficient χ in the plasma edge reduces until κ reaches the minimum value 0 κ θχ = . As a result, the plasma stored energy increases and reaches the saturation level. In order to find χ0 experimentally, we have to match the χ0 value in Equation (3) (1) is for a discharge without Ne injection. Despite the powerful heating, the stored energy is lower than for the ohmic discharge. Curves (2) and (3) show the improvement of the energy confinement for shots with stronger Ne injection.
As we need to use some gas injection for confinement increase, we have to know, should these ions accumulate in the hot core of plasma or it is possible to avoid this effect?
In ohmic discharges with low radiative losses in the plasma center, the profile of the input power has almost the same shape as the pressure profile, and the flux Γ 1 is close to zero. Thus, for those conditions, the best energy confinement is obtained (in this case, W d does not increase with P rad ). However, this also means that the pinch flux (second term in Equation (3)) is everywhere nearly comparable to the gradient flux, and impurities are able to penetrate into the hot plasma core. In such plasmas, the addition of a neoclassical flux may be important for ions. If their collision frequency is high, they may be not involved in turbulent movement and their confinement will be determined by usual collisional equations. In such a neoclassical case, ions will be distributed in plasma with a law: n z (r) ∝ n z p (r), where n z and n p are densities of ions with charge z and protons correspondingly. In the "clean" plasma at the beginning of the τ E (P rad ) dependence, Z eff = 1, but the energy confinement is poor. With increasing radiation losses at the edge, the pinch flow in Equation (3) will increase, the impurity ions will accumulate in the plasma core, and Z eff will increase. Its value depends on the intensity of the radiating gas injection. The highest Z eff will be reached when the plasma stored energy W d (P rad ) reaches saturation due to the insertion into the plasma of more radiating ions; then, the energy confinement reaches its maximum. In this case, the radial distribution of the impurity ions and the working gas will be the same due to the specifics of turbulent transport (plasma participates in the motion "as a whole"). So, the effective plasma charge Z e f f (r) = const. This is supported by experimental results [15].
This may be a problem for reactors, since we will be forced to inject radiating impurities into the plasma volume, but we do not want to have strong bremsstrahlung radiation.

Calculations of χ 0 and χ 1
During the impurity gas puffing, the increase of radiation losses at the plasma edge results in a decrease of the radial heat flux Γ 1 (up to Γ 1 = 0); hence, the coefficient χ in the plasma edge reduces until κ reaches the minimum value κ = θ χ 0 . As a result, the plasma stored energy increases and reaches the saturation level.
In order to find χ 0 experimentally, we have to match the χ 0 value in Equation (3)  not overlap with the area of the radiation losses, χ 0 can also be calculated with sufficient precision. However, for ohmic discharges, where the power is deposited over the whole plasma cross section, the accuracy may be lower due to the less precise measurement of the radiation losses. Figure 8 shows measured value χ 0 for T-10 experiments with Ne puffing, ECR heating, and ohmic heating. We see that within the experimental error, χ 0 is independent of the heating power input. the accuracy may be lower due to the less precise measurement of the radiation losses. Figure 8 shows measured value χ0 for T-10 experiments with Ne puffing, ECR heating, and ohmic heating. We see that within the experimental error, χ0 is independent of the heating power input. Γ = Γ + Γ , as χ 1 depends on Γ1. For each case, we have to calculate Γ 0, using the already known χ0, to find Γ1 and then match χ1(Γ1/Γ0) so that the given flux Γ corresponds to the measured value of Wd. This allows us to obtain the value ( ) k θ χ χ = + Γ Γ . If we do this for different heating fluxes, the dependence χ (Γ1/Γ0) may be obtained. The absolute values of χ0 and χ1 depend on the self-organized pressure profile pc, which we have to put into the Equations (2) and (3). As we do not have a reliable experimental profile, we used the theoretical profile from [8], shown in Figure  1.

Internal Transport Barriers (ITBs)
Local zones with enhanced confinement were discovered in tokamak plasmas at the beginning of the 1990s [16]. These zones are referred to as ITBs. Experimental investigation of the ITB has shown that in the ITB zone, the pressure gradients can largely exceed the gradient dictated by the profile consistency. So, existence of ITBs seems to be in contradiction with the self-consistent pressure profiles.
The following question naturally arises: what is the correspondence between profile consistency and ITBs? Are they two independent experimental facts, or are they in some way connected to each other?
From experimental results, it was concluded that a necessary condition for the appearance of the reduced transport zone is to have a low absolute value of / dq dr in the vicinity of the rational qsurface with low poloidal m and toroidal n numbers. The safety factor q characterizes a helicity of magnetic field lines, / q m n = . This result was found in experiments on many tokamaks: TCV [17], T-10 [18][19][20], and JET [21][22][23]. We know that rational surfaces have a specific feature. In principal, the magnetic surfaces fulfill all the space, but if we limit m by the upper number m1 (i.e., m < m1), then in the vicinity of the surface with low m and n, some area without magnetic surfaces exists, which is the so-called Gap. Inside this Gap, surfaces with m > m1 up to infinity exist, but there are no surfaces with m < m1. Simple analysis [20] shows that the radial width of the Gap is Determination of (χ 0 + χ 1 ) is more complicated. The value of κ = p 2 0 IR(χ 0 + χ 1 )/(B 3 0 a 2 ) itself depends on Γ = Γ 1 + Γ 0 , as χ 1 depends on Γ 1 . For each case, we have to calculate Γ 0 , using the already known χ 0 , to find Γ 1 and then match χ 1 (Γ 1 /Γ 0 ) so that the given flux Γ corresponds to the measured value of W d . This allows us to obtain the value k = θ(χ 0 + χ 1 (Γ 1 /Γ 0 )). If we do this for different heating fluxes, the dependence χ (Γ 1 /Γ 0 ) may be obtained. The absolute values of χ 0 and χ 1 depend on the self-organized pressure profile p c , which we have to put into the Equations (2) and (3). As we do not have a reliable experimental profile, we used the theoretical profile from [8], shown in Figure 1.

Internal Transport Barriers (ITBs)
Local zones with enhanced confinement were discovered in tokamak plasmas at the beginning of the 1990s [16]. These zones are referred to as ITBs. Experimental investigation of the ITB has shown that in the ITB zone, the pressure gradients can largely exceed the gradient dictated by the profile consistency. So, existence of ITBs seems to be in contradiction with the self-consistent pressure profiles.
The following question naturally arises: what is the correspondence between profile consistency and ITBs? Are they two independent experimental facts, or are they in some way connected to each other?
From experimental results, it was concluded that a necessary condition for the appearance of the reduced transport zone is to have a low absolute value of dq/dr in the vicinity of the rational q-surface with low poloidal m and toroidal n numbers. The safety factor q characterizes a helicity of magnetic field lines, q = m/n. This result was found in experiments on many tokamaks: TCV [17], T-10 [18][19][20], and JET [21][22][23].
We know that rational surfaces have a specific feature. In principal, the magnetic surfaces fulfill all the space, but if we limit m by the upper number m 1 (i.e., m < m 1 ), then in the vicinity of the surface with low m and n, some area without magnetic surfaces exists, which is the so-called Gap. Inside this Gap, surfaces with m > m 1 up to infinity exist, but there are no surfaces with m < m 1 . Simple analysis [20] shows that the radial width of the Gap is δ gap ∝ q/(m 1 dq/dr).
(7) Figure 9 shows the result of the rational surface density calculation for the MAST tokamak (low aspect ratio R/a = 1.5). We can see the electron temperature profile with ITB at q = 1 and rational surface density for two m 1 values: m 1 = 20 and m 1 = 30. The lower m 1 leads to a wider Gap. The turbulent heat flux experiences difficulties passing the Gap region. This means that the turbulent structure elements are linked with magnetic surfaces (may be with fluctuating current filaments). Cutting off the lower, the most powerful modes result in limiting the κ increase, more precisely, in χ increase, but κ, as we know, regulates the pressure profile near the self-consistent one for given Γ. If the external force distorts p(ρ), the brackets in Equation (3)  (7) Figure 9 shows the result of the rational surface density calculation for the MAST tokamak (low aspect ratio R/a = 1.5). We can see the electron temperature profile with ITB at q = 1 and rational surface density for two m1 values: m1 = 20 and m1 = 30. The lower m1 leads to a wider Gap. The turbulent heat flux experiences difficulties passing the Gap region. This means that the turbulent structure elements are linked with magnetic surfaces (may be with fluctuating current filaments). Cutting off the lower, the most powerful modes result in limiting the κ increase, more precisely, in χ increase, but κ, as we know, regulates the pressure profile near the self-consistent one for given Γ. If the external force distorts p(ρ), the brackets in Equation (3) must grow as a result of 1/p·dp/dr growth, and the free energy must grow too. In what form and in what structure will it appear in plasma? To understand the conditions for ITB appearance, it is necessary to know which lower turbulence mode number m1 is needed for the radial transfer of heat flux. Using Relation (7) and comparing the calculated Gap width with the barrier width for the large number of published results for different tokamaks, it was possible to find a rough dependence linking the radial dependence of the heat flux density with m1. The result is shown in Figure 10 [20]. Despite the roughness of such an estimate (we do not know the radial correlation length of a given mode), the points obtained at a variety of tokamaks were well grouped on a steep dependence: m1 should decrease with increasing the flux density Γ, that is, , where 3/2 < α < 5/2. This means that when we increase Γ, the lower mode in the turbulent spectrum can get into the Gap near some rational surface, and if the derivative dq/dr is sufficiently low, it leads to the formation of a barrier.  To understand the conditions for ITB appearance, it is necessary to know which lower turbulence mode number m 1 is needed for the radial transfer of heat flux. Using Relation (7) and comparing the calculated Gap width with the barrier width for the large number of published results for different tokamaks, it was possible to find a rough dependence linking the radial dependence of the heat flux density with m 1 . The result is shown in Figure 10 [20]. Despite the roughness of such an estimate (we do not know the radial correlation length of a given mode), the points obtained at a variety of tokamaks were well grouped on a steep dependence: m 1 should decrease with increasing the flux density Γ, that is, Γ ∼ m −α 1 , where 3/2 < α < 5/2. This means that when we increase Γ, the lower mode in the turbulent spectrum can get into the Gap near some rational surface, and if the derivative dq/dr is sufficiently low, it leads to the formation of a barrier. (7) Figure 9 shows the result of the rational surface density calculation for the MAST tokamak (low aspect ratio R/a = 1.5). We can see the electron temperature profile with ITB at q = 1 and rational surface density for two m1 values: m1 = 20 and m1 = 30. The lower m1 leads to a wider Gap. The turbulent heat flux experiences difficulties passing the Gap region. This means that the turbulent structure elements are linked with magnetic surfaces (may be with fluctuating current filaments). Cutting off the lower, the most powerful modes result in limiting the κ increase, more precisely, in χ increase, but κ, as we know, regulates the pressure profile near the self-consistent one for given Γ. If the external force distorts p(ρ), the brackets in Equation (3)  To understand the conditions for ITB appearance, it is necessary to know which lower turbulence mode number m1 is needed for the radial transfer of heat flux. Using Relation (7) and comparing the calculated Gap width with the barrier width for the large number of published results for different tokamaks, it was possible to find a rough dependence linking the radial dependence of the heat flux density with m1. The result is shown in Figure 10 [20]. Despite the roughness of such an estimate (we do not know the radial correlation length of a given mode), the points obtained at a variety of tokamaks were well grouped on a steep dependence: m1 should decrease with increasing the flux density Γ, that is, , where 3/2 < α < 5/2. This means that when we increase Γ, the lower mode in the turbulent spectrum can get into the Gap near some rational surface, and if the derivative dq/dr is sufficiently low, it leads to the formation of a barrier.  Taking into account that for ITB formation we need to organize a wide Gap, we performed experiments, in which the ITB was generated in a given place with a given q [24]. The barrier near the surface with q = 3/2 was generated using off-axis ECR heating and a rapidly increasing plasma current. In Figure 11, the time evolution of the electron temperature T e for three plasma radii during ITB formation is given with three different time scales.
Taking into account that for ITB formation we need to organize a wide Gap, we performed experiments, in which the ITB was generated in a given place with a given q [24]. The barrier near the surface with q = 3/2 was generated using off-axis ECR heating and a rapidly increasing plasma current. In Figure 11, the time evolution of the electron temperature Te for three plasma radii during ITB formation is given with three different time scales.
In Figure 12a, the scenario of the experiment is given. The ITB appeared after current ramp-up due to the reduction of dq/dr at r = 16 cm (Figure 12b). The temperature rose inside the zone with r = 16 cm. At t = 730 ms, the rise stopped and harmonic oscillations typical of the magnetic island with m/n = 3/2 appeared. The island grew, and at t = 763 ms, the first magnetic reconnection occurred with the release of energy outwards-this is the so-called "internal disruption". This had to decrease the level of free energy. The process was repeated until the current noticeably penetrated from the periphery of the plasma inside its core, destroying the barrier. In Figure 12a, the scenario of the experiment is given. The ITB appeared after current ramp-up due to the reduction of dq/dr at r = 16 cm (Figure 12b). The temperature rose inside the zone with r = 16 cm. At t = 730 ms, the rise stopped and harmonic oscillations typical of the magnetic island with m/n = 3/2 appeared. The island grew, and at t = 763 ms, the first magnetic reconnection occurred with the release of energy outwards-this is the so-called "internal disruption". This had to decrease the level of free energy. The process was repeated until the current noticeably penetrated from the periphery of the plasma inside its core, destroying the barrier.
So, we may conclude that the enhanced free energy (Equation (1)) at the top of the ITB leads to the formation of a new (fractal?) structure: the current island.
Our hypothesis is as follows: The heat flux coming from the inside increases ∇p in the Gap region. The heat flux supporting the self-consistent pressure profile should be based on pressure fluctuations extended along magnetic field lines. However, such fluctuations are associated with fluctuations in the bootstrap current, so we should expect the appearance of fluctuating current filaments. Overlapping in space, they can transfer the energy along the radius. The larger the flow, the larger the filaments (lower m) that should be generated; that is, the spectrum should expand toward low m. This is consistent with Figure 10. Deformation of the pressure profile in the Gap increases the free energy and decreases m, transforming the narrow filament into the current island.
We suppose that the radial flux is organized by small current filaments fluctuating in a radial direction. These filaments, accumulating free energy, appear due to deformation of the pressure profile in the Gap, forming current islands. Depending on the value of the stored free energy, they will have larger or smaller values of m and n, preserving the value q = m/n. This has been well demonstrated in beautiful KSTAR experiments ( Figure 13) [25]. In these experiments, the fine T e (and so, the current density) structure was registered from the image of the second EC harmonic emission. This gave the fine 2D temperature profile in some core plasma region with q ≤ 1. Note that the current, flowing inside this region, forms a spiral with m/n = 1/1 (kink mode), and the position of the EC resonance relative to the current center can be different. The geometry is complicated, and at various ECR heating positions relative to the current center, we may obtain q = 1 with m/n from 1/1 up to 3/3. surfaces calculated for m1 = 30 and q profiles for different time moments: in OH phase, before current ramp-up, and after current ramp-up with ECR heating. In the latter case, a large gap marked by ellipse is seen.
So, we may conclude that the enhanced free energy (Equation (1)) at the top of the ITB leads to the formation of a new (fractal?) structure: the current island.
Our hypothesis is as follows: The heat flux coming from the inside increases ∇ p in the Gap region. The heat flux supporting the self-consistent pressure profile should be based on pressure fluctuations extended along magnetic field lines. However, such fluctuations are associated with fluctuations in the bootstrap current, so we should expect the appearance of fluctuating current filaments. Overlapping in space, they can transfer the energy along the radius. The larger the flow, the larger the filaments (lower m) that should be generated; that is, the spectrum should expand toward low m. This is consistent with Figure 10. Deformation of the pressure profile in the Gap increases the free energy and decreases m, transforming the narrow filament into the current island.
We suppose that the radial flux is organized by small current filaments fluctuating in a radial direction. These filaments, accumulating free energy, appear due to deformation of the pressure profile in the Gap, forming current islands. Depending on the value of the stored free energy, they will have larger or smaller values of m and n, preserving the value / q m n = .
This has been well demonstrated in beautiful KSTAR experiments ( Figure 13) [25]. In these experiments, the fine Te (and so, the current density) structure was registered from the image of the second EC harmonic emission. This gave the fine 2D temperature profile in some core plasma region with 1 q ≤ . Note that the current, flowing inside this region, forms a spiral with / 1/1 m n = (kink mode), and the position of the EC resonance relative to the current center can be different. The geometry is complicated, and at various ECR heating positions relative to the current center, we may obtain q = 1 with m/n from 1/1 up to 3/3. Similar results have been obtained in KSTAR experiments for small islands generated at the top of the external barrier (H-mode), the so-called edge localized modes (ELMs). Their amplitude and m number depend on the radial heat flux [26].
Energy overload of the islands cannot reduce m below the minimum value for a given surface (m = 1 in the KSTAR case shown in Figure 13). Further energy growth should lead to internal disruption, which is the mechanism that expels the energy through the barrier (see Figures 11 and 12 as well as the KSTAR paper [25]). In this case, after the disruption, Γ(r) decreases in the region with q ≤ 1, and m/n changes from 1/1 to 8/8. Internal disruption is a new mechanism for maintaining the selforganized pressure profile. However, it is too rough and leads to the onset of large-scale MHD instabilities and global plasma disruption. Thus, the ITB, which initially leads to improvement in plasma energy confinement, is a precursor to its collapse. Figure 13. KSTAR experiment [25]. Islands with m/n = 3/3, 2/2, and 1/1 may exist interior to the region q = 1. Reproduced with permission IAEA. Figure 13. KSTAR experiment [25]. Islands with m/n = 3/3, 2/2, and 1/1 may exist interior to the region q = 1. Reproduced with permission IAEA.
Similar results have been obtained in KSTAR experiments for small islands generated at the top of the external barrier (H-mode), the so-called edge localized modes (ELMs). Their amplitude and m number depend on the radial heat flux [26].
Energy overload of the islands cannot reduce m below the minimum value for a given surface (m = 1 in the KSTAR case shown in Figure 13). Further energy growth should lead to internal disruption, which is the mechanism that expels the energy through the barrier (see Figures 11 and 12 as well as the KSTAR paper [25]). In this case, after the disruption, Γ(r) decreases in the region with q ≤ 1, and m/n changes from 1/1 to 8/8. Internal disruption is a new mechanism for maintaining the self-organized pressure profile. However, it is too rough and leads to the onset of large-scale MHD instabilities and global plasma disruption. Thus, the ITB, which initially leads to improvement in plasma energy confinement, is a precursor to its collapse.

Discussion
Plasma in the tokamak magnetic field is an open nonequilibrium system with many degrees of freedom. Such a system must be turbulent and self-organized, and it is not surprising that this is observed in experiments. The plasma system is well described from the position of nonequilibrium thermodynamics, giving the opportunity to move away from scalings and fitting models in fusion reactor design.
Following Dyabilin and Razumova [8,9], the energy balance equation was obtained, and the necessary transport coefficients were empirically found, which were determined only by external parameters and independent of the internal plasma parameters, such as temperature, density, effective charge, and so forth.
Dyabilin's equation allows us to calculate any regime without an ITB for any tokamak, including the fusion reactor; experiments show that the physics processes are the same for all tokamaks, from the small RTP (R = 0.5 m, a = 0.16 m) to the large JET (R = 3 m, a = 0.9 m). We qualitatively explained the processes with internal barriers, but for quantitative calculations, the correct theory, such as that created by Dyabilin for a self-consistent pressure profile, is desirable. In [8,9], the functional was minimized only for diamagnetic currents because the experiment showed that the self-consistent pressure profile was independent of the longitudinal current density profile. To solve the problem with barriers, it is necessary to do similar work for longitudinal currents. However, the calculation of regimes with barriers is more difficult, since it requires knowledge of the longitudinal current density profile, which can be different and vary over time.
As shown above, in a fusion reactor, it is necessary to radiate a significant part of the generated power but in the plasma volume (radiation in the divertor volume cannot help). How can this be done without creating a significant concentration of radiating ions in the hot core, which can lead to unacceptable intensity of bremsstrahlung? We have to analyze this problem.

Conclusions
In this paper, on the basis of the thermodynamic approach, the phenomena of improved confinement at impurity influx and saturation of confinement dependence on the density in tokamaks were explained. A hypothesis was proposed to explain the formation of internal transport barriers and magnetic islands.