A quantitative model for heat pulse propagation across large helical device plasmas

The transport of energy across magnetically confined fusion plasmas, and the storage of energy within them, reflects a wide range of turbulent and nonlinear phenomenology; see, for example, Refs. 1–32. There is extensive experimental evidence for transport phenomena that are non-diffusive and may be non-local. Examples have been found in many tokamak plasmas, for example, JET,^12,13 DIII-D,^14–16 JT-60U,¹⁷ HL-2A,¹⁸ Alcator C-Mod,¹⁹ TEXTOR,²⁰ TEXT,²¹ RTP,²² and TFTR,²³ as well as in the Large Helical Device (LHD) stellarator-torsatron.^17,24–27 A broad range of techniques for data analysis, see, for example, Refs. 12, 14, and 28, have been used to identify various forms of perturbation of heat and particle fluxes from their steady states, see the reviews.^29,30 Measurements of the spatio-temporal propagation of strongly nonlinear localised heat pulses provide a particularly interesting, and potentially fruitful, challenge to theoretical understanding and models. Here, we focus on cold pulse experiments, see, for example, Refs. 18, 19, and 24. These cold pulses are typically initiated by injection into the edge plasma of ice pellets or supersonic molecule beams, by gas puffing, and by laser ablation. The zero-dimensional model of Ref. 25, which incorporates only time dependence, is successful in quantitatively capturing the local time evolution of $δ∇Te$ and $δqe$ at a specific radius during two well diagnosed cold heat pulse experiments²⁴ using pellet injection in LHD; we refer, in particular, to Figs. 3 and 5 of Ref. 25. This motivates the following physical conjecture, which we test in the present paper. The zero-dimensional model is known to work at the best diagnosed spatial location, capturing the time evolution of the pulse there, where t = 0 is defined to be the local arrival time of the initial impulsive perturbation. Therefore, we conjecture that the zero-dimensional model ought to work at each location across the radial domain of the plasma within which the same physical processes determine the behaviour of the pulse. From this, we infer that the model ought to apply in a frame that is co-moving radially with the heat pulse across this region. This final step in the conjecture provides a simple path to construct a spatio-temporally dependent model from the tested time-dependent-only model of Ref. 25. We apply a travelling wave transformation by replacing t in Eqs. (1–3) of Ref. 25 by $ξ=x+v0t$⁠; here, v₀ is to be considered as a proxy of the pulse propagation velocity in the radial direction x, and the sign convention adopted for t assists consideration of inward propagation. Importantly, as we shall show, the mathematical structure of the resultant model equations, combined with the choice of physical model parameters that carries over from Ref. 25, yields a formula for v₀ that aligns with prior empirical expectations.

The normalised zero-dimensional model²⁵ examined below is constructed in terms of the three key physical quantities that were measured²⁴ so as to characterise pulse propagation in LHD. These quantities are the deviation from steady state of the electron temperature gradient $δ∇Te$⁠, the excess turbulent heat flux $δqe$⁠, and the deviation $δTe$ of electron temperature from its steady-state value. The dimensionless counterparts of these variables are denoted by x₁, x₂, and x₃, respectively, as defined in Ref. 25. The model of Ref. 25 shows quantitative agreement between its outputs and experimental measurements of the time evolution of these variables at fixed locations in the LHD plasma, after rapid cooling at the edge. It is successful both when core electron temperature rises on arrival of the pulse and when it drops. However, for given parameters and initial conditions, the model only simulates the time evolution of a passing heat pulse at a specific radius, for example, $r/a=0.19$ in Ref. 25. Spatial dependence is eliminated from the physical picture underlying the model in Ref. 25 by using the parameter $1/Lc$ as a proxy for divergence in the heat flux energy conservation equation; here, L_c is the characteristic scale-length of steady-state turbulent transport. The model equations (5) to (7) of Ref. 25 are as follows:

