Modelling the measured local time evolution of strongly nonlinear heat pulses in the Large Helical Device

In some magnetically confined plasmas, an applied pulse of rapid edge cooling can trigger either a positive or negative excursion in the core electron temperature from its steady state value. We present a new model which captures the time evolution of the transient, non-diffusive local dynamics in the core plasma. We show quantitative agreement between this model and recent spatially localized measurements (Inagaki et al 2010 Plasma Phys. Control. Fusion 52 075002) of the local time-evolving temperature pulse in cold pulse propagation experiments in the Large Helical Device.


Introduction
Understanding, prediction and control of energy transport is central to achieving nuclear fusion in magnetically confined plasmas, which are large scale physical systems, nonlinearly coupled across a broad range of spatio-temporal scales and far from equilibrium. They are typically turbulent and can exhibit energy transport phenomenology which is nonlinear and non-diffusive [1][2][3][4][5][6][7][8][9]. At-a-point measurements of plasma parameters during cold heat pulse experiments in tokamak and stellarator plasmas [10][11][12][13][14][15][16][17][18][19] probe the underlying transport processes in a strongly perturbed regime, and are an unresolved challenge to interpretation. We note that the nature and implications of a local measurement of transport may differ significantly between tokamak and stellarator plasmas. For example, in a large aspect ratio tokamak, local magnetic shear is a good proxy for flux-surface-averaged magnetic shear, but this is not the case in a stellarator such as the Large Helical Device (LHD). Insofar as magnetic shear affects transport, this may be important. In these experiments [10][11][12][13][14][15][16][17][18][19], the 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. electron temperature T e at the plasma edge is rapidly reduced by a transient local increase in radiation, typically induced by injection of a pellet which produces radiating impurity ions and cold electrons. The magnitude of the local gradient of electron temperature ∇T e at the edge is correspondingly increased. The subsequent measured local behaviour of the plasma, as the resultant cold pulse propagates rapidly inwards, cannot be understood in terms of either diffusive or convective transport. Paradoxically, an applied pulse of rapid edge cooling can trigger either a positive or negative excursion in the core electron temperature from its steady state value, depending on the confinement properties of the plasma [11][12][13]19]; see, for example, figure 1. Spatially non-local transport properties have been inferred from cold pulse experiments, and critical temperature gradient scale length models and empirical nonlocal models have been tested [20][21][22]. These models, which often incorporate large scale transport codes, contribute to explaining various transient transport phenomena. However successful analytical models and numerical simulations based on the fundamental equations of plasma physics remain elusive. These issues remain are topical, and their potential significance is growing. For example, Rice et al [23] have recently indicated that cold heat pulse effects observed in Alcator C-Mod may be closely connected to other seemingly unrelated transport effects, including core toroidal rotation reversals, energy confinement saturation and up/down impurity density asymmetries. Mantica et al [24] used cold pulses in JET to probe internal transport barriers. In recent reviews, cold pulse experiments are placed in the broader context of perturbative transport phenomena in fusion plasmas in [25], and theoretical approaches to nonlocal transport in tokamaks and stellarators are described and extended in [26].
Here we focus on understanding how the local core plasma responds, over time, to an applied pulse of rapid edge cooling, by exploiting newly available at-a-point measurements in the LHD [19]. In section 2, we propose a model comprising a set of three coupled nonlinear equations that depend on time, with spatial dependence entering implicitly through the local values of transport coefficients, which we infer from experimental measurements. Using experimental parameters as inputs, we show in section 3 that the system dynamics output from the model corresponds quantitatively to the observations. This suggests that we have captured the essential local time-evolving physics of the strongly nonlinear plasma response. The new physically motivated model is the first simple analytical model for local thermal evolution which captures both types of observed dynamics (central temperature rise or fall), which a successful model must predict. This approach may find wider application in understanding the pulsed nonlinear transport dynamics of other non-plasma macroscopic systems that are far from equilibrium.

