Skip to main content Accessibility help


  • Access
  • Open access


MathJax is a JavaScript display engine for mathematics. For more information see
      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Effect of plasma elongation on current dynamics during tokamak disruptions
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Effect of plasma elongation on current dynamics during tokamak disruptions
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Effect of plasma elongation on current dynamics during tokamak disruptions
        Available formats
Export citation


Plasma terminating disruptions in tokamaks may result in relativistic runaway electron beams with potentially serious consequences for future devices with large plasma currents. In this paper, we investigate the effect of plasma elongation on the coupled dynamics of runaway generation and resistive diffusion of the electric field. We find that elongated plasmas are less likely to produce large runaway currents, partly due to the lower induced electric fields associated with larger plasmas, and partly due to direct shaping effects, which mainly lead to a reduction in the runaway avalanche gain.

1 Introduction

Magnetic reconnection events in tokamaks often result in a sudden cooling of the plasma associated with an increase in the plasma resistivity, which in turn induces an electric field. If this electric field is larger than a certain critical electric field, electrons in the tail of the bulk Maxwellian distribution run away and acquire relativistic energies (Wilson 1925; Dreicer 1959). These runaway electrons can multiply by producing additional runaway electrons in close collisions with the thermal electrons, leading to an exponential increase of the number of runaways – an avalanche (Jayakumar, Fleischmann & Zweben 1993).

Plasma terminating disruptions in tokamaks often result in electric fields larger than the critical field. Since the avalanche production is exponentially sensitive to the initial plasma current, the problem with runaway generation is expected to be more serious in future tokamaks with larger plasma currents than in present experiments (Rosenbluth & Putvinski 1997). Uncontrolled loss of these high energy electrons may cause serious damage to plasma facing components.

Most studies of runaway current dynamics have been performed assuming circular magnetic flux surfaces, although both present and future devices often operate with elongated plasmas. In particular, future tokamaks with a large plasma current, such as ITER and SPARC (Greenwald et al. 2018), will have a significant elongation. Experimental observations indicate that runaway beams are more easily produced in limiter or low-elongation discharges than in more elongated ones (Izzo, Humphreys & Kornbluth 2012; Hollmann et al. 2013; Reux et al. 2015; Breizman et al. 2019). It is not clear if this is due to the elongation itself or to the difference of magnetic topology, i.e. limited versus divertor configuration. Magnetohydrodynamic (MHD) simulations show that differences in the MHD activity during the thermal quench phase produce better confinement of seed runaway electrons if the plasma is limited than if it is diverted (Izzo et al. 2011). Apart from these differences in MHD stability, there are also differences in the induced toroidal electric field and associated runaway current dynamics that depend on the plasma elongation directly, rather than indirectly via its effect on MHD stability.

In this paper, we focus on such effects of elongation on the coupled dynamics of runaway current generation and resistive electric field diffusion. We derive an equation governing the evolution of the toroidal electric field in a general magnetic geometry and consider the effect of elongation on the current dynamics in the case of an axisymmetric magnetic field with elliptical flux surfaces in the large aspect ratio limit. We show that the effect of elongation is to reduce both the maximum electric field, leading to a lower runaway generation rate, as well as the potential runaway avalanche multiplication. We then illustrate the effect of elongation in simulated idealized disruptions, with a pre-described exponentially decaying temperature evolution.

2 Evolution of toroidal electric field

The magnetic field in a general toroidal geometry can be written as

(2.1)$$\begin{eqnarray}\boldsymbol{B}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D712},\end{eqnarray}$$

where $\unicode[STIX]{x1D713}$ and $\unicode[STIX]{x1D712}$ are the toroidal and poloidal fluxes divided by $2\unicode[STIX]{x03C0}$, and $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D711}$ are magnetic coordinates corresponding to poloidal and toroidal angles, respectively. We are interested in the evolution of the flux-surface average of the toroidal component of the induced electric field $\boldsymbol{E}=-\unicode[STIX]{x2202}\boldsymbol{A}/\unicode[STIX]{x2202}t-\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$, where $\boldsymbol{A}$ and $\unicode[STIX]{x1D719}$ are the vector potential and the electrostatic potential, respectively. Following appendix A of Boozer (2018), we find that the electric field can be written as

(2.2)$$\begin{eqnarray}\boldsymbol{E}=\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D712}}{\unicode[STIX]{x2202}t}\right)_{\unicode[STIX]{x1D713}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}-\boldsymbol{u}\times \boldsymbol{B}-\unicode[STIX]{x1D735}(s+\unicode[STIX]{x1D719}),\end{eqnarray}$$

with $\boldsymbol{u}=(\unicode[STIX]{x2202}\boldsymbol{r}/\unicode[STIX]{x2202}t)_{\unicode[STIX]{x1D713},\unicode[STIX]{x1D703},\unicode[STIX]{x1D711}}$ the velocity of the canonical coordinates $(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})$ through space, $s=\unicode[STIX]{x1D713}\dot{\unicode[STIX]{x1D703}}-\unicode[STIX]{x1D712}\dot{\unicode[STIX]{x1D711}}$, and the dot denotes a time derivative taken at constant $\boldsymbol{r}$. The projections along $\boldsymbol{B}$ of the two last terms in (2.2) vanish upon flux-surface averaging, thus the evolution of the parallel electric field in an arbitrarily shaped tokamak or stellarator is given by

