It is known that rapid edge cooling of magnetically confined plasmas can trigger heat pulses that propagate rapidly inward. These can result in large excursion, either positive or negative, in the electron temperature at the core. A set of particularly detailed measurements was obtained in Large Helical Device (LHD) plasmas [S. Inagaki *et al*., Plasma Phys. Controlled Fusion **52**, 075002 (2010)], which are considered here. By applying a travelling wave transformation, we extend the model of Dendy *et al*., Plasma Phys. Controlled Fusion **55**, 115009 (2013), which successfully describes the local time-evolution of heat pulses in these plasmas, to include also spatial dependence. The new extended model comprises two coupled nonlinear first order differential equations for the (*x*, *t*) evolution of the deviation from steady state of two independent variables: the excess electron temperature gradient and the excess heat flux, both of which are measured in the LHD experiments. The mathematical structure of the model equations implies a formula for the pulse velocity, defined in terms of plasma quantities, which aligns with empirical expectations and is within a factor of two of the measured values. We thus model spatio-temporal pulse evolution, from first principles, in a way which yields as output the spatiotemporal evolution of the electron temperature, which is also measured in detail in the experiments. We compare the model results against LHD datasets using appropriate initial and boundary conditions. Sensitivity of this nonlinear model with respect to plasma parameters, initial conditions, and boundary conditions is also investigated. We conclude that this model is able to match experimental data for the spatio-temporal evolution of the temperature profiles of these pulses, and their propagation velocities, across a broad radial range from $r/a\u22430.5$ to the plasma core. The model further implies that the heat pulse may be related mathematically to soliton solutions of the Korteweg-de Vries-Burgers equation.

## I. INTRODUCTION

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,^{17} HL-2A,^{18} Alcator C-Mod,^{19} TEXTOR,^{20} TEXT,^{21} RTP,^{22} and TFTR,^{23} 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 $\delta \u2207Te$ and $\delta qe$ at a specific radius during two well diagnosed cold heat pulse experiments^{24} 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 $\xi =x+v0t$; here, *v*_{0} 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*_{0} that aligns with prior empirical expectations.

## II. MODEL DESCRIPTION

The normalised zero-dimensional model^{25} examined below is constructed in terms of the three key physical quantities that were measured^{24} so as to characterise pulse propagation in LHD. These quantities are the deviation from steady state of the electron temperature gradient $\delta \u2207Te$, the excess turbulent heat flux $\delta qe$, and the deviation $\delta Te$ of electron temperature from its steady-state value. The dimensionless counterparts of these variables are denoted by *x*_{1}, *x*_{2}, and *x*_{3}, 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 $\delta \u2207Te$ and $\delta 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*_{3} is slave to *x*_{2}. The change in the numerical value of *x*_{3} at each time step is calculated directly from Eq. (3) using the time-evolving value of *x*_{2} from the two-field system of Eqs. (1) and (2), in which *x*_{3} acts as an *x*_{2}-slaved coefficient. As noted, physically, the value of *x*_{3} 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*_{0} 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*_{3} is a dependent variable, we choose to simplify by neglecting all *y*_{3} 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*_{3} values are deduced directly from the time series of *y*_{2}, using Eq. (3) above. Let us now operate on Eq. (6) with $d/d\xi $ and multiply Eq. (7) by $\kappa T0$, then eliminate $dy2/d\xi $ 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,^{31} 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*_{1}, 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*_{1}.^{31} 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,^{33} which encompasses Eq. (10). An important avenue to explore is a potential link between the present approach and the theory of turbulence spreading,^{33} where the propagation velocity of a turbulence front scales as the square root of the product of a linear growth rate with heat diffusivity.

## III. COMPARISON OF MODEL RESULTS WITH LHD EXPERIMENTAL DATA

We recall that *y*_{1}, *y*_{2}, and *y*_{3} are the dimensionless counterparts of $\delta \u2207Te,\delta qe$, and $\delta 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 $\kappa T,\kappa Q$ and their derivatives, together with the numerical values of $\gamma L1,\gamma L2,\eta ,\tau c$, and *χ*_{0}, carry over identically from Ref. 25, to which we refer for the experimental motivation for these values; see Table I for detailed information.