Model equations
We propose here a model which we construct in terms of the three key observed thermal properties of the core plasma: the deviation from steady state of the heat flux, of the electron temperature gradient, and of the electron temperature. These three variables are nonlinearly coupled to each other. In our model, the steady state turbulent plasma transport processes contribute only damping of the strongly nonlinear evolution of the system. This starting point is different from that of transport models constructed in terms of fundamental quantities of plasma physics, such as perturbed probability density functions, and from transport models constructed in terms of corrections to linear transport coefficients. Our model is motivated by recent experiments [17][18][19] of cold pulse propagation from edge (ρ = r/a = 1) to deep core (ρ = 0.11) of LHD, where a denotes the local minor radius (on average a = 0.6 m) and r denotes distance from the axis of symmetry. Figure 1 encapsulates the results of these experiments, summarising the electron cyclotron emission (ECE) measurements of [19]. These track the simultaneous evolution in time, in response to the initial edge cooling, of two variables observed deep in the core at plasma radius ρ = 0.19. First, there is the measured excess of the local electron heat flux relative to its steady state value, normalized to local electron number density, represented by δq e (ρ, t)/n e defined by equation (1) of [19]. Second, there is the measured excess of the local electron temperature gradient relative to its steady state value, denoted by δ∇T e (ρ, t). The LHD experiments show three distinct types of dynamics in these variables. Below a threshold in ∇T e there is diffusive transport, in which δq e /n e and δ∇T e rise and fall together. Above this threshold, there is a bifurcation in dynamics; for similar ∇T e the core plasma executes a Lissajous figure in δq e /n e and δ∇T e which encompasses either a local fall, or rise, in T e . This bifurcation is conditioned by the magnitude of δq e /n e , and the phenomenology can also depend on target plasma conditions and on the magnitude of the edge cold pulse [19]. The nondiffusive nature of the energy transport is also indicated by measurements of the time evolution of T e at multiple ECE channels across the minor radius [19]. The edge cold pulse is observed to have no apparent effect on the value of T e or δ∇T e midway across the minor radius, while the edge and core responses are separated by a time interval of a few milliseconds. This dynamical richness suggests non-trivial underlying equations.
Here we will show that the key quantitative aspects of this observed local time-evolving phenomenology in the LHD core plasma can be captured by the following physically motivated mathematical model, embodied in a system of three coupled nonlinear differential equations for the experimentally measured variables δq e (ρ, t)/n e , δ∇T e , and T e − T e0 . We assume that the excess turbulent heat flux δq e (ρ, t)/n e acts to reduce the magnitude of the local electron temperature gradient by carrying away thermal energy. Since the sign of the spatial gradient of temperature is negative, this corresponds to reducing the magnitude of the negative excess δ∇T e , that is, driving it in the positive direction. In parallel, we assume that the steady state turbulent transport acts to damp any non-zero δ∇T e , at a rate γ L1 . Hence our first model equation is The vector heat flux has only one component, and similarly for the gradient operator. To normalize the dimensionally different variables δq e /n e and δ∇T e , we use the measured steady state turbulent heat diffusivity χ 0 , which is estimated by power balance analysis and is basically the same as the slope of figure 1(d) (the metric coefficient of magnetic flux surface is needed). The value of χ 0 also corresponds approximately to χ 0 = L 2 c /τ c , where L c is the characteristic scale-length of steady state turbulent transport (∼a) and τ c is the associated global energy confinement time. The coefficient κ T , which has dimension inverse time, depends on both the local electron temperature T e and on δ∇T e , which we treat as independent variables. The excess heat flux δq e is a proxy for excess turbulence, and this is driven by excess (negative) temperature gradient and damped by the steady state turbulent transport, so we assume We adopt a similar philosophy to [29] in our third model equation, for the local time derivative of the electron temperature difference T e − T e0 , in assuming that it approximately matches the divergence of the local excess turbulent heat flux, on energy conservation grounds. We again assume that steady state turbulent transport acts to damp the deviation-here the temperature excess: where η is the transient heat diffusivity whose value is taken to be identical to 2χ 0 /3, thus nonlinearity in heat diffusivity is not introduced in our model. In the divergence, we will replace ∇ by 1/L c . The coupling factor κ T (T e , ∇T e ) can be Taylor expanded (we will see that this approximation is supported by the experimental observations) about its steady state turbulence value κ T0 = κ T (T e = T e0 , δ∇T e = 0), giving to first order In the same way, κ Q (T e , ∇T e ) can be expanded about κ Q0 = κ Q (T e = T e0 , δ∇T e = 0). Using these expressions in equations (1) and (2), and defining dimensionless variables x 1 = L c δ∇T e /T 0 , x 2 = L c δq e /χ 0 n e T 0 and x 3 = (T e − T e0 )/T 0 , our system of equations becomes We note the strongly nonlinear character of this coupled system of equations, for which we know of no analytical solution. For these two reasons, the time evolution and final state of the system cannot be inferred from the initial conditions, other than by numerically solving the model as in the next section. In this respect our model differs fundamentally from approaches based on, for example, linear diffusion. To summarize briefly, the LHD experiments outlined in Section 1 find that if the initial heat flux jump on arrival at ρ = 0.19 is negative, the electron temperature rises, and vice versa. A challenge to theory is to find a model which reproduces this together with the time evolution of the three key experimental observables which are our model parameters, given the local initial conditions and the measured local plasma parameter values. In the next section, we show that our model in equations (5) to (7) achieves this.

