Hysteresis and fast timescales in transport relations of toroidal plasmas

This article assesses current understanding of hysteresis in transport relations, and its impact on the field. The rapid changes of fluxes compared to slow changes of plasma parameters are overviewed for both core and edge plasmas. The modulation ECH experiment is explained, in which the heating power cycles on-and-off periodically, revealing hysteresis and fast changes in the gradient–flux relation. The key finding is that hystereses were observed simultaneously in both the the gradient–flux and gradient–fluctuation relations. Hysteresis with rapid timescale exists in the channels of energy, electron and impurity densities, and plausibly in momentum. Advanced methods of data analysis are explained. Transport hysteresis can be studied by observing the higher harmonics of temperature perturbation δTm in heating modulation experiments. The hysteresis introduces the term δTm, which depends on the harmonic number m in an algebraic manner (not exponential decay). Next, the causes of hysteresis and its fast timescale are discussed. The nonlocal-in-space coupling works here, but does not suffice. One mechanism for ‘the heating heats turbulence’ is that the external source S in phase space for heating has its fluctuation in turbulent plasma. This coupling can induce the direct input of heating power into fluctuations. The height of the jump in transport hysteresis is smaller for heavier hydrogen isotopes, and could be one of the origins of isotope effects on confinement. Finally, the impacts of transport hysteresis on the control system are assessed. Control systems must be designed so as to protect the system from sudden plasma loss.

This article assesses current understanding of hysteresis in transport relations, and its impact on the field. The rapid changes of fluxes compared to slow changes of plasma parameters are overviewed for both core and edge plasmas. The modulation ECH experiment is explained, in which the heating power cycles on-and-off periodically, revealing hysteresis and fast changes in the gradient-flux relation. The key finding is that hystereses were observed simultaneously in both the the gradient-flux and gradient-fluctuation relations. Hysteresis with rapid timescale exists in the channels of energy, electron and impurity densities, and plausibly in momentum. Advanced methods of data analysis are explained. Transport hysteresis can be studied by observing the higher harmonics of temperature perturbation δT m in heating modulation experiments. The hysteresis introduces the term δT m , which depends on the harmonic number m in an algebraic manner (not exponential decay). Next, the causes of hysteresis and its fast timescale are discussed. The nonlocal-in-space coupling works here, but does not suffice. One mechanism for 'the heating heats turbulence' is that the external source S in phase space for heating has its fluctuation in turbulent plasma. This coupling can induce the direct input of heating power into fluctuations. The height of the jump in transport hysteresis is smaller for heavier hydrogen isotopes, and could be one of the origins of isotope effects on Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction and background of the problem
In the history of the study of confinement in toroidal plasmas, the approach of analyzing the response to periodic modulation of heating power has been applied routinely (for a review see e.g. [1]). A conventional approach has been to use the 'diffusive model', in which one assumes that the perturbed heat flux δq is proportional to the changes of local plasma parameters and their gradients. When the changes are dominant in the electron temperature and its gradient, the local diffusive model employs coefficients for thermal conductivity and heat pinch. In the analysis of modulation electron cyclotron heating (MECH) experiments, the 'heat pulse' thermal diffusivity, χ HP , has been deduced, but the discrepancy between this and the thermal diffusivity in the stationary power balance, χ PB , has long been recognized (see e.g. [2]). The physics process that causes the difference between χ HP and χ PB has been a mystery.
The understanding of the discrepancy between χ HP and χ PB has been advanced recently. There have been many experimental reports (e.g. on W7-AS [3], on DIII-D [4]) that the gradient-flux relation has a hysteresis (depending on switching on/off of ECH heating), and the heat flux may directly depend on the heating power (figure 1 [4,5]). This view shows a contrast to the conventional 'diffusive model' of the transport, in which the heat flux is expressed in terms of the local plasma parameters and their gradients. If the hysteresis in the gradient-flux relation exists, it is not surprising that the difference between χ HP and χ PB is quite large. An essential question in the pioneering work (such as [3,4]) was whether the profile of absorbed power by modulated heating was correctly predicted by theory or not. Measurement of the fluctuation intensity in the modulated heating has been performed on LHD [6]. The new, independent data of fluctuation intensity was found to change against the local temperature gradient with hysteresis on LHD [6,7]. Thus, the hysteresis is not an artifact due to the error in evaluating the power absorption profile. In such a situation, which is fully analyzed in [6], the ECH power of 2 MW is periodically applied on the plasma, which is sustained by the constant NBI heating power of 3 MW. The peak-to-peak modulation power, ∆P ECRH , is close to half the time-averaged heating power P total . In this case, it is reported that at the radius of r/a = 0.66 the peak-topeak change of heat flux per particle by MECH, ∆q e , amounts to ∼ 12 keVm s −1 , which is about one-half of the time-averaged heat flux per particle q e,total , i.e. ∆q e / q e,total ∼ 1/2 (thus is the ratio of ∆P ECRH / P total ). As is demonstrated in [6], the jump of the heat flux at the onset of ECH, δq jump , is close to 5 keVm s −1 . Thus, the jump of heat flux at on/off of heating is about 5/12 of the peak-to-peak change of heat flux, δq jump ∼ (5/12)∆q e . It was reported that when the modulation power of ECH increases, the ratio δq jump /∆q e increases. Therefore, the hysteresis in the gradient-flux relation is unambiguous, and plays essential roles in the dynamic response. Developments in identifying transport hysteresis have been assessed in the literature [5,8]. If the hysteresis actually exists, and the heat flux directly depends on the heating power, a discrepancy between 'power-balance' and 'heat-pulse' conductivities [1] arises. Transport hysteresis can also explain the disparity of heat pulse between inward-and outward-propagations [9], which cannot be explained by a local diffusive model. In addition, this hysteresis can have a profound impact on our evaluation of the dynamical response of burning plasmas. The understanding of these observations is necessary for the accurate prediction of burning plasma, and in designing the temporal control system to prevent dynamic changes of core plasma.

