Determination of the Dynamic ITER Energy Confinement Time Scalings

We derive the differential equation, which is satisfied by the ITER scalings for the dynamic energy confinement time. We show that this differential equation can also be obtained from the differential equation for the energy confinement time, derived from the energy balance equation, when the plasma is near the steady state. We find that the values of the scaling parameters are linked to the second derivative of the power loss, estimated at the steady state. As an example of an application, the solution of the differential equation for the energy confinement time is compared with the profile obtained by solving numerically the balance equations (closed by a transport model) for a concrete Tokamak-plasma.


I. INTRODUCTION
Global scaling expressions for the energy confinement time, τ E , or the stored energy, W , are powerful tools for predicting the confinement performance of burning plasmas [1], [2], [3]. The fusion performance of ITER is predicted using three different techniques: statistical analysis of the global energy confinement data in the parameters (simple (multivariate) linear regression tools can be used to determine the parameters from a set of data) [4], [5], a dimensionless scaling analysis, based on dimensionless physics parameters [6], [5], [7], and theory -based on transport models and modelling the plasma profiles [8], [9] and [10]. Although the three methods give overlapping predictions for the performance of ITER, the confidence interval of all of the techniques is still quite wide [11]. The Confinement Database and Modelling Expert Group recommended for ITER design the so-called IT ERH − 98P (y, 2) confinement scaling [5], [12]: This expression is valid for the ELMy H-mode thermal energy confinement time. In Eq. (1), the parameters are the plasma current I p , the major radius R, the inverse aspect ratio = a/R (with a denoting the minor radius of the Tokamak), the elongation κ, the toroidal magnetic field (at the major radius R) B 0φ , the central line averaged electron densityn e , the loss power P , and the ion mass number M , respectively. The point prediction for the thermal energy confinement time in ITER is τ E = 3.6 sec. The 2 log-linear interval was determined to be 20%. By recent analyzing the enlarged IT ERH.DB3 dataset, the practical reliability of the IT ERH − 98(y, 2) scaling was confirmed and 2 log-linear interval was reduced to 14% [13].
In this work, we shall show that Eq. (1) [or Eq. (2)] satisfies a nonlinear differential equation of second order in time. This equation can also be derived from the energy balance equation. This will allow a linking of the scaling coefficients with the (measurable) second time derivatives of the heat power loss, estimated at the stationary state. These tasks will be accomplished in Sections (II) and (III) respectively. As an example of an application, in the Section (IV), we compare the solution obtained by solving the confinement time differential equation with the numerical simulations obtained using the code JETTO [23], for the specific case of IGNITOR-plasmas. Concluding remarks can be found in Section (V).

II. DIFFERENTIAL EQUATION SATISFIED BY THE ITER SCALINGS
The expression for the energy confinement time, obtained by scaling laws, raises several questions. Firstly, Eq. (1) applies quite well to a large number of Tokamaks (ASDEX, JET, DIII-D, ALCATOR C-Mod, COMPASS, etc.) and it is currently used for predicting the energy confinement time for Tokamaks, which are presently in construction (ITER) or will be constructed in the future (DEMO). Hence, the first objective of this work is to understand the main reason for such a "universal" validity. Secondly, it is legitimate to ask "where does this expression ?". More concretely, "Is it possible to determine the (minimal) differential equation which is satisfied by expression (1) ? ". In case of a positive answer, "Is it possible to re-obtain this (minimal) differential equation from the balance equations and, in particular, from the energy balance equations ?". Finally, "How can we estimate the values of the scaling coefficients α i ?". In this Section, we shall determine the (minimal) differential equation satisfied by Eq. (1). In the next Section we shall prove that, near the stationary state, this differential equation can be re-obtained from the energy balance equation.
The equations of one-dimensional plasma dynamics, in toroidal geometry, assuming the validity of the standard model, can be brought into the form (see, for example, [20]) with r and q(r) denoting the radial coordinate and the safety factor, respectively. p, n e , T e and Z are the total plasma pressure, the electron density, the electron temperature and the ion charge number, respectively. Here, < · · · > denotes the surface-average operation.
< q ζ > and < γ e r > are the averaged radial heat flux of species ζ (ζ = e for electrons and ζ = i for ions) and the averaged electron flux, respectively. c and E 0 are light speed and the external electric field, respectively, and S gain−loss is the source term, i.e. the loss and energy with where the "dot" over the variables stands for the (total) time derivative (d/dt).
The energy confinement time is defined as From definition (6), we findτ Note that the stationary state is reached when P Q = P tot . Hence, at the steady state At the steady state, we find where W 0 e and P 0 tot indicate the values of W e and P tot , estimated at the steady state, respectively.
Eqs (1) and (2) can be re-written in the generic form: where X 1 , X 2 , · · · are a positive and independent system of variables X i , and α i the scaling parameters, respectively. Note that C is a (dimensional) constant satisfying the condition Unless stated otherwise, in the sequel we shall adopt the summation convention on the repeated indexes. By taking the logarithm of Eq. (10) we find with y ≡ log τ E and ξ i ≡ log X i (with i = 1, · · · n). The first and the second derivatives of y, with respect to variable ξ i , read respectively In terms of variable τ E , instead of y, we get The differential equation with respect to time is easily obtained by tacking into account the By multiplying the second equation of Eqs (14) byξ iξj and by summing over indexes, we finally obtain the differential equation satisfied by the ITER scaling laws Eq. (16) should be solved with the initial conditions (9): with We have derived two differential equations for the time derivatives of τ , the first equation of Eq. (15) which is first order and also Eq. (16) which is second order. It may appear hopeless to solve these equations, as they depend on α iξi (t) and χ(t) = α iξi (t) respectively, which in turn depend on the full dynamics of the system.
The critical fact which makes our approach useful is that the second time derivatives of the logarithm of X i are generally weakly dependent on time. As a result, one may approximate χ(t) to be a constant, χ 0 . In this sense, all of the dependence on the machine is reduced to just a number, which can be determined. The evolution of τ E can then be obtained uniquely by integrating Eq. (16) with the initial conditions (9). Such an approach would not work for the first equation of Eq. (15) as α iξi (t) depends strongly on time, indeed it vanishes at the initial stationary state and then becomes nonzero as the state evolves.
It is not difficult to check that Eq. (17) is the "minimal" differential equation, in the sense that Eq. (17) admits one, and only one, solution given by Eq. (10). For this, it is sufficient to note that Eq. (17) can be brought into the form By integrating this equation twice, from t 0 to t, and taking into account thatτ 0 Recalling that ξ i ≡ log X i , we find where Eq. (11) has been taken into account. In brief, to the expression Eq. (10) we can associate the differential equation (17). This equation admits one, and only one, solution given by expression (10).
We conclude this Section by mentioning that in all the cases examined by the authors, χ(t) is very well approximated (numerically) by a linear function in time where the coefficients c i (with i = 1, 2) are estimated at the steady state.