Parameter . | . | . | . | . | . | . | . | . |
---|---|---|---|---|---|---|---|---|

Case . | $\kappa T0$ . | $\u2202\kappa T/\u2202x1=\u2202\kappa T/\u2202x3$ . | $\kappa Q0$ . | $\u2202\kappa Q/\u2202x1=\u2202\kappa Q/\u2202x3$ . | $\gamma L1=\gamma L2$ . | χ_{0}
. | $\eta =2\chi 0/3$ . | $\eta /\tau c\chi 0$ . |

T rise (R) _{e} | 15 | 1.5 | 225 | 22.5 | 35 | 3.2 | 2.13 | 10.5 |

T drop (D) _{e} | 20 | 2.0 | 400 | 40.0 | 35 | 3.2 | 2.13 | 10.5 |

Parameter . | . | . | . | . | . | . | . | . |
---|---|---|---|---|---|---|---|---|

Case . | $\kappa T0$ . | $\u2202\kappa T/\u2202x1=\u2202\kappa T/\u2202x3$ . | $\kappa Q0$ . | $\u2202\kappa Q/\u2202x1=\u2202\kappa Q/\u2202x3$ . | $\gamma L1=\gamma L2$ . | χ_{0}
. | $\eta =2\chi 0/3$ . | $\eta /\tau c\chi 0$ . |

T rise (R) _{e} | 15 | 1.5 | 225 | 22.5 | 35 | 3.2 | 2.13 | 10.5 |

T drop (D) _{e} | 20 | 2.0 | 400 | 40.0 | 35 | 3.2 | 2.13 | 10.5 |

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\u2032(\xi =0)$. Second, there is the value of the pseudo-velocity *v*_{0}. 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*_{1}, *y*_{2}, and *y*_{3}) 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(\xi =0)=y1(\xi =10)=0$; $y2(\xi =0)=\u22121.5,y2(\xi =10)=0$, and $y3(\xi =0)=0.01$ for the core electron temperature Rise (R) case in LHD plasma 49708. We take $y1(\xi =0)=y1(\xi =10)=0$; $y2(\xi =0)=1,y2(\xi =10)=0$, and $y3(\xi =0)=\u22120.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\u2032(\xi =0)=\u22121.496$ and $y2\u2032(\xi =0)=3.552$ for R case; $y1\u2032(\xi =0)=0.668$ and $y2\u2032(\xi =0)=\u22121.156$ for D case. Table II lists all boundary conditions.

Boundary condition . | . | . | . | . | . | . |
---|---|---|---|---|---|---|

Case . | $y1(\xi =0)=y1(\xi =10)$ . | $y2(\xi =0)$ . | $y2(\xi =10)$ . | $y3(\xi =0)$ . | $y1\u2032(\xi =0)$ . | $y2\u2032(\xi =0)$ . |

T rise (R) _{e} | 0 | −1.5 | 0 | 0.01 | −1.496 | 3.552 |

T drop (D) _{e} | 0 | 1.0 | 0 | −0.01 | 0.668 | −1.156 |

Boundary condition . | . | . | . | . | . | . |
---|---|---|---|---|---|---|

Case . | $y1(\xi =0)=y1(\xi =10)$ . | $y2(\xi =0)$ . | $y2(\xi =10)$ . | $y3(\xi =0)$ . | $y1\u2032(\xi =0)$ . | $y2\u2032(\xi =0)$ . |

T rise (R) _{e} | 0 | −1.5 | 0 | 0.01 | −1.496 | 3.552 |

T drop (D) _{e} | 0 | 1.0 | 0 | −0.01 | 0.668 | −1.156 |

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 $\rho =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 $\rho =0.546$ outward to $\rho =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 $\rho =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*_{0} 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^{−1} and (53.50 ± 20.97) ms^{−1}, 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\u223c(\kappa TD\kappa QD/\kappa TR\kappa QR)1/2=1.54\u22432$, 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*

_{0}in Eq. (4).

Mathematically, we may treat *v*_{0} 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 $(\xi =5;\xi =15;\xi =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.

## IV. CONCLUSIONS

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^{24} 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 $\u2243$ 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*_{0}, 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\u22430.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.

## ACKNOWLEDGMENTS

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.