(2.3)$$\begin{eqnarray}\langle \boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B}\rangle =\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D712}}{\unicode[STIX]{x2202}t}\right)_{\unicode[STIX]{x1D713}}\langle \boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}\rangle ,\end{eqnarray}$$

where the flux-surface average is defined as

(2.4)$$\begin{eqnarray}\langle \unicode[STIX]{x1D701}\rangle (\unicode[STIX]{x1D712})=\left.\int \unicode[STIX]{x1D701}(\unicode[STIX]{x1D712},\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})\,\text{d}V^{\prime }\!\right/\!\int \,\text{d}V^{\prime },\end{eqnarray}$$

where the integrals are taken over the volume between two neighbouring flux surfaces and the volume element has been written as $\text{d}V^{\prime }=(\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703})^{-1}\,\text{d}\unicode[STIX]{x1D713}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D711}$.

In an axisymmetric magnetic field, we can write (Helander & Sigmar 2005)

(2.5)$$\begin{eqnarray}\boldsymbol{B}=F(\unicode[STIX]{x1D712})\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D712},\end{eqnarray}$$

which simplifies the toroidal electric field to

(2.6)$$\begin{eqnarray}\langle \boldsymbol{E}\boldsymbol{\cdot }\boldsymbol{B}\rangle =F\langle R^{-2}\rangle \left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D712}}{\unicode[STIX]{x2202}t}\right)_{\unicode[STIX]{x1D713}},\end{eqnarray}$$

where $R=|\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}|^{-1}$ denotes the major radius.

If we neglect the contribution of the plasma pressure in the equilibrium equation $\boldsymbol{j}\times \boldsymbol{B}=\unicode[STIX]{x1D735}p$ the current has the form $\unicode[STIX]{x1D707}_{0}\boldsymbol{j}=K(\unicode[STIX]{x1D712},t)\boldsymbol{B}$, so that, using Ampère’s law $\unicode[STIX]{x1D707}_{0}\boldsymbol{j}=\unicode[STIX]{x1D735}\times \boldsymbol{B}$, we find

(2.7)$$\begin{eqnarray}\unicode[STIX]{x1D707}_{0}\boldsymbol{j}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}=KF/R^{2}=(\unicode[STIX]{x1D735}\times \boldsymbol{B})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{B}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D711})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(R^{-2}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}).\end{eqnarray}$$

Upon flux-surface averaging we then have

(2.8)$$\begin{eqnarray}K=\frac{1}{F\langle R^{-2}\rangle }\left\langle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}}{R^{2}}\right)\right\rangle =\frac{1}{F\langle R^{-2}\rangle V^{\prime }}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}V^{\prime }\left\langle \left|\frac{\unicode[STIX]{x1D735}r}{R}\right|^{2}\right\rangle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D712}}{\unicode[STIX]{x2202}r},\end{eqnarray}$$

where we have introduced a coordinate $r$ labelling flux surfaces, the volume $V(r)$ enclosed by such surfaces and the prime denotes a derivative with respect to $r$.

We now specialize further to the limit of large aspect ratio, in which the current density $j_{\Vert }/R=\boldsymbol{j}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}$ and the electric field $E\simeq R^{-1}\unicode[STIX]{x2202}\unicode[STIX]{x1D712}/\unicode[STIX]{x2202}t$ are approximately constant over each magnetic surface. Taking a time derivative of (2.7) and using (2.6), we find that these quantities satisfy the following equation:

(2.9)$$\begin{eqnarray}\unicode[STIX]{x1D707}_{0}\frac{\unicode[STIX]{x2202}j_{\Vert }}{\unicode[STIX]{x2202}t}=\frac{1}{V^{\prime }}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}V^{\prime }\langle |\unicode[STIX]{x1D735}r|^{2}\rangle \frac{\unicode[STIX]{x2202}E}{\unicode[STIX]{x2202}r}.\end{eqnarray}$$

If the magnetic surfaces have elliptical cross-section with elongation $\unicode[STIX]{x1D705}$ (which can depend on radius) and we let $r$ denote the radius along the minor axis, then $V(r)=2\unicode[STIX]{x03C0}^{2}R\unicode[STIX]{x1D705}r^{2}$ so that $V^{\prime }(r)=2\unicode[STIX]{x03C0}^{2}R(\unicode[STIX]{x1D705}^{\prime }r^{2}+2\unicode[STIX]{x1D705}r)$. Since $r^{2}=x^{2}+y^{2}/\unicode[STIX]{x1D705}^{2}$ we have $r\unicode[STIX]{x1D735}r=x\unicode[STIX]{x1D735}x+y\unicode[STIX]{x1D735}y/\unicode[STIX]{x1D705}^{2}-(y^{2}\unicode[STIX]{x1D705}^{\prime }/\unicode[STIX]{x1D705}^{3})\unicode[STIX]{x1D735}r$, which in terms of the polar angle $\unicode[STIX]{x1D703}$ ($\tan \unicode[STIX]{x1D703}=\unicode[STIX]{x1D705}^{-1}y/x$) takes the form

(2.10)$$\begin{eqnarray}|\unicode[STIX]{x1D735}r|^{2}=\frac{\cos ^{2}\unicode[STIX]{x1D703}+\unicode[STIX]{x1D705}^{-2}\sin ^{2}\unicode[STIX]{x1D703}}{\left(1+{\displaystyle \frac{\unicode[STIX]{x1D705}^{\prime }r}{\unicode[STIX]{x1D705}}}\sin ^{2}\unicode[STIX]{x1D703}\right)^{2}}.\end{eqnarray}$$