Physically, Eqs. (1) to (3) embody on assumption of dominant strong coupling between $δ∇Te$ and $δqe$⁠, while the background transport acts to dissipate the pulse. The physical significance of the various coefficients is described in the discussion of equations (1)–(4) of Ref. 25. We note that the mathematical structure of Eq. (3) above implies that x₃ is slave to x₂. The change in the numerical value of x₃ at each time step is calculated directly from Eq. (3) using the time-evolving value of x₂ from the two-field system of Eqs. (1) and (2), in which x₃ acts as an x₂-slaved coefficient. As noted, physically, the value of x₃ maps to the deviation in electron temperature over time, which is a primary experimentally measured time-dependent variable for the model comparisons reported below.

Central to the present paper is the adoption of travelling wave transformations as a method for generating (x, t)-dependence from the t-dependent model of Eqs. (1) to (3). We assume that

Here, x and t are considered to be independent variables, and we refer to v₀ as the pseudo-velocity of the pulse. This pseudo-velocity is expected to be similar, but not identical, to the real measured velocity of the pulse. It then follows that

Due to the fact that y₃ is a dependent variable, we choose to simplify by neglecting all y₃ related terms in Eqs. (1) to (3), and substitute Eq. (5) into Eqs. (1) and (2), yielding

Our new model [Eqs. (6) and (7)] comprises only two independent variables: y₃ values are deduced directly from the time series of y₂, using Eq. (3) above. Let us now operate on Eq. (6) with $d/dξ$ and multiply Eq. (7) by $κT0$⁠, then eliminate $dy2/dξ$ by substitutions. This yields the following equation after leading order approximation and transpositions:

The second coupled nonlinear ordinary differential equation is derived by applying similar procedures

Equations (8) and (9) comprise our mathematical model for the propagating heat pulse. Its physical motivation is that of Ref. 25, combined with the co-moving conjecture outlined above and expressed in Eq. (4). The mathematical structure of Eq. (8) can be linked, in certain approximations, to the Korteweg-de Vries-Burgers (KdV-Burgers) equation,³¹ and Eq. (9) can be linked to a damped wave equation in the co-moving frame. In particular, the left hand side of Eq. (8) corresponds to that in the formulation of the KdV-Burgers equation in Eq. (5.3) of Ref. 31. This nonlinear differential operator, acting on y₁, is known to support soliton pulses, and furthermore these pulses move at a speed determined by the coefficient of the term which multiplies the linear term in y₁.³¹ The expression for this coefficient in Eq. (8) motivates the following conjecture regarding heat pulse velocity in our model:

This scaling has previously been noted empirically from heat pulse experiments.^29,32 It is possible that there will in future arise a convergence between the mathematical structure of the present model, empirical scalings,^29,32 and other theoretical work,³³ which encompasses Eq. (10). An important avenue to explore is a potential link between the present approach and the theory of turbulence spreading,³³ where the propagation velocity of a turbulence front scales as the square root of the product of a linear growth rate with heat diffusivity.

We recall that y₁, y₂, and y₃ are the dimensionless counterparts of $δ∇Te,δqe$⁠, and $δTe$ respectively. Solution of our model, embodied in the two coupled nonlinear ordinary differential equations (8) and (9), proceeds as follows. The numerical values for $κT,κQ$ and their derivatives, together with the numerical values of $γL1,γL2,η,τc$⁠, and χ₀, carry over identically from Ref. 25, to which we refer for the experimental motivation for these values; see Table I for detailed information.

The new model embodied in Eqs. (8) and (9) is strongly nonlinear. In solving it, we have three kinds of degrees of freedom. First, there are the initial values of $yi′(ξ=0)$⁠. Second, there is the value of the pseudo-velocity v₀. Third, it is possible to apply fixed horizontal and vertical shifts in the values of the electron temperature and time, provided these shifts are applied uniformly to all outputs of a simulation, as a way of reducing systematic errors introduced by parameter choices. We optimise model outputs, for comparison with the experimental measurements, in these three ways. These model outputs (time traces of y₁, y₂, and y₃) reflect the underlying phase space structure of the solutions of the nonlinear system of equations. This structure is robust, in the sense discussed in Ref. 25, and as demonstrated for the present system of equation in Fig. 7 below, for example. We recall also that the great majority of model parameter values, see, e.g., Table I, are determined by experimental measurements.