Comparison of model outputs with LHD results
In our approach to modelling local plasma evolution in time but not space, the arrival of the pulse at a particular radial location is treated as an impulsive change in the local heat flux, local electron temperature gradient, and local electron temperature, all defined relative to their steady state values. This determines the initial conditions for our implementation of equations (5) to (7). We now solve numerically the system equations (5) to (7), using values for coefficients that are inferred from the experimental measurements [19] on LHD plasmas. Figure 2 compares our model results to the corresponding experimental results in figure 1(b) of [19], where the temperature rises and then declines to the steady state value. The correspondence between the model output and experimental time series is good given that these are not fitted curves: both the model parameters, and the initial conditions, are determined independently from the plasma conditions in LHD. Measured values of χ 0 , L c and T e0 are employed for normalization of the time traces. The interaction between core heat flux and edge temperature gradient discussed in [19], which results in a local heat flux jump without change of local temperature gradient, provides the initial values of (x 1 , x 2 , x 3 ) = (0, −1.5, 0). Coefficients in the model are estimated from direct experimental measurements, together with other physical considerations, as follows. Coefficients κ T0 and κ Q0 are the primary determinants of the characteristic time scale of variation of the measured variables, for which the observed value is ∼60 ms, and √ κ T0 κ Q0 ∼ 100 s −1 .
Gyrokinetic-like dependence of thermal diffusivity on T e and ∇T e is indicated by transient transport experiments in LHD, and κ T and κ Q are assumed to depend similarly on T e and ∇T e [27,28]. Hence we take ∂κ T /∂x 1 = ∂κ T /∂x 3 = 0.1κ T0 , ∂κ Q /∂x 1 = ∂κ Q /∂x 3 = 0.1κ Q0 . There is thus empirical support for the truncation, to first order, of the Taylor expansion in equation (4)   parametric coordinate, as in [19]. Figure 3(a) shows the model evolution following the initial local negative perturbation in heat flux which is anticlockwise. This is overlaid with the experimental data in figure 3(b), which shows that both the direction of rotation and the approximate timing are consistent between model and experiment. Figure 1 further shows that in the non-diffusive transport regime associated with larger heat fluxes, there is a second type of core temperature response with respect to the edge perturbation of LHD, in which the core temperature initially falls rather than rises. An important test of our model is that, with its coefficient values changed only in the way determined by the changed experimental conditions, it should reconstruct the time evolution of heat flux and T e gradient in this case also. Encouragingly, the correspondence between model and observations is equally close in both cases. The results are displayed in figures 4 and 5, which show the overlaid model and experiment timeseries and phase portrait in the same format as figures 2 and 3. Again, these are not fitted curves. In figures 4 and 5 the initial conditions differ substantially from figures 2 and 3 in that the initial heat flux jump is positive (x 2 = 0.5). Slightly different coefficients κ T0 = 20 and κ Q0 = 400, which reflect the differences in T e0 and ∇T e0 , are used, but the same T e and ∇T e dependences of κ T (κ Q ) as above are assumed: ∂κ T /∂x 1 = ∂κ T /∂x 3 = 0.1κ T0 , ∂κ Q /∂x 1 = ∂κ Q /∂x 3 = 0.1κ Q0 . Damping rates γ L1 , γ L2 and η/τ c χ 0 identical to those in figure 2 are employed. Both forms of strongly nonlinear response of the experimental system are thus captured by our model, by using as initial condition the locally observed impulsive change in heat flux, together with values for model coefficients that are inferred from the experimental measurements [19] on LHD.