We thus have, with the Jacobian $|\unicode[STIX]{x2202}(x,y)/\unicode[STIX]{x2202}(r,\unicode[STIX]{x1D703})|=\unicode[STIX]{x1D705}r+\unicode[STIX]{x1D705}^{\prime }r^{2}\sin ^{2}\unicode[STIX]{x1D703}$,

(2.11)$$\begin{eqnarray}\displaystyle \langle |\unicode[STIX]{x1D735}r|^{2}\rangle & = & \displaystyle \frac{\unicode[STIX]{x1D705}r}{\unicode[STIX]{x03C0}(2\unicode[STIX]{x1D705}r+\unicode[STIX]{x1D705}^{\prime }r^{2})}\int _{0}^{2\unicode[STIX]{x03C0}}\frac{1-(1-\unicode[STIX]{x1D705}^{-2})\sin ^{2}\unicode[STIX]{x1D703}}{1+{\displaystyle \frac{\unicode[STIX]{x1D705}^{\prime }r}{\unicode[STIX]{x1D705}}}\sin ^{2}\unicode[STIX]{x1D703}}\,\text{d}\unicode[STIX]{x1D703}\nonumber\\ \displaystyle & = & \displaystyle \frac{\unicode[STIX]{x1D705}^{-2}+\sqrt{1+{\displaystyle \frac{\unicode[STIX]{x1D705}^{\prime }r}{\unicode[STIX]{x1D705}}}}}{\left(1+{\displaystyle \frac{\unicode[STIX]{x1D705}^{\prime }r}{2\unicode[STIX]{x1D705}}}\right)\left(1+\sqrt{1+{\displaystyle \frac{\unicode[STIX]{x1D705}^{\prime }r}{\unicode[STIX]{x1D705}}}}\right)\sqrt{1+{\displaystyle \frac{\unicode[STIX]{x1D705}^{\prime }r}{\unicode[STIX]{x1D705}}}}},\end{eqnarray}$$

which can be substituted into (2.9) to give an equation for the current density evolution in the case of an arbitrary elongation. As most of the runaway generation occurs close to the magnetic axis, we can neglect finite aspect ratio effects here, which would give corrections only of order $(r/R)^{2}$.

If the elongation is constant, this equation simplifies to

(2.12)$$\begin{eqnarray}\unicode[STIX]{x1D707}_{0}\frac{\unicode[STIX]{x2202}j_{\Vert }}{\unicode[STIX]{x2202}t}=\frac{1+\unicode[STIX]{x1D705}^{-2}}{2}\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}r\frac{\unicode[STIX]{x2202}E}{\unicode[STIX]{x2202}r}.\end{eqnarray}$$

Hence, it is apparent that elongation affects the resistive diffusion of the electric field. Moreover, at fixed total plasma current and minor semi-radius of the elliptical plasma cross-section, the current density is inversely proportional to $\unicode[STIX]{x1D705}$. The induced electric field before and immediately after the thermal quench depends similarly on $\unicode[STIX]{x1D705}$. The Dreicer runaway production rate (Dreicer 1959; Connor & Hastie 1975), which is exponentially sensitive to the electric field, can therefore be reduced significantly by finite elongation.

3 Runaway generation and evolution of plasma current

The evolution of the current density is governed by a balance between the generation of the runaway electrons and the resistive diffusion of the electric field accelerating them (Eriksson et al. 2004). In the cylindrical case, a one-dimensional model of these processes is numerically evolved by the go code (Smith et al. 2006), in which the current density is assumed to be the sum of the ohmic and runaway current densities. The runaways are assumed to travel at the speed of light, so that $j_{\Vert }=\unicode[STIX]{x1D70E}E+ecn_{\text{RE}}$, where $n_{\text{RE}}$ is the number density of runaways. The go code has been used extensively for evaluating pellet- and gas-injection scenarios in JET and ITER-like plasmas (Gál et al. 2008; Fehér et al. 2011; Hollmann et al. 2015; Hesslow et al. 2019a), and for interpretative modelling of experiments, e.g. the effect of the wall material on runaway electron beam formation in the JET tokamak (Papp et al. 2013).

As the resistivity increases in connection with the thermal quench, an electric field is induced, which gives rise to a runaway seed population by velocity space diffusion into the runaway region due to small angle collisions (Dreicer 1959). To determine the Dreicer runaway growth rate accurately, we use a neural network1 (Hesslow et al. 2019b) trained on a large number of kinetic simulations by the code kinetic equation solver (Landreman, Stahl & Fülöp 2014). In the case of fully ionized plasmas and constant Coulomb logarithm, the primary runaway growth rate given by the neural network agrees with the analytical formulas for the runaway growth rate of Connor & Hastie (1975). However, due to the velocity dependence of the Coulomb logarithm, in certain regions of parameter space, in particular for low temperatures and electric fields, the growth rate can significantly differ from the analytical formulas, even in fully ionized plasmas (Hesslow et al. 2019b). This leads to substantial changes in the final runaway current, as we will see in the next section.

The primary effect of the Dreicer mechanism is to generate a ‘seed’ of runaways which is amplified by the avalanche, but there are also other ways in which seeding occurs. For instance, tritium decay produces a seed which can be modelled as (Martín-Solís, Loarte & Lehnen 2017)