III. DIFFERENTIAL EQUATION FOR THE ENERGY CONFINEMENT TIME
The aim of this Section is to obtain the differential equation for the energy confinement time from the balance equations. In analogy with Eq. (17), the coefficients of this differential equation should be expressed only in terms of the internal energy W e and the total power P tot . To this end, let us reconsider the energy balance equation Eq. (4) and the definition of the energy confinement time, Eq. (6). Taking the derivative of Eq. (7) with respect to time, after a little algebra, we get with Note that the dimensions of f (t) and g(t) are [t] −2 and [t] −1 , respectively. Finally, the differential equation for the energy confinement time reads Note that Eqs (24) may be re-written in the more convenient form showing that the differential equation for the energy confinement time may be expressed as two quasi-decoupled differential equations of first order in time derivative. The general solution of Eqs (25) may be brought into the form By taking into account that f (t) = − n i=1 α iξi (with ξ i = log X i ), solution (26) generalizes the ITER scaling laws out of the steady state, reducing to Eq. (1) [or Eq. (2)] close to the stationary state.

• Differential Equation for the Energy Confinement Time Near the Steady State
The term P tot is specified as follows where P α (T ) is the alpha power, P b (T ) is the power radiation loss (Bremsstrahlung) and P aux is the external heating power density supplied to the system (e.g. ohmic heating power or external RF), respectively. The alpha power and the Bremsstrahlung power loss depend explicitly on the temperature of the plasma. The auxiliary heating power is operational during both the transient and steady states. This is the dominant source of external heating power, and it is assumed to be deposited in the plasma with a known profile, independent of p and T . Hence, P Aux = P Aux (r, t). The time derivative of P tot readṡ At the steady stateṪ (t 0 ) = 0 andṖ Q (t 0 ) =Ṗ Aux (t 0 ) = 0. Consequently, from the energy balance equation we find that alsoẄ e (t 0 ) = 0. By taking into account Eqs (8) and (23), we get g(t) → 0 as the system approaches the steady state. Hence, near the stationary state, where Eq. provides the desired relation between the exponent coefficients α i and the macroscopic quantities P tot and W e . If we have n free exponent coefficients α i , we can set the following n Eq. (31) link the exponent coefficients with variables which, at least in principle, are under the control of the experimental physicist.
• Analysis in the "Physics" Variables Looking at Eqs (1) and (2) we realize that several independent variables are timeindependent (e.g., major and minor radii, elongation etc. where, now, the independent variables X α i i (t)/X α i i (t 0 ) are dimensionless. Note that in this case variables ξ i are defined as ξ i = log(X i /X 0 i ) (no summation convention over the repeated indexes). Of course, this operation reduces the number of independent variables. However, this number may be reduced further if, instead of "engineering variables", the confinement time is expressed in terms of "physics" parameters such as ρ (normalized Larmor radius), β (normalized pressure), ν (collisionality), etc. Indeed, according to the observation of Kadomtsev, the transport in the plasma core should be fundamentally governed by three physical dimensionless plasma parameters ρ , β and ν [21]. In this respect, an interesting paper is Ref. [22]. In [22] the authors show that, due to the Kadomtsev constraint, the final expression for the ELMy H-mode thermal confinement time has only one free exponent coefficient, according to the law: where Eq. (32) has been taken into account. We find

IV. COMPARISON WITH THE NUMERICAL SIMULATION OF THE BAL-ANCE EQUATIONS FOR AN L-MODE TOKAMAK-PLASMA
As an example application, we consider in this Section the case of one of the simplest L-mode Tokamak-plasma where the evolution of the energy confinement time has been estimated by solving numerically the balance equations, completed with a transport model.
In [23] we find the profile of τ E against time for Ignitor-plasma. The numerical solution has been obtained by using the code JETTO. To compare this profile with the numerical solution of Eq. (17), we should firstly estimate t 0 , τ 0 E , c 1 =P Q (t 0 )/P tot (t 0 ) and c 2 = − d dt (P Q (t)/P tot (t)) | t=t 0 [see Eqs (21) and (22)]. In [24], we have estimated the values of these parameters for Ignitor subject to ICRH power (i.e., P Aux = P ICRH ). The scenario is considered where IGNITOR is led to operate in a slightly sub-critical regime by adding a small fraction of 3 He to the nominal 50 − 50 Deuterium-Tritium mixture. The difference between power lost and alpha heating is compensated by an additional ICRH power equal to 1.46 MW, which should be able to increase the global plasma temperature via collisions between 3 He minority and the background D − T ions. The analytical expression for the ICRH power profiles inside the plasma has been deduced by fitting the numerical results giving an expression for P Aux = P ICRH (r), which is essentially independent of the bulk temperature. Denoting the ICRH power-density as P ICRH , we have with P 0ICRH = 6.59126 10 −6 M W/m 3 ,α 2 = 15.3478 and ∆ = 0.0477032, respectively.
The value of τ 0 E has been estimated by the expression [24] τ 0 E = 12n e T E α n 2 e < σv > D−T −C B n 2 e T 1/2 + 4P ICRH  in case B and the highest in case D, as shown in Fig. 8.
The total energy confinement time evolution is plotted in Fig. 9 for case C together with the ITER97L scaling evaluated by using P input , P net and a further assessment with a reduced power P red = P net P RadTot . Note that, until ignition, the energy confinement time remains below the scaling expression evaluated by using P red . Figures 10-13 display the normalized pressure profiles at fixed times along the simulations and at ignition. The plots again highlight the di↵erent evolutions; but at the ignition time, specific for each case, the plasma pressure assumes a consistent shape. This 'format' does not depend on the path followed for obtaining ignition but is biased by the time spent. In fact as the ignition time becomes longer, the pressure profile becomes more peaked, as shown in cases B, C and D compared with case A. It should be noted that density and temperature profiles, separately, do not exhibit the same characterization. The simple power balance between the alpha source (assuming P ↵ / n 2 T 2 ) and the plasma losses suggests that, given a ⌧ E value, the leading parameter turns out to be the plasma pressure. The consistency of the pressure profile is stressed in Fig. 14, where all the final shapes are drawn. The privileged profile can be fitted by the following analytic formulas: P tot − P RadT ot , whereas in our work we use P tot , which includes the Bremsstrahlung radiation loss. This may explain the little difference between the numerical [23] and the analytical slopes.

V. CONCLUSIONS
A large database on plasma energy confinement in Tokamaks can be summarized in single empirical value of τ E , referred to as ITER-scalings. These expressions are "Universal", in the sense that they apply to a large number of Tokamaks. Scalings are expressed in terms of product of powers of independent variables [see Eq. (10)] and correspond to the L-mode as well as the H-mode confinements. The recommended scaling for ITER operation remains the IP B98 scaling law, while this issue is further investigated. In this work we have shown that the ITER scalings satisfy a general non-linear differential equation of second order in time.
For given initial conditions, this differential equation admits one, and only one, solution, which coincides exactly with the ITER scalings. In parallel, we determined the differential equation for energy confinement time from the energy balance equation. This differential equation contains a nonlinear extra term, which behaves as ∼ τ EτE . However, when plasma is near the steady state, this differential equation reduces to the one admitting the ITER scaling as solution. We have also seen that the scaling coefficients may be linked to the variables which, at least in principle, are under the control of the experimental physicist.
The validity of our approach has been tested by analyzing a concrete example of Tokamakplasma where the profile of the energy confinement time has been previously determined by solving the balance equations (with the auxilium of a transport model). The solution of the differential equation for the ITER scaling is in a fairly agreement with the numerical finding.

VI. ACKNOLWEDGMENTS
One of us (G.Sonnino) is very grateful to Dr Fulvio Zonca and Dr Alessandro Cardinali, from the ENEA for EUROfusion in Frascati (Italy), for the useful discussions and suggestions. JE is supported by NSFC MianShang grant 11375201.