Received 19 June 2016; accepted 8 August 2016; published 11 August 2016
Global scaling expressions for the energy confinement time, , or the stored energy, W, are powerful tools for predicting the confinement performance of burning plasmas  -  . 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)   , a dimensionless scaling analysis, based on dimensionless physics parameters  -  , and theory- based on transport models and modelling the plasma profiles  -  . Although the three methods give overlapping predictions for the performance of ITER, the confidence interval of all of the techniques is still quite wide  . The Confinement Database and Modelling Expert Group recommended for ITER design the so-called confinement scaling   :
Here, the parameters are the plasma current, the major radius R, the inverse aspect ratio (with a denoting the minor radius of the Tokamak), the elongation, the toroidal magnetic field (at the major radius R), the central line averaged electron density, the loss power P, and the ion mass number M, respectively. The expression (1) is valid for the ELMy H-mode thermal energy confinement time. The 2log-linear interval was determined to be 20%. By recent analyzing the enlarged dataset, the practical reliability of the scaling was confirmed and 2log-linear interval was reduced to 14%  . Tables showing some of the most generally used sets of scaling parameters for the ELMy H-mode and L-mode can be found in Refs   -  .
For stellarators, a similar scaling has been obtained  
where is the rotational transform (or the field line pitch).
The confinement time is defined as
where, and are the internal energy, the power loss and the power source, respectively. From Equation (3) results that when the tokamak is not in the steady state the quantity is a time dependent quantity. Hence, , given by Equations (1) and (2), is viewed as a time-dependent variable, which depends on a collection of variables dependent on time (e.g., , P, etc.). The value of at the steady state condition, attained at some time moment, corresponds to the numerical value provided by the database. For example, the point prediction for the thermal energy confinement time in ITER is.
The main objective of this work is to estimate the energy confinement time, close to the steady state. at the steady state condition is calculated by using the expression
where and are obtained by solving the stationary balance equations. An example of calculation can be found in Ref.  . To estimate the dynamic confinement time we should solve the evolutive balance equations. However, this is a very complex task. An alternative strategy (which is the one that we shall adopt here) consists in deriving the time differential equation for the energy confinement time, with, estimated by using Equation (4), playing the role of the initial condition. We show that is the solution of a nonlinear differential equation of second order in time, obtained by combining Equation (3) with the (dynamic) balance equations. The critical fact which makes our approach useful is that in the vicinity of the stationary state, this differential equation depends only on one coefficient which varies very slowing in time
where, and is a numerical coefficient estimated at the steady state. Hence, at the “leading order”, all of the dependence on the machine is reduced to just a number, , which can be determined. This is the real advantage of this approach. As an example of calculation, we have considered the simplest case of IGNITOR-plasmas. In this case, we solved the time differential equation for where the parameters (i.e., the initial condition as well as the coefficient appearing in the differential equation) have been estimated at the steady state. The solution of this equation is in agreement with the one obtained by solving numerically the dynamic balance equations, with the aid of a transport model  .
In this work, we shall also justify the dynamic scaling laws, like
where C is a constant and M is the effective mass, respectively (note that when the plasma is a mixture, due to the dependence of particle transport properties on particle mass and charge, M is also time dependent). In particular, we shall prove that the dynamic expression for the energy confinement time, like Equation (6), is solution of the differential equation for, which can be obtained by combining Equation (3) with the energy balance equation.
The paper is organized as follows. In Section (2), we show that Equation (6) satisfies a nonlinear differential equation of the second order in time, tacking into account the (experimentally established) slow variation in time of the coefficient entering in this equation. Successively, we show that this equation can also be derived from the energy balance equation, combined with definition (3). This will allow a linking of the scaling coefficients with the (measurable) second time derivatives of the heat power loss, which at the leading order may also be estimated at the stationary state. These tasks will be accomplished in the Section (3). As an example of an application, in the Section (4), we compare the solution obtained by solving the differential equation for the energy confinement time with the numerical simulations obtained using the code JETTO  , for the specific case of IGNITOR-plasmas. Concluding remarks can be found in Section (5).
2. Differential Equation Satisfied by the ITER Scalings
The expression for the energy confinement time, obtained by scaling laws, raises several questions. Firstly, Equation (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 originate from ?”. More concretely, “Is it possible to determine the (minimal) differential equation which is satisfied by expression (6)?”. 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?”. In this Section, we shall determine the (minimal) differential equation satisfied by Equation (6). 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,  )
with r and denoting the radial coordinate and the safety factor, respectively. p, , 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. and are the averaged radial heat flux of species (for electrons and for ions) and the averaged electron flux, respectively. c and are light speed and the external electric field, respectively, and is the source term, i.e. the loss and energy gain. Equation (7) must be completed with the transport equations, i.e. with the thermodynamic flux- force relations, in order to close the plasma dynamical equations. The power balance equation is now derived as follows. Equation (7) is integrated over the volume of the plasma and then divided by the plasma volume V. We obtain
where the “dot” over the variables stands for the (total) time derivative ().
The energy confinement time is defined as
From definition (10), we find
Note that the stationary state is reached when. Hence, at the steady state (corresponding to) we have
At the steady state, we find
where and indicate the values of and, estimated at the steady state, respectively.
Equation (6) may be re-written in the generic form:
where are a positive and independent system of variables, and the scaling parameters, respectively. For simplicity, we firstly suppose that in Equation (14) all the variables are time-dependent. The case whereby is a collection of variables dependent on time, as well as variables not-dependent on time, will be treated in the following sub-Section Analysis in the Physics Variables. 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 Equation (14) we find
with and (with). The first and the second derivatives of y, with respect to variable, read respectively
In terms of variable, instead of y, we get
The differential equation with respect to time is easily obtained by tacking into account the identities
By multiplying the second equation of Equation (18) by and by summing over indexes, we finally obtain the differential equation satisfied by the ITER scaling laws
Equation (20) should be solved with the initial conditions (13):
with. We have derived two differential equations for the time derivatives of, the first
equation of Equation (19) which is first order and also Equation (20) which is second order. It may appear hopeless to solve these equations, as they depend on and 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 are generally weakly dependent on time. As a result, one may approximate to be a constant,. In this sense, all of the dependence on the machine is reduced to just a number, which can be determined. The evolution of can then be obtained uniquely by integrating Equation (20) with the initial conditions (13). Such an approach would not work for the first equation of Equation (19) as 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 the nonlinear Equation (21) is the “minimal” differential equation, in the sense that Equation (21) admits one, and only one, solution (i.e., the nonlinear differential Equation (21) does not generate additional solutions).
It may appear hopeless to solve Equation (21), as it depends on the coefficient, 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 are generally weakly time-dependent. In all the cases examined by the authors, is very well approximated (numerically) by a linear function in time
Hence, all of the dependence on the machine is reduced to just a number, , which can be estimated at the steady state.
3. 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 Equation (21), the coefficients of this differential equation should be expressed only in terms of the internal energy and the total power. To this end, let us reconsider the energy balance equation Equation (8) and the definition of the energy confinement time, Equation (10). Taking the derivative of Equation (11) with respect to time, after a little algebra, we get
Note that the dimensions of and are and, respectively. Finally, the differential equation for the energy confinement time reads
We might object that the previous equation has the same degree of difficulty as the initial expression, Equation (10). However, as we shall see more in detail in the next Subsection, the coefficients and possess special properties: close to the steady state tends to vanish and is a function varying very slowing in time. So, at the leading order, and may be estimated at the stationary state [see Equation (22) and the discussion after Equation (21)]. This is the real advantage of Equation (25) with respect to Equation (10): Equation (25) allows determining the dynamic behaviour of the energy confinement time when the system is close to the steady state, solely by the knowledge of one coefficient estimated at the stationary state. Moreover, from the previous Section we know that this equation admits one (and only one) solution corresponding to the ITER scalings. A concrete application of Equation (25) can be found in the Section (4). Note that Equation (25) 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 Equations (26) may be brought into the form
By taking into account that (with), solution (27) generalizes the ITER scaling
laws out of the steady state, reducing to Equation (14) close to the stationary state. Equation (27) shows that close to the steady state, the leading contribution to the mathematical expression for the energy confinement time is provided by the power laws. However, when we deviate from the steady state, supplementary con- tributions, which are different from the power ones, may modify the mathematical form of the power laws significantly. Generally, for ITER, these contributions tend to lower the numerical value of the energy con- finement time. This can be easily checked by setting in Equation (27), with and denoting a very small parameter and a positive constant having dimension, respectively. By developing expression (27) up to the first order in, we find the power laws at the leading order, corrected by a (small) negative expression at the first order in.
・ Differential Equation for the Energy Confinement Time Near the Steady State
The term is specified as follows
where is the alpha power, is the power radiation loss (Bremsstrahlung) and 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,. The time derivative of reads
At the steady state and. Consequently, from the energy balance equation we find that also. By taking into account Equations (12) and (24), we get as the system approaches the steady state. Hence, near the stationary state, we find
where Equation (22) has been used. As shown in the Section (2), Equation (30) admits one (and only one) solution, corresponding to the ITER scalings Equation (6). Note that Equation (31) provides the desired relation between the exponent coefficients and the macroscopic quantities and. If we have n free ex- ponent coefficients, we can set the following n relations
Equation (32) link the exponent coefficients with variables which, at least in principle, are under the control of the experimental physicist.
・ Analysis in the “Physics” Variables
As mentioned, Equations (1) and (2) are composed by several variables independent of time (e.g., major and minor radii, elongation etc.). In this case, it is more convenient to express the energy confinement time only in terms of the time-dependent variables. Let us suppose that m variables are time-dependent and the remaining not. In this case, the energy confinement time takes the form [see Equation (15)]
where, now, the independent variables are dimensionless. Note that in this case variables are defined as (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  . In this respect, an interesting paper is Ref.  . In  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:
with denoting the density of the power loss (i.e.,). With the choice, in “physics” variables, scaling (34) goes as (i.e. a gyro-Bohm-like scaling), and. This choice may be tested by using Equation (32) which, in this particular case, reads
where Equation (33) has been taken into account. We find
4. Comparison with the Numerical Simulation of the Balance 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  we find the profile of 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 Equation (21), we should firstly estimate, and [see Equation (22)
and (24)]. In  , we have estimated the values of these parameters for Ignitor subject to ICRH power (i.e.,). The scenario is considered where IGNITOR is led to operate in a slightly sub-critical regime by adding a small fraction of 3He 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 3He minority and the background ions. The analytical expression for the ICRH power profiles inside the plasma has been deduced by fitting the numerical results giving an expression for, which is essentially independent of the bulk temperature. Denoting the ICRH power-density as, we have
with, and, respectively.
The value of has been estimated by the expression 
with and denoting the energy at which the alpha particles are created (3.5 MeV), and the Brems- strahlung constant, respectively. is the reaction cross section giving a measure of the probability of a fusion reaction as a function of the relative velocity of the two reactant nuclei. provides an average over the dis- tributions of the product of cross section and velocity v. In the core of the plasma we found  , and. Figure 1 reports on the energy confinement time, , against time for Ignitor-plasmas in the above mentioned conditions. The profiles have been obtained by solving (with the code JETTO) the balance equations and refer to the ITER scalings ITER97L (full dots), (open dots) and ITER97L  . Figure 2 shows the solutions of the differential equation for the ITER scalings, Equation (21), at the three values of ():(ITER97L-blue line), (-green line) and = (ITER97L(P_red)-brown line).
Note that in  the authors evaluate the ITER scalings by using the reduced power, whereas in our work we use, which includes the Bremsstrahlung radiation loss. This may explain the little difference between the numerical  and the analytical slopes.
A large database on plasma energy confinement in Tokamaks can be summarized in single empirical value of, referred to as the 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 Equation (14)] and correspond to the L-mode as well as the H-mode confinements. The recommended scaling for ITER operation remains the IPB98 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. The value provided by the database for ITER scaling laws, coincides with, estimated by Equation (4), with and evaluated by solving the stationary balance equations. To estimate the dynamic confinement
Figure 1. Solutions of Equation (21) at the three values of. Blue line: (0.35 sec, 0.43 sec) (ITER97L), green line: () and Brown line: = (ITER97L(P_red)).
Figure 2. Energy confinement time evolution estimated in  by solving with JETTO the balance equations (completed with a transport model): ITER97L scaling (full dots-blue line), scaling (open dots-green line) and ITER97L(P_red) scaling (brown line).
time, we determined the differential equation for by combining the energy balance equation with definition (3). We found Equations (25). We have solved this equation by taking into account that, in vicinity of the steady state, the coefficient tends to vanish and, at the leading order, is (almost) a constant independent of time, which may be evaluated at the stationary state. This is the real advantage of the proposed approach: close to the steady state, the differential equation for the energy confinement time reduces to
where “at the leading order” is a numerical constant, which may be estimated at the stationary state. As a result, one may approximate to be a constant, or better, by a linear function. In this sense, all of the dependence on the machine is reduced to just a number, , which can be estimated at the steady state. Far from the stationary state the differential equation for contains a nonlinear extra term, which behaves as. This extra term tends to modify the mathematical form of the power laws. For ITER, the main effect of this nonlinear extra term is to lower the numerical value of the energy confinement time. The general solution is given by Equation (27), which reduces to the one admitting the ITER scaling power laws as the system approaches the steady state. 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 Tokamak-plasma 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.
Giorgio Sonnino is very grateful to Dr Philippe Peeters for his useful suggestions. Jarah Evslin is supported by NSFC MianShang grant 11375201. The reproduction of Figure 2, reported in Ref.  , has been authorized by the review Nuclear Fusion (IoP).