(3.1)$$\begin{eqnarray}\left(\frac{\unicode[STIX]{x2202}n_{\text{RE}}}{\unicode[STIX]{x2202}t}\right)^{\text{tritium}}=\ln (2)\frac{n_{\text{T}}}{\unicode[STIX]{x1D70F}_{\text{T}}}f(W_{\text{crit}}),\end{eqnarray}$$

where $n_{\text{T}}$ is the tritium density, $\unicode[STIX]{x1D70F}_{\text{T}}\approx 4500$ days is the half-life of tritium and $f(W_{\text{crit}})$ is the fraction of the electron spectrum of the tritium decay above the critical runaway energy $W_{\text{crit}}$. If we consider the emitted electron as a free particle, i.e. neglect the effect of the Coulomb interaction between the electron and the tritium nucleus on the spectrum, the tritium energy spectrum can be shown to fulfil (Martin & Shaw 2019) $I(W)\propto (Q-W)^{2}\sqrt{W}$, where $Q=18.6~\text{keV}$. The fraction of the electron spectrum from tritium decay that lies within the runaway region can then be calculated analytically as

(3.2)$$\begin{eqnarray}f(W_{\text{crit}})=\frac{\displaystyle \int _{W_{\text{crit}}}^{Q}I(W)\,\text{d}W}{\displaystyle \int _{\text{0}}^{Q}I(W)\,\text{d}W}=1-\frac{35}{8}\left(\frac{W_{\text{crit}}}{Q}\right)^{3/2}+\frac{21}{4}\left(\frac{W_{\text{crit}}}{Q}\right)^{5/2}-\frac{15}{8}\left(\frac{W_{\text{crit}}}{Q}\right)^{7/2}\!,\end{eqnarray}$$

where $W_{\text{crit}}$ is the critical runaway energy.

The seed runaways are amplified due to close collisions. In fully ionized plasmas, the avalanche growth rate is given by (Rosenbluth & Putvinski 1997)

(3.3)$$\begin{eqnarray}\displaystyle \left(\frac{\unicode[STIX]{x2202}n_{\text{RE}}}{\unicode[STIX]{x2202}t}\right)^{\text{aval}} & = & \displaystyle \frac{e(E_{\Vert }-E_{c})n_{\text{RE}}}{m_{\text{e}}c\ln \unicode[STIX]{x1D6EC}_{\text{c}}}\sqrt{\frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D716}}}{3(Z_{\text{eff}}+5)}}\nonumber\\ \displaystyle & & \displaystyle \times \left(1-\frac{E_{c}}{E_{\Vert }}+\frac{4\unicode[STIX]{x03C0}(Z_{\text{eff}}+1)^{2}}{3\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D716}}(Z_{\text{eff}}+5)(E_{\Vert }^{2}/E_{c}^{2}+4/\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D716}}^{2}-1)}\right)^{-1/2},\end{eqnarray}$$

where $E_{\text{c}}=m_{\text{e}}c/(e\unicode[STIX]{x1D70F}_{\text{c}})$ is the Connor–Hastie critical electric field, $\unicode[STIX]{x1D70F}_{\text{c}}=4\unicode[STIX]{x03C0}\unicode[STIX]{x1D716}_{0}^{2}m_{\text{e}}^{2}c^{3}/(n_{\text{e}}e^{4}\ln \unicode[STIX]{x1D6EC}_{\text{c}})$ is the relativistic collision time, $\ln \unicode[STIX]{x1D6EC}_{\text{c}}\approx 14.6+0.5\ln (T_{\text{eV}}/n_{\text{e20}})$ is the Coulomb logarithm for relativistic electrons, with $T_{\text{eV}}$ the electron temperature in electronvolts, $n_{\text{e20}}$ the density of the background electrons in units of $10^{20}~\text{m}^{-3}$, $\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D716}}=(1+1.46\unicode[STIX]{x1D716}^{1/2}+1.72\unicode[STIX]{x1D716})^{-1}$ is the neoclassical factor and $\unicode[STIX]{x1D716}=r/R$ denotes the inverse aspect ratio. The effect of elongation on the neoclassical factor $\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D716}}$ is negligible. Note that, in the presence of partially ionized impurities, the avalanche growth rate will no longer be directly proportional to the electric field, and the multiplication factor will instead become sensitive to the details of the electric field evolution (Martín-Solís, Loarte & Lehnen 2015; Hesslow et al. 2019a); however, in fully ionized plasmas which we consider here, these effects can be ignored.

4 Dependence of avalanche multiplication on plasma elongation

When primary runaway generation is negligible, the plasma current evolution is governed by the diffusion of the electric field and runaway avalanche multiplication according to (2.9) and (3.3),

(4.1)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}_{0}\frac{\unicode[STIX]{x2202}j_{\Vert }}{\unicode[STIX]{x2202}t}=\frac{1}{V^{\prime }}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}V^{\prime }\langle |\unicode[STIX]{x1D735}r|^{2}\rangle \frac{\unicode[STIX]{x2202}E_{\Vert }}{\unicode[STIX]{x2202}r} & \displaystyle\end{eqnarray}$$
(4.2)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\ln n_{\text{RE}}}{\unicode[STIX]{x2202}t}=\frac{E_{\Vert }}{E_{\text{c}}\unicode[STIX]{x1D70F}_{\text{a}}}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70F}_{\text{a}}=\unicode[STIX]{x1D70F}_{\text{c}}\ln \unicode[STIX]{x1D6EC}_{\text{c}}\sqrt{5+Z_{\text{e ff }}}$ and we have assumed $E_{\Vert }\gg E_{\text{c}}$.

