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 Reference Wilson1925; Dreicer Reference Dreicer1959). 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 Reference Jayakumar, Fleischmann and Zweben1993).
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 Reference Rosenbluth and Putvinski1997). 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. Reference Greenwald, Whyte, Bonoli, Hartwig, Irby, Labombard, Marmar, Minervini, Takayasu and Terry2018), 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 Reference Izzo, Humphreys and Kornbluth2012; Hollmann et al. Reference Hollmann, Austin, Boedo, Brooks, Commaux, Eidietis, Humphreys, Izzo, James and Jernigan2013; Reux et al. Reference Reux, Plyusnin, Alper, Alves, Bazylev, Belonohy, Boboc, Brezinsek, Coffey and Decker2015; Breizman et al. Reference Breizman, Aleynikov, Hollmann and Lehnen2019). 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. Reference Izzo, Hollmann, James, Yu, Humphreys, Lao, Parks, Sieck, Wesley and Granetz2011). 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
$$\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 (Reference Boozer2018), we find that the electric field can be written as
$$\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
$$\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
$$\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 Reference Helander and Sigmar2005)
$$\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
$$\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
$$\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
$$\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:
$$\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
$$\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}$,
$$\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
$$\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 Reference Dreicer1959; Connor & Hastie Reference Connor and Hastie1975), 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. Reference Eriksson, Helander, Andersson, Anderson and Lisak2004). In the cylindrical case, a one-dimensional model of these processes is numerically evolved by the go code (Smith et al. Reference Smith, Helander, Eriksson, Anderson, Lisak and Andersson2006), 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. Reference Gál, Fehér, Smith, Fülöp and Helander2008; Fehér et al. Reference Fehér, Smith, Fülöp and Gál2011; Hollmann et al. Reference Hollmann, Aleynikov, Fülöp, Humphreys, Izzo, Lehnen, Lukash, Papp, Pautasso and Saint-Laurent2015; Hesslow et al. Reference Hesslow, Embréus, Vallhagen and Fülöp2019a), 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. Reference Papp, Fülöp, Fehér, de Vries, Riccardo, Reux, Lehnen, Kiptily, Plyusnin and Alper2013).
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 Reference Dreicer1959). To determine the Dreicer runaway growth rate accurately, we use a neural networkFootnote 1 (Hesslow et al. Reference Hesslow, Unnerfelt, Vallhagen, Embreus, Hoppe, Papp and Fülöp2019b) trained on a large number of kinetic simulations by the code kinetic equation solver (Landreman, Stahl & Fülöp Reference Landreman, Stahl and Fülöp2014). 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 (Reference Connor and Hastie1975). 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. Reference Hesslow, Unnerfelt, Vallhagen, Embreus, Hoppe, Papp and Fülöp2019b). 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 Reference Martín-Solís, Loarte and Lehnen2017)
$$\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 Reference Martin and Shaw2019) 
$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
$$\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 Reference Rosenbluth and Putvinski1997)
$$\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 Reference Martín-Solís, Loarte and Lehnen2015; Hesslow et al. Reference Hesslow, Embréus, Vallhagen and Fülöp2019a); 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),
$$\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}$$
$$\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,
$$\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
$$\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,
$$\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
$$\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 Reference Rosenbluth and Putvinski1997), 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. Reference Breslau, Gorelenkova, Poli, Sachdev and Yuan2018), 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. (Reference Martín-Solís, Loarte and Lehnen2017), 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. (Reference Martín-Solís, Loarte and Lehnen2017).
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. Reference Chiu, Rosenbluth, Harvey and Chan1998; Helander et al. Reference Helander, Smith, Fülöp and Eriksson2004; Aleynikov & Breizman Reference Aleynikov and Breizman2017; Martín-Solís et al. Reference Martín-Solís, Loarte and Lehnen2017). Hot-tail generation is very sensitive to the details of the cooling process following the magnetic reconnection, which is not well understood (Breizman & Aleynikov Reference Breizman and Aleynikov2017). 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. Reference Hesslow, Embréus, Vallhagen and Fülöp2019a). 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.
Acknowledgements
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.
Open access
