Effect of plasma elongation on current dynamics during tokamak disruptions

Abstract

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 network¹ (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.