If the density, effective charge and the Coulomb logarithm are constant in time, we can integrate the diffusion equation (4.1) from $t=0$ to $t=\infty$ when the entire current is carried by runaways,

(4.3)$$\begin{eqnarray}-\unicode[STIX]{x1D707}_{0}[j_{\text{0}}-j_{\text{RE}}]=\frac{1}{V^{\prime }}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}V^{\prime }\langle |\unicode[STIX]{x1D735}r|^{2}\rangle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left[E_{\text{c}}\unicode[STIX]{x1D70F}_{\text{a}}\ln \frac{j_{\text{RE}}(r)}{ecn_{\text{seed}}(r)}\right].\end{eqnarray}$$

To find the maximum possible avalanche multiplication in the trace runaway limit, $j_{\text{RE}}\ll j_{\text{0}}$, we integrate in radius from $0$ to $r$ and obtain

(4.4)$$\begin{eqnarray}2\unicode[STIX]{x03C0}R\unicode[STIX]{x1D707}_{0}I_{\text{0}}(r)=-V^{\prime }(r)\langle |\unicode[STIX]{x1D735}r|^{2}\rangle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left[E_{\text{c}}\unicode[STIX]{x1D70F}_{\text{a}}\ln \frac{j_{\text{RE}}(r)}{ecn_{\text{seed}}(r)}\right],\end{eqnarray}$$

where $I_{\text{0}}=(2\unicode[STIX]{x03C0}R)^{-1}\int _{0}^{r}\,\text{d}r^{\prime }V^{\prime }(r^{\prime })j(r^{\prime })$ denotes the total initial current enclosed by a flux surface of minor radius $r$. Integrating again from $r$ to $a$, where we assume a perfectly conducting wall, we obtain the avalanche multiplication factor that accounts for the radial diffusion of the electric field,

(4.5)$$\begin{eqnarray}\ln \frac{j_{\text{RE}}(r)}{ecn_{\text{seed}}(r)}=\frac{\unicode[STIX]{x1D707}_{0}}{E_{\text{c}}(r)\unicode[STIX]{x1D70F}_{\text{a}}(r)}\int _{r}^{a}\frac{2\unicode[STIX]{x03C0}RI_{\text{0}}(r^{\prime })\,\text{d}r^{\prime }}{V^{\prime }(r^{\prime })\langle |\unicode[STIX]{x1D735}r|^{2}\rangle (r^{\prime })}.\end{eqnarray}$$

In a geometry with constant elongation, $\langle |\unicode[STIX]{x1D735}r|^{2}\rangle =(1+\unicode[STIX]{x1D705}^{-2})/2$ and $V^{\prime }(r)=4\unicode[STIX]{x03C0}^{2}\unicode[STIX]{x1D705}rR$, and assuming that the plasma current $I$ is independent of the elongation, equation (4.5) gives

(4.6)$$\begin{eqnarray}\ln \frac{j_{\text{RE}}}{ecn_{\text{seed}}}=\frac{\unicode[STIX]{x1D707}_{0}G(\unicode[STIX]{x1D705})}{2\unicode[STIX]{x03C0}E_{\text{c}}\unicode[STIX]{x1D70F}_{a}}\int _{r}^{a}\frac{\text{d}r^{\prime }}{r^{\prime }}I_{\text{0}}(r^{\prime }),\end{eqnarray}$$

where $G(\unicode[STIX]{x1D705})=2/(\unicode[STIX]{x1D705}+\unicode[STIX]{x1D705}^{-1})$. Thus the avalanche generated runaway current will be reduced by a factor of $\exp [N_{\text{exp}}^{\unicode[STIX]{x1D705}=1}(1-G(\unicode[STIX]{x1D705}))]$, where $N_{\text{exp}}^{\unicode[STIX]{x1D705}=1}$ is the number of exponentiations for $\unicode[STIX]{x1D705}=1$. In the case of $\unicode[STIX]{x1D705}=1.45$ used in the examples in the next section, we have $G(1.45)=0.93$. With a representative value of $N_{\text{exp}}\approx 50$, for ITER-like parameters in a fully ionized plasma (Rosenbluth & Putvinski 1997), this would correspond to a reduction factor of approximately 26. However, this argument applies only in the case when the final runaway current is much smaller than the initial one. For higher runaway currents, the effect of the runaway current on the electric field will cause the runaway current to saturate and hence the number of exponentiations can be much lower.

5 Numerical results for high-current devices

To illustrate the effect of elongation, in the following, we present numerical solutions of (2.12), with the runaway growth rate given by the sum of the primary (Dreicer+tritium decay) and avalanche growth rates, for parameters characteristic of a SPARC V0 and an ITER discharge. We take the temperature to decay exponentially as $T_{\text{e}}(t,x)=T_{\text{f}}(x)+[T_{\text{0}}(x)-T_{\text{f}}(x)]\text{e}^{-t/t_{\text{0}}}$, where $T_{\text{0}}$ and $T_{\text{f}}$ are the initial and final electron temperatures, respectively, and $t_{\text{0}}$ is the thermal quench (TQ) time. For simplicity, we consider pure plasmas and neglect hot-tail generation and runaways produced by Compton scattering of $\unicode[STIX]{x1D6FE}$-rays due to the activated wall in the nuclear phase of operation. In view of all these assumptions, the results can only be used as an illustration of the effect of elongation, and not to draw conclusions on the final runaway current in any future tokamak.