Modulation ECH (MECH) experiment
The key experimental discovery is the identification of the difference between the time scale of response of heat flux and that of global parameters after the sudden change of heating power. The modulation ECH (MECH) experiment, in which the heating power cycles on-and-off following the step function has revealed fast changes in the gradient-flux relation. The switching time in the modulational EC heating is shorter than one millisecond. The global plasma parameters change in the order of energy confinement time (a few 10 ms to a few 100 ms). The change of the heat flux was found to occur with the former short time, not with those for global parameters [6,7]. Thus, the heat flux cannot be expressed by a function of local plasma parameters. The decisive progress is that both the hystereses-in the gradient-flux and gradientfluctuation relations-were observed simultaneously, as is illustrated in figure 2 [8]. The hysteresis in the gradient-flux relation is considered to indicate the influence of heating on transport. This conclusion is based on the comparison of time scales. Analyses of observations that can be interpreted as evidence of hysteresis have been undertaken on various experiments [9,10]. Discussions were made on observations on LHD, DIII-D, HL-2A, JFT-2M, JT-60U, KSTAR, TJ-II, and W7-AS.

Advanced methods of data analysis
The discrimination of two time scales in changes of heat flux and of plasma temperature is clarified by improving the statistical accuracy in observing dynamical changes. By introducing convolution analysis, in which the time is segmented by time periods t j < t < t j + 2π/Ω (t j = 2πj/Ω) and data in this segment, X(t) at t = t j + s, is relabeled as X(s; j) so that an average is taken over N periods as By this convolution method, the periodic change of temporal evolution of T(r, t) at the switching on and off of heating power is measured accurately. In the case of [6], the time resolution of T(r, t) at on-off of heating reaches sub-milliseconds, and the discrimination between the two time scales in flux and gradient of temperature becomes unambiguous. In terms of Fourier components, the higher harmonics in T(r, t), e.g. ninth or higher harmonics, can be observed accurately enough compared to the noise level [5]. This perspective has become available by the introduction of the convolution method to the analysis. Transport hysteresis is searched for most simply by observing the higher harmonics of temperature perturbation in the MECH experiment. We study the case where the gradientflux relation (in the radial region where the heating power is absent) has a 'jump' in the hysteresis, q jump , in addition to conventional diffusive terms, where q r is the radial heat flux, n is the plasma density and χ and V are conventional transport coefficients. Let us write the contribution of the jump to the energy balance equation as −∇ · q jump = nQH(t). Its dependence on H(t) models the fact that the jump in the hysteresis occurs in a short time, in comparison with the period of heating modulation. One has the equation for the perturbed temperature δT as We assume that Q is a smooth function in radius, for the sake of analytic transparency, so that |∆Q/Q| is much smaller than the wavenumber of the temperature perturbation. (This condition is satisfied across a wide region of the plasma column in the example of [6].) The special solution, which varies slowly in radius following Q, is given as [11] The homogeneous solution, which is used in the conventional analysis with the diffusive model without transport hysteresis (i.e. q jump = 0), obeys the dispersion relation for the waveform x is the distance from the reference radius. In using equation (3), geometrical effects in the cylinder or torus are often neglected. (The range of validity of this geometrical simplification is discussed in [9].) The wave number is predicted from equation (6a) as a function of the frequency ω as For the components with higher frequencies, (2m − 1) V 2 HP /4χ HP Ω, one has, from equation (6c),  [4], and reproduced from [5] with permission from JSPF) . for the (2m − 1)th harmonics. That is, for the homogeneous solution (the case where the heat flux is given by the diffusion model), the profile of temperature perturbation of the nth harmonics δT 2m−1 under the MECH equation (1) satisfies the relation One comes to the conclusion that, while the diffusion model of flux predicts the exponential decay of perturbation amplitude at higher harmonics, the jump in the gradient-flux relation with hysteresis induces algebraic decay with respect to the harmonic number. Comparing equations (5) and (7b), one sees that the radial profile of the higher harmonic is a key to investigating the properties of the gradient-flux relation. We here explain the comparison of the transport hysteresis model with experimental observations. Figure 3 summarizes experimental observations on LHD, TJ-II, KSTAR and DIII-D. The ECH frequency is chosen as fundamental heating (LHD) and second harmonic heating (TJ-II, KSTAR and DIII-D). The diagnostics view lines are on the low-field side in LHD and DIII-D, and on the high-field side in TJ-II and KSTAR. Thus the conclusion of this article is not specific to particular heating schemes or to limited diagnostic configurations.
(See [9] for further description of experiments.) The dependence of amplitude on the harmonic number has been tested on LHD and TJ-II (figures 3(a) and (b)). It has been shown that the amplitude of higher harmonics decays in radius with decay index similar to that of the fundamental mode. In figure 3, the prediction of the local diffusive model (i.e. with no jump in the hysteresis), equation (7b), is illustrated by red and blue dotted lines for 3rd and 5th harmonics. The red and blue dotted lines are plotted using the radial profile of the fundamental harmonics and equation (7b). Comparing the observed amplitude with these predictions, we see that the amplitude of higher harmonics has much weaker dependence on m than the prediction of equation (7b). The exponential dependence on the harmonic number, which should appear if the system obeys the local diffusive model (i.e. has no jump in the hysteresis), is not observed. We here note that the experimental observation for the case of LHD was reproduced by the transport model with jump in the heat flux [6]. A similar analysis has been performed on the MECH experiment on KSTAR and DIII-D [9]. Figures 3(c) and (d) are the cases of KSTAR and DIII-D respectively. In the propagation region of 0.3 < r/a < 0.5, the radial shapes of amplitudes of third and fifth harmonics are similar to that for the fundamental mode. In the conventional analysis, the slow radial decay of perturbation amplitudes of higher harmonics leads to the fitting to larger values of transport coefficients for higher harmonics. In the conventional analysis, in which the hysteresis is neglected, one uses the homogeneous solution to obtain the interpretation from equation (6a) as The characteristic dependence of 'effective transport coefficients' (χ HP and V HP ) on the harmonic number is deduced in the presence of transport hysteresis. When the contribution of the hysteresis is substantial, equation (5) shows that the dependence of amplitude on radius is common to higher harmonics (as in, e.g. figure 3(a)). If one fits the response equation (5) to the 'wave form', equation (6b), one obtains that the fitted wavenumbers (k r and k i ) do not depend on the harmonic number. Thus, the relation equation (8) leads to the properties (i) χ HP and V HP increase in proportion to the harmonic number, and (ii) the sign of fitted V HP depends on whether the temperature perturbation propagates outward or inward. The fitting of observations in figure 3 to the effective transport coefficients was performed in [9]. It is reproduced in figure 4, where the vertical shade-band indicates the frequency at which the ratio between the frequency and change rate 1/τ E(scaling) is three. The parameter 1/τ E(scaling) is chosen as a characteristic rate of change that is induced by the diffusive process. In the data analysis [9], τ E(scaling) is evaluated for each experimental condition by use of the formula for ITER89-P [13] or ISS95 [14]. If the oscillation frequency of the temperature perturbation is much higher than this criterion, the contribution of the homogeneous solution of equation (4)-which decays exponentially in radius-becomes small, and does not significantly influence the fitted values of wavenumber. In such circumstances, the fitted 'effective transport coefficients' show the properties that (i) χ HP and V HP are proportional to the harmonic number, and (ii) the sign of fitted V HP depends on the propagation direction. The results in figure 4 show these dependencies. It is highly plausible that the transport hysteresis exists, to a greater or lesser extent, in wider experimental conditions. If the diffusion model is correct, the transport coefficient must be recovered as a constant value for various harmonic numbers.
In the case of NBI heating, if the source is modulated similarly to equation (1), δT 2m−1 is predicted to be proportional to Q(2m − 1) −3 in the large-m limit. Hysteresis has also been reported for the case of NBI heating in the literature [4]. Thus, evidence of the transport hysteresis will be found in much wider circumstances. (If one could measure the response of electron temperature with very high time resolution, one might be able to study the influence of the supply of relatively cold electrons by the NBI. When the high-energy neutral particle beam is ionized in the plasma, cold electrons (the temperature of which is usually lower than the electron temperature of target plasma) are generated together with hot ions. The impact of this electron cooling is smaller by the factor O(T e /E b ) compared to the heating by NBI, where E b is the kinetic energy of beam ions. The impact of electron cooling remains a small correction to the main heating. Based on the preceding analysis, one gains an estimate that this impact of cooling may induce the additional component in higher harmonics that has a dependence on m as (2m − 1) −2 , but is smaller by the factor O(T e /E b ) compared to the dominant part. Thus, the prediction of m-dependence of (2m We here make a short note on the influence of the outward motion of plasmas (which is induced by the change of Shafranov shift as a result of the increment/reduction of the plasma pressure gradient) on the observation of transport hysteresis. In principle, the change of Shafranov shift can introduce an additional change of temperature at points fixed for measurement. (In the [15], the low-field side and highfield side measurements were performed, and an average was taken.) However, this process affects the conclusion of this article very little. As is explained quantitatively in [6], the jump of the heat flux at on-off of heating was found to occur in a very short time interval, i.e. within 1 ms, which is limited by the time resolution of the evaluation. This is much smaller than the characteristic time scale (∼ 50 ms) over which gradient (and Shafranov shift) changes. During this short time, the change of plasma pressure remains very small, so that the substantial change of Shafranov shift does not occur. The jump of the heat flux (that in the hysteresis of gradient-flux relation) is influenced very little by the periodic change of Shafranov shift. The LFS case (LHD and DIII-D) and HFS case (TJ-II and KSTAR) are analyzed here, and support a common view.

Possible hysteresis in other transport channels
Hysteresis with rapid timescale exists in core plasma transport in tokamaks and helical systems, in the channels of energy, electron and impurity densities, and plausibly in momentum transport. The mechanisms that induce nonlocal-in-space coupling work here, as has been assessed in [8], but do not suffice. The very fast time scale in change of flux indicates that the heating can directly enhance the turbulence and transport, i.e. 'the heating heats turbulence'. Rapid response of particle flux at the onset of heating [16], sudden jumps in the gradient-flux relation for carbon impurity flux [17], and dynamic hysteresis loops in velocities [18] have been observed. Interference between the particle flux and heat flux in the MECH has been investigated [19].

Physics models for hysteresis in transport relation and hydrogen isotope effect
Theories predict that 'the heating heats turbulence'. One possibility was pointed out by considering the coupling between fluctuations and the external source S in phase space associated with the heating [20]. The velocity moments of S with the weighting factor mv 2 /2 is the heating density. In turbulent plasmas, the source of heating, S, has its fluctuation part: S[ f ] = S[ f 0 ] + (dS/df )δf , where δf is the perturbation of distribution function. Thus, in the kinetic equation that describes the evolution of the fluctuation δf , a new source term (dS/df )δf appears, in addition to the conventional terms, which represent destabilizing and nonlinear-damping effects. Considering the competition between this new effect and the nonlinear-damping rate by ambient turbulence, the fluctuation intensity I is given as where I 0 is the intensity without heating effect. The heating effect is represented by the parameter Γ h = (dP/dp)/γ damp , where γ damp = χ N k 2 ⊥ is the damping rate without heating effect, k is the relevant wavenumber, P is the heating power, p is the plasma pressure, and χ N is the diffusion coefficient by ambient microscopic turbulence [20]. (The physics basis for the parameter Γ h is added here. In the new term (dS/df )δf of the kinetic equation, the time rate (dS/df ) appears, which, with an appropriate weighting function, yields the destabilizing rate by this new mechanism. For a fixed heating scheme, the coefficient (dS/df ) is proportional to (dP/dp) as is demonstrated in [20]. As a result of this, the correction of turbulence intensity by this new effect is evaluated by use of the ratio between (dP/dp) and γ damp .) The difference I − I 0 at the onset of heating may correspond to the jump in hysteresis. When P is stronger, Γ h is larger and I − I 0 becomes larger. This explains the observation in [6].
The 'isotope effect' on the direct heating effect [21] is also discussed. The local turbulent diffusivity χ N is assumed to be an increasing function of mass number A of hydrogen isotope. Then Γ h is a decreasing function of A, if other parameters are common. For instance, Γ h ∼ A −0.5 , if χ N ∼ χ gB ∼ A 0.5 and k does not depend upon ρ i . Equation (10) shows that the height of the jump in transport hysteresis is smaller for larger A. The change of the temperature gradient, against the change of power per particle, ∆q r /n, is calculated in figure 5, where χ is the local thermal diffusivity. Here the parameter α represents the relative contribution of the jump in the hysteresis, and is interrelated with Γ h as α ∼ Γ h /(1 − Γ h ). Comparing the cases of heavy (2) and light (1) hydrogen isotopes, one has where the local gyro-Bohm relation is assumed for χ N . The hysteresis can offset the isotope-dependence of local diffusion. If the condition α 2 < (A 2 /A 1 ) 0.5 α 1 + 1 − (A 2 /A 1 ) 0.5 is satisfied, the heavier hydrogen isotope can show reduced transport (better confinement) in comparison with the light hydrogen isotope.

Impacts on control system in fusion device
The hysteresis in the core transport has a strong impact on the control system against transient change of plasmas. For instance, the short interval of heating has strong influence on the edge plasma conditions. A short pulse-like change of heating in the core is smoothed, if the propagation of variation X(x, t) is a diffusive process (where the diffusivity is denoted by χ): At the distance of the order of plasma radius a, the half width of a peak in time is of the order of a 2 /χ. By this broadening, the peak height of the impact becomes very small when the pulse arrives at the edge. In contrast, when hysteresis exists, as in figures 1 and 2, in the gradient-flux relation, the sudden change propagates rapidly, without being smoothed by diffusive processes. The transport hysteresis can cause a large change in heat flux at the edge, so as to cause rapid H-L back transition and sudden heat load at the divertor. Observations in [22] indicate that the edge plasma parameter starts to change very shortly when the core heating is turned off. The transition from the H-mode with small ELMs to the ELMfree H-mode happens (almost immediately) when the core heating power is turned off. Considering the fact that the abrupt change of transport occurs when the core heating is turned on or off in W7-AS plasmas [22], it is highly plausible that the abrupt change of heat flux propagates rapidly (as was observed on DIII-D and LHD), and that rapid transmission of the change of heat flux to the edge occurs. The mechanisms of such rapid propagation have been discussed extensively, as is assessed in [7], including the influence of multiple-scale nonlinear interactions [23]. Combining studies of nonlocal dynamics and of transport hysteresis, will allow essential progress to be made. This problem of rapid transition of change might have an important influence on the dynamic properties of burning plasmas in experimental reactors. The possibility that transport hysteresis causes giant ELMs, burn through of detached plasmas, etc will be a serious issue of the edge and scrapeoff layer. The influence becomes more serious as the plasma becomes larger. The dynamics of plasma response entails a much faster time scale than diffusive change; the control system to protect the system from sudden plasma loss must be designed accordingly.

Summary
In short: (1) hysteresis in the gradient-flux relation can exist unambiguously in core plasma transport, for channels of energy, electron and impurity densities, and others; (2) theoretical modeling is in progress, and (3) the knowledge and understanding of transport hysteresis is critical for control systems of fusion devices. This might also provide a clue to understanding the hydrogen-isotope effect in plasma confinement. Figure 5. Illustration of the isotope effect in the presence of the hysteresis of the gradient flux relation. The cases for light isotopes (thick dashed line, and suffix 1) and heavy isotopes (solid line, and suffix 2) are shown. Change of gradient against variation of the heat flux per particle ∆q r /n, χ∆|∇T| = (1 − α)∆q r /n, is shown by the thick arrow on the horizontal axis.