The boundary conditions on the y_i, i = 1, 2, 3, and their derivatives, evaluated at ξ = 0 and ξ = 10 which define the solution domain, coincide with the experimentally motivated values in Ref. 25 where these carry over. Specifically, $y1(ξ=0)=y1(ξ=10)=0$⁠; $y2(ξ=0)=−1.5,y2(ξ=10)=0$⁠, and $y3(ξ=0)=0.01$ for the core electron temperature Rise (R) case in LHD plasma 49708. We take $y1(ξ=0)=y1(ξ=10)=0$⁠; $y2(ξ=0)=1,y2(ξ=10)=0$⁠, and $y3(ξ=0)=−0.01$ for the core electron temperature Drop (D) case in LHD plasma 49719. Boundary conditions unspecified in Ref. 25, which are needed in order to satisfy the boundary conditions above, are assumed to be as follows here: $y1′(ξ=0)=−1.496$ and $y2′(ξ=0)=3.552$ for R case; $y1′(ξ=0)=0.668$ and $y2′(ξ=0)=−1.156$ for D case. Table II lists all boundary conditions.

Figure 1 compares time traces of the evolving electron temperature at multiple radial locations, obtained from the model and from experimental data(# Te49708) for the R case. Several representative radii are marked by arrows on the right hand side. The model results are able to match experimental data from $ρ=0.450$ inward to the core, if we uniformly apply horizontal ( + 0.01) and vertical (+0.20) shifts, suggesting that the model applies over this broad radial range. It is also clear that electron temperature profiles from $ρ=0.546$ outward to $ρ=0.703$ are not simulated by the model. This suggests that different physics dominates heat pulse propagation in the outer region of this plasma. Figure 2 demonstrates the comparison of model results and experimental data(# Te49719) for the D case. In common with the R case shown in Fig. 1, the model results are good from $ρ=0.450$ inward to the core, with uniformly applied horizontal (+0.04) and vertical (+0.07) shifts, but not in the outer region of this plasma. Figs. 3 and 4 show zoomed-in electron temperature pulse plots at three specific radii, selected from Figs. 1 and 2, respectively.

Empirically, we approach v₀ by identifying a best fit straight line linking the positions of the peak or trough of the heat pulse temperature profile in the two cases, see Figs. 5 and 6. The effect of data noise on peak location is minimised by taking a five per cent running window. The width of the filter window determines the horizontal error bars. We assume, for simplicity and to match our travelling-wave model, that the inward radial propagation velocity of the heat pulse is invariant across the plasma volume of interest. A straight line fit to the radial location of the pulse peak versus time is thus applied. The error of this fitting may come from the radial dependence of the pulse velocity, which tends to increase in the central region. We find from Figs. 5 and 6 that $v0=15$ for the temperature rise case and $v0=30$ for the temperature drop case. These dimensionless pseudo-velocity values equate to the dimensional values (32.62 ± 9.89) ms⁻¹ and (53.50 ± 20.97) ms⁻¹, respectively, which are comparable in magnitude to measured heat pulse propagation velocities. The ratio 2 between these two velocities is broadly consistent with the value inferred by substituting experimentally measured values for transport coefficients into Eq. (10). We have $v0D/v0R∼(κTDκQD/κTRκQR)1/2=1.54≃2$⁠, referring also to Table I, where superscripts R and D denote the core T_e rise and drop cases, respectively. We note that the empirically determined velocities $v0R$ and $v0D$ of heat pulse propagation need not necessarily coincide with the optimal value of the mathematical transformation parameter v₀ in Eq. (4).

Mathematically, we may treat v₀ as a free parameter in Eqs. (8) and (9), which we have labelled the pseudo-velocity. We solve Eqs. (8) and (9) repeatedly for different values of this pseudo-velocity, see, for example, Fig. 7 for the R case, and identify the best fit value. To test sensitivity with respect to ξ, three other cases $(ξ=5;ξ=15;ξ=20)$ have been examined. All three of these test options exhibit the same properties as the case of ξ = 10. Various combinations of initial conditions and pseudo-velocities were also tested, showing that the combination specified above provides the best fit to the experimental data. For the phase plots, see, for example, the right-hand panel of Fig. 7, circulation directions are identical with those of the experimental data for both the R and D cases.

The fitting operation can be assisted by shifting the origin of co-ordinates. The robustness, with respect to variation of initial conditions, of this mathematical approach to identifying pseudo-velocity has been tested. This robustness reflects the structure of the underlying attractor in phase space, projected in $(y1,y2)$ coordinates in these two figures; compare also Figs. 6 and 7 of Ref. 25.

We have derived from first principles a time-dependent model in one spatial dimension, which is able to describe quantitatively the radial inward propagation of heat pulses in the core of two plasmas²⁴ in the LHD. In one plasma the central electron temperature rises, in the other it falls. This new model is derived from a travelling wave transformation of the zero-dimensional model of Ref. 25, which is known to capture the time-evolution of the heat pulse as it passes through a fixed radial location in these two plasmas. From the experimental data, we infer that the velocity of the propagating pulse is constant in both the electron temperature rise and drop cases. The pulse velocity in the electron temperature rise case is smaller than in the drop case by a factor $≃$ 2. This aligns with Eq.(10), which reaches back into the mathematical structure of our model, and also coincides with empirical expectations given the values of the heat conduction parameters for these two plasmas. A pseudo-velocity parameter is introduced in the travelling wave transformations, in order to model heat pulse propagation across spatial location as well as in time. From numerical tests, we discover that real pulse velocity is about two times the best estimate of the travelling-wave transformation parameter v₀, referred to as the pseudo-velocity. Comparison between model outputs and raw experimental data suggests that our model is able to describe heat pulse propagation well, within a broad radial range of the LHD core plasma from $r/a≃0.5$ to the centre.

The results of the present paper provide additional support to the physical proposals, described in Sec. 2 of Ref. 25, which motivate the simple model equations reproduced as Eqs. (1) to (3) above. Central to these proposals is the conjecture that heat pulses are structures which involve plasma physics processes which are so strongly nonlinear that heat pulse evolution is primarily determined by the reactions of the perturbed heat flux and the perturbed temperature gradient on each other. In the present model, this is the dominant element of nonlinear physics, whereas conventional turbulent transport plays a relatively minor dissipative role. It is this mutually coupled interaction between perturbed heat flux and temperature gradient that governs the local plasma dynamics of the heat pulse in space, equivalent to the local up-and-down dynamics of a water wave under gravity. We have shown in this paper that this coupling model lends itself readily to a travelling wave transformation, yielding spatio-temporal pulse propagation, and that the pulse velocity that emerges mathematically provides an adequate match to empirical results and expectations. This aspect of the analysis also provides guidance on a previously unanswered question, namely, the generic character of the heat pulse: we have shown that it may be closely related to a Korteweg-de Vries-Burgers soliton. Two avenues of investigation would repay immediate attention. First, this model has been tested on, and motivated by, measurements from only two plasmas in LHD—albeit plasmas with exceptionally high quality measurements. The simplicity of the model encourages one to hope that it may be more widely applicable and clearly it should now be tested on a broader range of heat pulse experimental datasets, provided that they possess the required spatial and temporal resolution and that the other relevant plasma parameters are well diagnosed, as in LHD. Second, while gyrokinetic or other computationally intensive transport simulations have not yet (to the authors' knowledge) been applied to heat pulse experiments, the outputs of any future simulations of heat pulses could be tested directly against the analytical model presented here, which is constructed in terms of variables that are directly measurable, as distinct from the first-principles particle and field distributions of gyrokinetic theory.

This work was part-funded by the RCUK Energy Programme under Grant No. EP/I501045. This project has received funding from the European Union's Horizon 2020 research and innovation programme under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This work was also supported in part by the grant-in-aid for scientific research of JSPF, Japan (23244113) and KAKENHI (21224014) from JSPS, Japan. This work was also supported in part by the Grant-in-Aid for Scientific Research of JSPF, Japan (23360414) and by the collaboration programs of the National Institute for Fusion Science (NIFS) (NIFS13KOCT001) and the Research Institute for Applied Mechanics of Kyushu University.