Figure 1. Radial profiles of plasma parameters for SPARC V0. (a) Initial current density. (b) Electron density. (c) Electron temperature. (d) Elongation.

Figure 2. Plasma current and electric field evolution in a simulated SPARC V0 thermal quench. (a) Total plasma current as function of time. Dotted lines correspond to circular plasma ($\unicode[STIX]{x1D705}=1$), and dashed lines are for ($\unicode[STIX]{x1D705}=1.45$). (b) Contour plot of the current conversion $I_{\text{RE}}/I_{\text{tot}}$ as a function of TQ time $T_{\text{0}}$ and final temperature $T_{f}$ for $\unicode[STIX]{x1D705}=1$ (dashed) and $\unicode[STIX]{x1D705}=1.45$ (solid). (c,d) Electric field and current density evolution for circular (dashed) and elongated (solid) plasmas. The parameters are $t_{\text{0}}=1~\text{ms}$ and $T_{\text{f}}=20~\text{eV}$, except in (b), where they are varied.

Figure 3. Plasma current and electric field evolution in a simulated ITER thermal quench. (a) Total plasma current as a function of time for circular (dashed) and elongated (solid) plasmas. (b) Contour plot of the current conversion $I_{\text{RE}}/I_{\text{tot}}$ as a function of TQ time $t_{\text{0}}$ and final temperature $T_{f}$ for $\unicode[STIX]{x1D705}=1.45$. (c,d) Electric field and current density evolution. The parameters are $t_{\text{0}}=1~\text{ms}$ and $T_{f}=20~\text{eV}$ except in (b) where they vary.

For SPARC V0, the initial plasma current is $I_{\text{0}}=7.5~\text{MA}$, major radius $R=1.65~\text{m}$ and minor radius $a=0.5~\text{m}$. The initial current density, electron density, temperature and elongation profiles are shown in figure 1. Profiles were obtained from predictive modelling using the transp code (Breslau et al. 2018), which takes into account sources, sinks and transport projected to be present in SPARC V0. In the simulations presented below we keep the total plasma current constant when changing the elongation, which means that the local current density is rescaled correspondingly.

After a disruption, the final temperature profile is usually flatter than the initial one, and we will assume it to be constant $T_{\text{f}}=20~\text{e}V$. The final density is normally larger than before the disruption, due to an influx of impurities from the wall or intentional injection of gas to mitigate the effects of the disruption. Here, for simplicity we take the density to be constant in time with the same radial profile as the initial density.

Figure 2 shows the electric field and current evolution in a simulated SPARC V0 thermal quench, for both elongated and circular plasma shapes. As figure 1(d) shows, the elongation varies radially, but our simulations show that what matters for runaway generation is the value of the elongation in the central part of the plasma. For SPARC V0 we find that the runaway current evolution is the same when we take into account the radial dependence of the elongation as when we take the value $\unicode[STIX]{x1D705}=1.45$. As the difference to the radially varying case is insignificant, figure 2 only shows the electric field and current evolution for a constant value of $\unicode[STIX]{x1D705}=1.45$.

Figure 2(a) shows that, if the plasma is elongated, only a negligible part of the initial plasma current is converted to a runaway current, compared to the circular case. Figure 2(b) shows the current conversion as a function of TQ time and final temperature in the elongated case. Clearly, runaway production is not significant unless the cooling is extremely rapid and the final temperature reaches less than a few eV.

The simulations show that the main reason for the reduction in the current conversion in elongated plasmas is the lower maximum electric field compared to circular plasmas, as shown in figure 2(c). The change in electric field significantly affects the Dreicer generation. The avalanche multiplication factor is also reduced for $\unicode[STIX]{x1D705}>1$, as shown in the previous section. The maximum of the electric field and consequently the runaway production is mainly on axis, as shown in figure 2(d).

Including tritium seed generation leads to a runaway current conversion of 4.4 % in the cylindrical case and 1.2 % in the elongated case for TQ time $t_{\text{0}}=1~\text{ms}$ and final temperature $T_{\text{f}}=20~\text{eV}$. The corresponding numbers without tritium are 4.1 % in the cylindrical case and $10^{-6}$ in the elongated case. For $T_{\text{f}}=20~\text{eV}$, Dreicer and tritium runaway seeds are of the same order of magnitude in the cylindrical case, but the tritium seed dominates in the elongated case.

For the ITER scenario we consider, the initial plasma current is $I_{\text{0}}=15~\text{MA}$ and the major and minor radii are $R=6~\text{m}$ and $a=2~\text{m}$, respectively. The pre-disruption average temperature is $\langle T_{\text{e}}\rangle =10~\text{keV}$ and density $\langle n_{\text{e}}\rangle =10^{20}~\text{m}^{-3}$. The temperature profile is taken as $T_{\text{e}}=T_{0}[1-(r/a)^{2}]$, with $T_{0}=2\langle T_{\text{e}}\rangle$, and the electron density $n_{\text{e}}$ profile is assumed to be flat. The initial current density profile is assumed to be $j(r)=j_{\text{0}}[1-(r/a)^{0.41}]$, corresponding to an internal inductance of $l_{\text{i}}=0.7$; $j_{\text{0}}$ is a normalization parameter chosen so that the total plasma current integrates to $I_{p}$. For $\unicode[STIX]{x1D705}=1$ it is $j_{\text{0}}=1.69~\text{MA}~\text{m}^{-2}$ and for $\unicode[STIX]{x1D705}=1.45$, $j_{\text{0}}=1.16~\text{MA}~\text{m}^{-2}$. These parameters and initial profiles are similar to the ones used by Martín-Solís et al. (2017), except for the effect of elongation, which was not considered there.