Discussion
In our model, non-diffusive nonlinear coupling between heat flux and temperature gradient is essential. This coupling is embodied in the coupling parameters, κ T and κ Q . Whereas values for the model parameters associated with quasiequilibrium transport properties, such as χ 0 and γ L1 , can be evaluated directly from experimental measurements, this is more difficult for the newly introduced coupling parameters. We have calculated the values of the coupling parameters as follows. First, we note that if the coupling parameters were independent of T e and ∇T e , a damped oscillatory solution of equations (1) and (2) would exist, proportional to exp(−t/τ ) sin(ωt), where ω 2 = κ T0 κ Q0 and 1/τ = γ L1 . From the experimental measurements of the oscillation period and duration of the pulses in δq e and δ∇T e , we can therefore estimate the possible range of values of κ T0 κ Q0 and γ L1 in this simplified limit. Figure 6 demonstrates that the output of our model depends strongly on the values of κ T0 and κ Q0 . We can minimize the difference between experimental and model results by choosing optimal values for κ T0 and κ Q0 from within this range. These values for κ T0 and κ Q0 are found to lie between ∼1/τ c and ∼10/τ c .
In addition, we expect the dependence of coupling parameters on local thermodynamic variables to play a role. Taking the gyro-Bohm dependence of thermal diffusivity as a conceptual starting point, we assume that the coupling parameters κ T and κ Q depend on the local temperature and its spatial gradient. An example of the relative insensitivity of our model output to the T e and ∇T e dependence of the coupling parameters is shown in figure 7. This plots the normalized fluxgradient diagram inferred from the model, using different T e and ∇T e dependence for the coupling parameters. Importantly, we find that even if there is strong T e and ∇T e dependence, the output of the model is qualitatively similar to the case where there is no dependence. We note that dependence of the coupling parameters on the remaining thermodynamic variable in the present context, namely the heat flux, would imply significantly different physics, in particular a requirement for non-local dependence, as indicated in [17].

Conclusions
The physical model represented by equations (5) to (7) is mathematically simple in its structure, and is constructed in terms of variables that are directly observed and inferred experimentally. It is motivated by the strongly nonlinear plasma response to perturbations acting on an already turbulent system seen, in particular, in LHD. There is quantitative agreement between our new analytical model and local measurements of time-evolving heat fluxes, temperature gradients, and temperature excursions. So far as we are aware, this is a novel result. Our model may have interpretive and predictive potential for fusion plasmas beyond LHD, insofar as the recent LHD results [17][18][19] represent state-ofart measurements of heat pulse phenomenology that may share features in common with nonlinear heat pulse measurements on other experiments [10][11][12][13][14][15].
Modelling the response to a nonlinear heat pulse is important for the development of fundamental understanding of the equilibrium and non-equilibrium transport processes which govern energy confinement in magnetically confined fusion plasmas.
These transport processes arise from phenomenology that is also found beyond the plasma state: multi-scale coupled physics, incorporating turbulence and coherent nonlinear structures such as zonal flows, in a spatially inhomogeneous macroscopic driven-dissipative system that has many degrees of freedom. The approach taken here to modelling strongly nonlinear perturbed transport phenomenology in LHD may therefore be of wider physical interest, in at least two respects. First, our model embodies a separation of timescales between the slow pulse evolution and the fast initial change in local temperature gradient. This separation of timescales is characteristic of anomalous or bursty transport. Second, we model the dynamics in terms of deviations from steady state, not necessarily from equilibrium. This may be helpful in modelling the dynamics of other systems that are far from equilibrium.