Figure 3 shows the current and electric field evolution in a simulated ITER disruption, comparing a circular and an elongated plasma, with constant elongation $\unicode[STIX]{x1D705}=1.45$. In the considered case, the maximum electric field is much lower in ITER than in SPARC. Therefore, the Dreicer part of the runaway seed is negligibly small, and the tritium seed dominates by orders of magnitude. This was noted also in previous work, e.g. by Martín-Solís et al. (2017).

The tritium seed generation is primarily determined by the time interval during which the electric field is high enough for the critical energy for runaway generation $W_{\text{crit}}$ to be lower than the maximum energy of the emitted electrons during tritium decay $Q=18.6~\text{keV}$. This interval increases with elongation, and consequently the tritium seed slightly increases with $\unicode[STIX]{x1D705}$. However, the final runaway current is still reduced, but only marginally. The reason for the reduction is that the avalanche multiplication is weaker in the elongated case. The sensitivity of the current conversion to the TQ time and final temperature is shown in figure 3(b).

We do not consider the effect of shaping on the MHD stability of the discharge, which might give rise to radial transport of energetic runaways, or any other loss processes. We have also ignored several processes that would lead to higher runaway currents. Perhaps one of the most important of these is hot-tail generation of runaways, which is expected to be significant in large tokamak disruptions (Chiu et al. 1998; Helander et al. 2004; Aleynikov & Breizman 2017; Martín-Solís et al. 2017). Hot-tail generation is very sensitive to the details of the cooling process following the magnetic reconnection, which is not well understood (Breizman & Aleynikov 2017). In addition, as the hot-tail runaways are produced in the early phase of the TQ, their transport is likely to be significantly affected by the high level of magnetic fluctuations following the magnetic reconnection. Therefore, it is difficult to obtain accurate predictions regarding the effect of hot-tail generation.

An additional limitation of the analysis above is the assumption of a pre-described exponentially decaying temperature evolution, with a flat final temperature profile. In a more realistic scenario, the temperature is determined by a balance between ohmic heating and impurity radiation, and it will have an important effect on the evolving current density profile. To investigate whether or not the effect of elongation changes if we use a different temperature profile, we have simulated a case where the radial temperature profile is kept constant (i.e. the same as the initial). We find that the effect of elongation is not affected by this change. For ITER, the change in the final runaway current is negligible; the only effect is that the current quench becomes somewhat shorter and the electric field larger, but during a shorter time, leading to the same final runaway current. For SPARC, where Dreicer generation dominates, the difference in the final current is somewhat larger but the trend of how the elongation reduces the current is the same.

Currently envisaged disruption mitigation methods involve injection of massive amounts of material, and that will also change the current dynamics substantially, and may lead to higher runaway currents (Hesslow et al. 2019a). On the other hand, in connection with material injection, the injected density required to raise the critical field $E_{c}$ for runaway generation (or the threshold field for tritium seed generation) above the maximum induced field will be lower by a factor of $\unicode[STIX]{x1D705}$ in elongated plasmas.

6 Conclusions

We show that elongated plasmas are less prone to runaway electron generation in tokamak disruptions. Since the current density is approximately inversely proportional to the vertical elongation $\unicode[STIX]{x1D705}$, the maximum induced electric field is reduced by a similar factor, which has a significant effect on the primary runaway generation, which is exponentially sensitive to the electric field. In addition, shaping reduces the maximum avalanche gain by a factor of $2/(\unicode[STIX]{x1D705}+\unicode[STIX]{x1D705}^{-1})$. Numerical solution of the coupled equations of runaway generation and resistive diffusion of electric field in simulated disruptions in high-current devices show that the final runaway current is expected to be reduced considerably in tokamaks where the primary runaway generation is dominated by the Dreicer process. When the primary generation is dominated by other processes, such as tritium decay, we expect the elongation to cause only a marginal reduction of the final runaway current.


The authors are grateful to S. Newton, M. Hoppe, L. Unnerfelt, A. Tinguely and G. Papp for fruitful discussions. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 647121. The work was also supported by the Swedish Research Council (Dnr. 2018-03911), the EUROfusion – Theory and Advanced Simulation Coordination (E-TASC) and was carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 and 2019–2020 under grant agreement no. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.


Aleynikov, P. & Breizman, B. N.2017 Generation of runaway electrons during the thermal quench in tokamaks. Nucl. Fusion 57 (4), 046009.
Boozer, A. H.2018 Pivotal issues on relativistic electrons in ITER. Nucl. Fusion 58 (3), 036006.
Breizman, B. & Aleynikov, P.2017 Kinetics of relativistic runaway electrons. Nucl. Fusion 57 (12), 125002.
Breizman, B. N., Aleynikov, P., Hollmann, E. M. & Lehnen, M.2019 Physics of runaway electrons in tokamaks. Nucl. Fusion 59 (8), 083001.
Breslau, J., Gorelenkova, M., Poli, F., Sachdev, J. & Yuan, X.TRANSP. (Computer Software) USDOE Office of Science (SC), Fusion Energy Sciences (FES) (SC-24). 27 Jun. 2018. doi:10.11578/dc.20180627.4.
Chiu, S., Rosenbluth, M., Harvey, R. & Chan, V.1998 Fokker–Planck simulations mylb of knock-on electron runaway avalanche and bursts in tokamaks. Nucl. Fusion 38 (11), 17111721.
Connor, J. & Hastie, R.1975 Relativistic limitations on runaway electrons. Nucl. Fusion 15, 415.
Dreicer, H.1959 Electron and ion runaway in a fully ionized gas. I. Phys. Rev. 115, 238.
Eriksson, L.-G., Helander, P., Andersson, F., Anderson, D. & Lisak, M.2004 Current dynamics during disruptions in large tokamaks. Phys. Rev. Lett. 92, 205004.
Fehér, T., Smith, H. M., Fülöp, T. & Gál, K.2011 Simulation of runaway electron generation during plasma shutdown by impurity injection in ITER. Plasma Phys. Control. Fusion 53 (3), 035014.
Gál, K., Fehér, T., Smith, H. M., Fülöp, T. & Helander, P.2008 Runaway electron generation during plasma shutdown by killer pellet injection. Plasma Phys. Control. Fusion 50, 055006.
Greenwald, M., Whyte, D., Bonoli, P., Hartwig, Z., Irby, J., Labombard, B., Marmar, E., Minervini, J., Takayasu, M., Terry, J. et al. 2018 The high-field path to practical fusion energy.
Helander, P. & Sigmar, D.2005 Collisional Transport in Magnetized Plasmas. Cambridge University Press.
Helander, P., Smith, H. M., Fülöp, T. & Eriksson, L. G.2004 Electron kinetics in a cooling plasma. Phys. Plasmas 11, 5704.
Hesslow, L., Embréus, O., Vallhagen, O. & Fülöp, T.2019a Influence of massive material injection on avalanche runaway generation during tokamak disruptions. Nucl. Fusion 59 (8), 084004.
Hesslow, L., Unnerfelt, L., Vallhagen, O., Embreus, O., Hoppe, M., Papp, G. & Fülöp, T.2019b Evaluation of the Dreicer runaway generation rate in the presence of high-$Z$ impurities using a neural network. J. Plasma Phys. 85 (6), 475850601.
Hollmann, E., Austin, M., Boedo, J., Brooks, N., Commaux, N., Eidietis, N., Humphreys, D., Izzo, V., James, A., Jernigan, T. et al. Control and dissipation of runaway electron beams created during rapid shutdown experiments in DIII-D. Nucl. Fusion 53 (8), 083004.
Hollmann, E. M., Aleynikov, P. B., Fülöp, T., Humphreys, D. A., Izzo, V. A., Lehnen, M., Lukash, V. E., Papp, G., Pautasso, G., Saint-Laurent, F. et al. 2015 Status of research toward the ITER disruption mitigation system. Phys. Plasmas 22 (2), 021802.
Izzo, V., Hollmann, E., James, A., Yu, J., Humphreys, D., Lao, L., Parks, P., Sieck, P., Wesley, J., Granetz, R. et al. 2011 Runaway electron confinement modelling for rapid shutdown scenarios in DIII-D, Alcator C-Mod and ITER. Nucl. Fusion 51 (6), 063032.
Izzo, V. A., Humphreys, D. A. & Kornbluth, M.2012 Analysis of shot-to-shot variability in post-disruption runaway electron currents for diverted DIII-D discharges. Plasma Phys. Control. Fusion 54 (9), 095002.
Jayakumar, R., Fleischmann, H. & Zweben, S.1993 Collisional avalanche exponentiation of runaway electrons in electrified plasmas. Phys. Lett. A 172, 447451.
Landreman, M., Stahl, A. & Fülöp, T.2014 Numerical calculation of the runaway electron distribution function and associated synchrotron emission. Comput. Phys. Commun. 185 (3), 847.
Martin, B. R. & Shaw, G.2019 Nuclear and Particle Physics: An Introduction. John Wiley & Sons, Inc.
Martín-Solís, J. R., Loarte, A. & Lehnen, M.2015 Runaway electron dynamics in tokamak plasmas with high impurity content. Phys. Plasmas 22, 092512.
Martín-Solís, J. R., Loarte, A. & Lehnen, M.2017 Formation and termination of runaway beams in ITER disruptions. Nucl. Fusion 57 (6), 066025.
Papp, G., Fülöp, T., Fehér, T., de Vries, P., Riccardo, V., Reux, C., Lehnen, M., Kiptily, V., Plyusnin, V., Alper, B.& JET EFDA contributors 2013 The effect of ITER-like wall on runaway electron generation in JET. Nucl. Fusion 53 (12), 123017.
Reux, C., Plyusnin, V., Alper, B., Alves, D., Bazylev, B., Belonohy, E., Boboc, A., Brezinsek, S., Coffey, I., Decker, J. et al. 2015 Runaway electron beam generation and mitigation during disruptions at JET-ILW. Nucl. Fusion 55 (9), 093013.
Rosenbluth, M. & Putvinski, S.1997 Theory for avalanche of runaway electrons in tokamaks. Nucl. Fusion 37, 13551362.
Smith, H., Helander, P., Eriksson, L.-G., Anderson, D., Lisak, M. & Andersson, F.2006 Runaway electrons and the evolution of the plasma current in tokamak disruptions. Phys. Plasmas 13 (10), 102502.
Wilson, C. T. R.1925 The acceleration of $\unicode[STIX]{x1D6FD}$-particles in strong electric fields such as those of thunderclouds. Math. Proc. Cambridge Philos. Soc. 22, 534.

1 The neural network is available at