DREAM: a fluid-kinetic framework for tokamak disruption runaway electron simulations

Avoidance of the harmful effects of runaway electrons (REs) in plasma-terminating disruptions is pivotal in the design of safety systems for magnetic fusion devices. Here, we describe a computationally efficient numerical tool, that allows for self-consistent simulations of plasma cooling and associated RE dynamics during disruptions. It solves flux-surface averaged transport equations for the plasma density, temperature and poloidal flux, using a bounce-averaged kinetic equation to self-consistently provide the electron current, heat, density and RE evolution, as well as the electron distribution function. As an example, we consider disruption scenarios with material injection and compare the electron dynamics resolved with different levels of complexity, from fully kinetic to fluid modes.


Introduction
Disruptions of tokamak plasmas involve a partial loss of magnetic confinement and a sudden cooling of the plasma [1]. This thermal quench leads to an increase in the plasma resistivity, causing the plasma current to decay over a period termed the current quench. The toroidal current cannot change significantly on the short thermal quench timescale, and therefore an inductive electric field is produced that can lead to electron runaway. The main runaway generation processes are the Dreicer [2,3], hot-tail [4,5,6] and avalanche [7,8,9] mechanisms. In the nuclear phase of operations, these are complemented with electrons generated by the beta decay of tritium and Compton scattering of γ-rays emitted by the activated wall.
During the current quench, a large part of the plasma current can be converted to a beam of energetic electrons which has the potential to cause severe damage to plasma facing components [10,11]. As the runaway generation is exponentially sensitive to the plasma current [8], this problem is ex-pected to be unacceptable in future high current tokamaks, such as ITER [1] and SPARC [12]. Development of control methods for avoidance or mitigation of disruptions is therefore of critical importance and urgency in fusion physics.
The most discussed disruption mitigation method is material injection. However, in certain cases the injection of impurities and the associated radiative cooling can lead to the generation of even larger runaway currents [13]. There are a large number of degrees of freedom associated with proposed mitigation methods [14]. However, only a small part of the parameter space is accessible to existing experiments and the extrapolation of the results of existing experiments to next-generation tokamaks is not straightforward [15,16]. Therefore theoretical modelling and computationally efficient numerical simulations of the disruption and associated runaway electron generation are essential.
Nonlinear magnetohydrodynamic (MHD) simulation tools, such as Jorek [17,18] and Nimrod [19] have many of the necessary components to simulate mitigated disruptions. However, they so far only allow studies of the dynamics of runaway electrons as test particles [20,21] or via a fluid model [22,23,24], and are computationally expensive. To self-consistently simulate disruption dynamics, an integrated tool that can simulate situations when relativistic electrons comprise a significant part of the electron distribution is required.
Simplified fluid codes such as go [13] and the 1.5dimensional transport toolkit astra-strahl [25] include kinetically benchmarked models for Dreicer and avalanche runaway electron generation, but have simplified models for hot-tail generation [26]. Fully kinetic tools modelling the momentum-space dynamics of relativistic electrons, such as code [27,28], norse [29], or the bounce-averaged codes luke [30] and cql3d [4,31] are suitable for capturing the hot-tail generation, however, to couple these codes to a self-consistent model of the global disruption dynamics leads to prohibitively expensive simulations.
In this paper we describe a new integrated tool for self-consistently simulating the evolution of temperature, poloidal flux, and impurity densities, along with the generation and transport of runaway electrons in tokamaks: Dream (Disruption Runaway Electron Analysis Model). The fullyimplicit tool solves a nonlinear set of coupled equations describing the evolution of temperature, den-sity, current density and electric field, as well as the full electron distribution function in arbitrary axisymmetric geometry. It employs a combination of fluid models for background plasma parameters, including the toroidal electric field, electron and ion temperatures, ion densities and charge states, as well as various models for runaway electrons, ranging from fluid to fully kinetic. The most complete drift-kinetic model includes a fully relativistic Fokker-Planck test-particle operator for electronelectron collisions, synchrotron radiation reaction, an avalanche operator, bremsstrahlung and effects of screening in a partially ionized plasma. Neglecting the field-particle part of the collision operator means that the conductivity is underestimated [32], therefore the ohmic current is amended with a conductivity correction to capture the correct Spitzer response to an electric field. A distinguishing feature of Dream is the possibility of choosing reduced kinetic modes, which allow parts of the electron phase space to be modelled kinetically, and the remainder to be described by fluid equations.
The equations that describe the electron kinetics and background plasma evolution in Dream are discussed in Sections 2-3. The numerical implementation is then outlined in Section 4, and benchmarked in Section 5 through a comparison to other numerical tools in various limits. Finally, in Section 6, we use Dream to investigate a disruption simulation in a toroidal plasma, through which we highlight the difference between the hierarchy of electron models implemented.

Electron kinetics
Electrons in Dream are primarily modelled with a bounce-averaged drift-kinetic equation of the form where f (t, z) denotes the electron distribution, z i denotes phase space coordinates and the contravariant components A m = (∂z m /∂z) · A and (in dyadic notation) D mn = (∂z m /∂z)(∂z n /∂z) : D, where A and D represent the underlying advection vector and diffusion tensor, driving electron phasespace flows. We consider the zero-orbit-width limit, wherein the Larmor radius and cross field drifts are neglected, so that electrons exactly follow magnetic field lines, and choose z as the constants of motion z = (r, p, ξ 0 ), where the flux-surface label r is the distance from the magnetic axis when the particle orbit passes the point of minimum magnetic field strength B min ; ξ 0 = B · p/(Bp)| B=Bmin is the particle pitch with respect to the magnetic field B at this point and p is the magnitude of the particle momentum. The kinetic equation (1) has been averaged over the remaining three coordinates, using the bounce average denoted with curly brackets with toroidal angle φ, gyrophase ζ and where the angle θ is an arbitrary poloidal angle which parametrizes the flux surface, assumed to be 2πperiodic. The metric √ g of the (r, θ, φ, p, ξ, ζ) coordinates and the spatial Jacobian J are defined by where ξ = sgn(ξ 0 ) 1 − (1 − ξ 2 0 )B/B min . The poloidal integral along the particle orbit is taken as dθ X(r, θ, φ, p, ξ, ζ) =          π −π dθ X(r, θ, φ, p, ξ, ζ), |ξ 0 | > ξ T θ b2 θ b1 dθ [X(r, θ, φ, p, ξ, ζ) +X(r, θ, φ, p, −ξ, ζ)], 0 where the bounce points θ b1 and θ b2 are the two separate poloidal angles defined by ξ(ξ 0 , θ b1 ) = ξ(ξ 0 , θ b2 ) = 0, and the trapped-passing boundary is denoted ξ T = 1 − B min /B max . In the positive trapped region, 0 < ξ 0 ≤ ξ T , the full contributions from both co-moving and counter-moving trapped particles are accounted for, and therefore all bounce averages are set to zero in −ξ T ≤ ξ 0 ≤ 0 in order to avoid double counting (in this region the solution satisfies f (ξ 0 ) = f (−ξ 0 )). We will also utilize the spatial flux surface average The geometry of the magnetic surfaces enters into J and into the dependence on poloidal angle of B(θ); in Dream we evolve the equations of motion in a static magnetic geometry where each flux surface is parametrized by its major radius, elongation, Shafranov shift, triangularity and a reference poloidal flux profile, as described in Appendix A. The magnetic field is represented by the mixed form with ψ the poloidal flux-labeling magnetic surfaces-defined as the magnetic flux through a horizontal disc centred on the tokamak axis of symmetry, as in Ref. [33].

Kinetic equation
Dream implements the following model for the phase-space advection and diffusion: where A E describes acceleration in an electric field, A C collisional friction, A B the bremsstrahlung radiation reaction force, A S the synchrotron radiation reaction force and A T radial transport. Similarly, D C denotes collisional momentum-space diffusion and D T radial diffusion. The particle source is modelled as where C ava denotes the knock-on collision operator responsible for the runaway avalanche, and S p a particle source to model electron density variations, for example due to ion transport or ionizationrecombination.
In the following, we give explicit expressions for these terms in the bounce-averaged electron driftkinetic equation, as implemented in Dream. For the remainder of the paper, p = γv/c denotes the relativistic momentum normalized to m e c, with γ = (1 − v 2 /c 2 ) −1/2 denoting the Lorentz factor and v the speed. Other quantities are given in SI units, except temperatures which have the dimension of energy (the Boltzmann constant is k B = 1).

Electric field
The acceleration in the parallel electric field is described by the advection term with the non-vanishing bounce-averaged components where E = E · B/B, e is the elementary charge, and the average can be written

Collision operator
Collisions are modelled with a test-particle Fokker-Planck operator consisting of the advection and diffusion terms D C = m e T cold γν spp + m 2 e c 2 ν D 2 p 2 (I −pp), with the momentum unit vectorp = p/p. These have the non-vanishing components The collision frequencies ν s and ν D describe slowing down and pitch-angle scattering and, with the collision operator written in this form, support a Maxwell-Jüttner equilibrium distribution at temperature T cold independently of their value. These test-particle collision frequencies are given by the sum of the contributions from different sources.
For collisions with free electrons, the collision frequencies are taken from the relativistic Coulomb Fokker-Planck operator [34] ν ee where K 2 is the second-order modified Bessel function of the second kind and r 0 = e 2 /(4πε 0 m e c 2 ) the classical electron radius. Here, n cold and T cold refer to the density and temperature of the Maxwellian component of the electron distribution, around which the collision operator has been linearized. Ion collisions, assumed to be against infinitely massive targets, are accounted for using the model presented in Ref. [35], yielding the contributions where the sum is taken over all ion species (and charge states) in the plasma, Z i denotes their atomic number and Z 0i the charge number. The ad-hoc matching parameter k = 5, andā i is an effective-size parameter calculated and tabulated in Ref. [35], or evaluated using the recommended for- ln Λ 0 = 14.9 + ln T cold 1 keV − 0.5 ln n free 10 20 m −3 , with the free-electron density n free = i Z 0i n i , p T e = 2T cold /m e c 2 is the normalized thermal momentum and k = 5 is a matching parameter analogous to the one appearing in (17).

Radiation reaction force
Radiation losses due to bremsstrahlung are captured using a mean-force model [37] where effects of straggling [38] and screening are ignored. Synchrotron emission due to the electron gyromotion around the magnetic field line is accounted for using the advection term where p ⊥ = p − B · p/B denotes the perpendicular momentum. This synchrotron advection term has the components where τ S,min = τ S (B = B min ).

Avalanche source
Avalanche generation is modelled using the Rosenbluth-Putvinski [8] source term where n re is the density of runaway electrons and n tot = i Z i n i denotes the total density of electrons (free and bound). The inclusion of bound electrons in this source term provides an approximation for the energy spectrum of electrons created via ionization, which is valid when the energy of the created electron is much greater than the binding energy of the ion [39]. The conservative discretization of this source is detailed in section 4.2.

Radial transport
Radial transport is captured by prescribing the contravariant radial coefficients A r = (∂x/∂r) · A T and D rr = (∂x/∂r)(∂x/∂r) : D T . They can either be specified as arbitrary functions of (t, r, p, ξ 0 ), or via a Rechester-Rosenbluth model representing diffusion in fully stochastic magnetic field line regions [40] where (δB/B)(r, t) represents a normalized radial magnetic field fluctuation amplitude on the flux surface, R m is the major radius of the magnetic axis, q = q(r) is the safety factor calculated from the dynamically evolved plasma current density, and the step function H equals unity for passing particles and zero for trapped: with ξ T = 1 − B min /B max denoting the trappedpassing boundary. Since the Rechester-Rosenbluth model describes parallel transport along open field lines, trapped particles which bounce back and forth, will not undergo any net radial transport. The bounce average of (23) can concisely be formulated as {D rr } = πqR m c(δB/B) 2 v|ξ 0 | {ξ/ξ 0 }. The heat transport D W associated with this transport model is derived in Appendix B.1.

Particle source
In order to describe a time-dependent electron density due to ionization dynamics, an ad-hoc source term is added, of the form where electrons are created (or removed) with zero kinetic energy in the ion rest frame. The same source term is used for replacing radially transported fast electrons with cold electrons, required for quasi-neutrality to be maintained. The source amplitude S p (t, r) is treated as an additional unknown quantity which is solved for non-linearly in the Dream equation system, under the constraint of quasineutrality

Three electron populations
Dream allows the electron distribution to be evolved using the full kinetic equation (1), which is the most accurate but computationally expensive approach. The code also supports the solution of simplified equations for the electron dynamics which often provide adequately accurate results at significantly reduced computational cost. These approximations will be described in this section, and are based on the fact that electron dynamics are qualitatively different on three, typically well separated momentum scales; the corresponding characteristic momenta are 1. The thermal momentum p T e ∼ 0.01 at which the ohmic current, joule heating and many atomic-physics processes are dominant. 2. The runaway critical momentum p c ∼ E c /E ∼ 0.1 at which the runaway generation rate is determined. In this region, electrons that will be thermalized are separated from those runaway electrons that are accelerated towards ultra-relativistic energies. 3. The acceleration region p ∼ eE dt/m e c ∼ 100 where the energy spectrum (and pitch distribution) of the runaway tail is set. The dynamics in this region determine the synchrotron and bremsstrahlung radiation emitted by a runaway beam, as well as the current decay rate during the runaway plateau.
In order to resolve each of these regions efficiently, as well as to allow flexibility in approximating parts Figure 1: Illustration of the three kinetic regions defined in the fully kinetic mode. In this mode, the distribution function f hot describes both the cold (p < p hot ) and hot (p hot ≤ p < pre) electrons, while fre describes the runaway (p ≥ pre) electrons on a separate simulation grid.
of the electron dynamics (such as neglecting to resolve the energy spectrum of the runaway tail if only the generation rate is of interest), electrons in Dream are split into three separate populations: cold, hot and runaway (re) electrons which are, respectively, characterized by the densities n cold , n hot and n re . These three momentum regions, as well as the two associated electron distribution functions f hot and f re , are illustrated in figure 1.
In the general case, these are distinguished only by labels on different phase-space regions, and the dynamics are purely governed by (1) for all momenta.
In the subsections that follow, we describe approximations to the electron kinetics in each of these three regions.

Cold electrons
The Maxwellian component of the electrons, the cold population, is characterized by the density n cold , temperature T cold and the parallel (ohmic) current density j Ω . The density represents the density of all free electrons that are not labelled hot or runaway: where n free = i Z 0i n i is the total number density of free electrons, defined by the ion composition of the plasma, the sum taken over all ion charge states. The temperature is evolved according to a transport equation, to be described further in section 3.
The cold thermal component of the postdisruption plasma is not well-described by the bounce-averaged equation (1) due to its typically high collisionality, and in addition the test-particle collision operator employed in Dream is not adequate for resolving the ohmic current due to the lack of field-particle collisions. Therefore, the ohmic current is modelled according to where σ(n i , T cold ) denotes the parallel electric conductivity, and we implement the Sauter-Redl model [41] which accounts for neoclassical effects at arbitrary collisionality. The conductivity correction δj corr is an addition that is used to correct for transient currents when the thermal population is modelled kinetically (see section 2.2.2 below), and is given by where the step function H is defined in (24), p hot defines a threshold momentum (typically temperature dependent, p hot = 7p T e ) separating cold from hot electrons, and σ kineq represents the conductivity supported by the kinetic equation (1) in a fully ionized plasma; in such plasmas, δj corr = 0 in steady state. The δj corr term therefore provides a correction to the neoclassical conductivity in a rapidly time-varying plasma and accounts for the effect of partially ionized impurities. By matching the Ohmic conductivity σ kineq calculated with the test-particle collision operator for Z eff ∈ [0, 50] to the plasma conductivity calculated in Ref. [42], it was found that the formula provides an accurate approximation, where σ 0 is the Sauter-Redl formula in the collisionless limit.

Hot electrons
The hot electrons, described by a distribution function f hot , are modelled kinetically according to equation (1). They contribute density and parallel current density according to where the two threshold momenta p hot and p re differentiate hot electrons from cold and runaway electrons, respectively. There are two main modes of treating hot electrons in Dream: the fully kinetic mode where the full distribution functionincluding the thermal Maxwellian-is resolved kinetically on the grid, and superthermal where only the non-Maxwellian part of the distribution is followed: Fully kinetic. The full equation (1) is solved, in the form presented in section 2.1. Therefore, thermal electrons are resolved kinetically alongside hot electrons in the distribution function f hot as illustrated in figure 1, since the collision operator is density conserving and drives the solution towards a Maxwell-Jüttner distribution at temperature T cold . In this case, the cold density is given by corresponding to the density of electrons having momentum less than the hot-threshold p hot .
Superthermal. Only the superthermal electrons are resolved by taking the electron-electron collision frequencies as their T cold → 0 limits in which case the kinetic equation will no longer support a Maxwellian equilibrium distribution, but instead acquires a particle sink at p = 0 since the advective particle flux lim p→0 V pν s is finite. In this case, the hot-electron threshold is chosen to be p hot = 0, meaning that all electrons in f hot are considered hot, as illustrated in figure 2. By the cold-density equation (27), the electrons that are lost at p = 0 will be added to n cold . In this mode, the conductivity correction δj corr in equation (29) is dropped since ohmic current is not resolved kinetically. This approach was pioneered in Ref. [6] for thermal quench simulations.

Isotropic.
A further approximate form of the superthermal mode is supported by Dream, where the kinetic equation is analytically pitch-angle averaged based on an asymptotic expansion in which pitch-angle scattering is assumed to dominate the electron dynamics. The procedure mirrors the approximations employed in the calculation of the avalanche growth rate in ref. [8], but is generalized here to time-dependent situations which covers hot-tail formation. Details on the derivation of this reduced kinetic equation are provided in Appendix B.2, and the resulting expression is given by (B.15).

Runaway electrons
Runaways are described by the distribution function f re , which satisfies the same kinetic equation as the hot electrons f hot , but is defined on a separate grid spanning the interval p ∈ [p re , p max ], where the runaway threshold is typically chosen as p re ∼ 1, and the maximum resolved momentum is p max ∼ 100. Note that the boundary p re for the runaway grid is distinct from the usual critical momentum for runaway p c sometimes used to define a runaway electron in other codes. The runaway boundary p re used in Dream is constant in the simulation and must be appropriately chosen by the user so that both the hot and runaway electron distribution functions can be well approximated numerically. As long as the numerical grid is sufficiently resolved, the choice of the boundary p re has no effect on the simulation results.
The reason for separating hot and runaway electrons is that the electron distribution in the hot generation region (p ∼ p c ) is nearly isotropic, whereas in the runaway tail it can be extremely anisotropic. By utilizing a separate grid, resolution and discretization methods can be tailored to more efficiently resolve the runaway tail. Although the runaway electron distribution couples relatively weakly to the plasma evolution since they move with parallel velocities near the speed of light, the runaway distribution is essential when coupling to synthetic radiation diagnostic tools [43,44] or when considering kinetic instabilities driven by the runaways [45,46]. Particle conservation is ensured by enforcing a boundary condition on f re at the lower p boundary p=pre using a method described in section 4.5, where F hot/re denotes the particle flux through the shared p = p re boundary on the respective kinetic grids, and (∂ n re /∂t) fluid denotes the runaway generation rate due to sources that are not modelled kinetically on the hot grid. An example of such a source is Dreicer generation when using the superthermal mode, in which case the thermal Maxwellian is not resolved so that no electrons are available for Dreicer acceleration. The runaway density, corresponding to the number density of f re , therefore evolves according to If the runaway electron distribution function is solved for, the transport term in equation (36) is replaced with a source term which is the momentumspace integral of the transport term in the kinetic Table 1: Summary of the approximations used in the four main operating modes of Dream. In the fluid, isotropic and superthermal modes, the cold electrons are represented with a fluid model (here indicated with the density n cold and temperature T cold ). Notably, the representation of the hot population varies between the models, being represented either as part of the same distribution function f hot as the thermal electrons (fully kinetic), with a distribution function in the superthermal limit (superthermal), with an energy distribution in the superthermal limit (isotropic), or not at all (fluid).

Cold/thermal electrons Hot electrons
Fluid (1): In case the runaway electron distribution function is not explicitly evolved, the transport coefficients can instead be directly prescribed as functions of only radius, or be integrated over momentum space in the manner described in Ref. [47]. The parallel runaway current density is defined as The fluid runaway rate (∂ n re /∂t) fluid depends on the kinetic model employed, and support is implemented for modelling runaway generation due to Dreicer, hot-tail, Compton, tritium-decay and avalanche [13,48] depending on which equation terms are enabled. The fluid runaway models are described in Appendix C. If the hot grid is disabled altogether, the runaway distribution f re can still be evolved using fluid models for the generation via the boundary condition (35). Likewise, if the runaway grid is disabled, runaway generation is still captured via (36) where transport coefficients A re and D re can be imposed. In this case, the runaway current is instead determined by If neither hot or runaway electrons are resolved kinetically, a fully fluid-like system is obtained, corresponding approximately to the model contained in the Go code [13,49,50], generalized here to account for effects of arbitrary axisymmetric toroidal geometry.

Background plasma evolution
Dream utilizes a test-particle collision operator for the electron dynamics, which refers to the Maxwellian electron component of the plasma as well as the ion composition. In order to close the equation system, equations governing the evolution of these quantities as well as the electric field must be introduced. In this section we describe the evolution of the background plasma and how it couples to the electrons.
Most equations for the background plasma can be written in the one-dimensional transport equation form which has a structure very similar to the kinetic equation (1), but only evolves the quantity X in the radial coordinate r. In addition, it has the spatial Jacobian V instead of the full phase-space Jacobian V , and uses flux surface averages for its coefficients, defined by (5), instead of bounce averages. The similar structure is utilized when implementing the equations, as the discretized forms of the equations are near identical with only coefficients differing.

Ions
Ions are modelled by the densities n (j) i of each species i with atomic number Z i , and charge state j with charge number Z 0j , which are assumed to be uniformly distributed on flux surfaces. They are evolved according to where I and R denote ionization and recombination rate coefficients, respectively, which are extracted from the OpenADAS database [51]. The kinetic ionization rates are modelled by where the ionization cross-sections σ (j) ion,i are taken as in Ref. [52], which extended the validity of the Burgess-Chidichimo model [53] to relativistic energies, with the momentum integration limits taken as in section 2.2. Unlike the cited studies where model parameters were chosen from atomic data, Dream uses parameters which have been fitted to the OpenADAS ionization coefficients in the lowdensity (coronal) limit when the distribution is taken to be a Maxwellian for a range of temperatures, assuming that a single shell dominates the ionization. This method provides a smooth transition between kinetic and thermal (fluid) ionization rate coefficients, improving numerical stability of the solver. If the hot electrons f hot or the runaways f re are not resolved kinetically, the contribution to the ionization rate from such electrons is neglected.

Temperature
The background electron temperature T cold , which is assumed to be uniform on flux surfaces, is modelled via the evolution of the thermal energy W cold = 3 n cold T cold /2. The time evolution is governed by where the advection coefficient A W and diffusion coefficient D W are either prescribed functions of time and radius, or derived from the particle transport model, as described in Appendix B. The cold electron density has been pulled out of the radial derivative, consistent with the assumption that electrons cannot be transported independently of ions, in order to maintain quasineutrality. The collisional heat transfer to the cold population from hot and runaway electrons, as well as ions, is given by where we assume that hot and runaway electrons only deposit energy via elastic collisions with free cold electrons, thereby neglecting energy transfer by electron impact ionization. The sum i is taken over all ion species in the plasma, and it has been assumed that different charge states of the same ion species have the same temperature so that only the total density n i = j n (j) i and weighted i Z 2 0j of each species i appears. The ion temperature evolves according to the thermal-energy equation where W i = 3T i n i /2. The rate of energy loss L (j) i by inelastic atomic processes is modelled via (46) where L line is the radiated power by line radiation, L free by recombination radiation and bremsstrahlung, and the last terms represent the change in potential energy due to excitation and recombination, with the rate coefficients as in equation (41). The ionization threshold ∆W (j) i is retrieved from the NIST database [54], and the other rate coefficients are taken from OpenADAS.

Current, electric field and poloidal flux
The current is evolved via a mean-field equation for the poloidal flux [16], where the loop voltage is given by V loop = 2π E · B / B · ∇φ , and the generalized Ohm's law for the system is given by j tot = j Ω + j hot + j re , where the total parallel current density j tot satisfies Ampère's law Note that ψ here is distinct from the poloidal flux appearing in the definition of the magnetic field (7). The latter is taken as a static parameter which specifies the magnetic geometry and enters into all bounce and flux surface averages, while the ψ of equation (47) is a dynamically evolved quantity that sets the current profile evolution. The magnetic equilibrium (i.e. shape of the flux surfaces) is therefore kept constant during a simulation. The second term on the right-hand-side of equation (47), the hyperresistive term, acts to flatten the current profile due to magnetic field line breaking, while conserving the helicity content of the plasma [55]. The derivatives with respect to toroidal flux ψ t are taken at constant radius r. The helicity transport coefficient Λ(t, r), which is allowed to be arbitrarily prescribed, can be estimated e.g. by experimental observations of the temporal current spike and drop in the internal inductance in disruptions [56].
The boundary condition for the poloidal flux is modelled by introducing ψ edge = ψ(r = a) and ψ wall = ψ(r = b), where a is the (low-field side midplane) minor radius of the plasma and b represents the radius of the conducting wall. The flux at the edge couples to that at the wall via an approximate edge-wall mutual flux inductance M we = µ 0 R m ln b a with R m the major radius of the plasma, such that ψ edge = ψ wall − M we I p . The total toroidal plasma current is The poloidal flux ψ wall at the tokamak wall is evolved by where the loop voltage V (wall) loop at the tokamak wall can be modelled in two different ways, given as follows: Poloidal flux boundary condition 1: Prescribed. The wall loop voltage is provided as a prescribed time-dependent input, and the initial condition for the wall poloidal flux is chosen as ψ wall (t = 0) = 0.
Poloidal flux boundary condition 2: Self-consistent with circuit equation. The wall loop voltage is modelled by assuming an external inductance L ext = µ 0 R m ln(R m /b) and a wall resistivity R wall such that the characteristic wall time is where the wall current I wall is given in terms of the external inductance as

Numerics and implementation
The main part of the codebase constituting Dream is written in C++, with a rich interface written in Python. The C++ code is divided into two libraries that contain all routines necessary for setting up and running simulations, and one executable, which merely provides a thin commandline interface to the library routines. The difference between the two libraries is that one library contains lower-level routines for building finite volume stencils and interacting with the sparse linear algebra library PETSc [57,58], while the other library contains code for the physical models available, for initializing simulations, as well as for evolving the equation system in time.
In this section we will describe the details of the implementation, particularly the time, space and momentum discretizations used, but also the special boundary condition used to connect the hot and runaway electron kinetic grids. We end the section with a discussion about the performance of the code as well as some special methods that are used to improve performance of the code.

Finite volume discretization
To preserve the integral of the evolved quantities, we discretize equation (1) using a finite volume method [59]. This involves dividing the computational grid into N r ×N p ×N ξ cells with center values denoted using integer indices (i, j, k) (henceforth referred to as the distribution grid ), and points on the cell faces denoted with half integers, e.g.
(i − 1/2, j, k) (referred to as the flux grid ). In Dream, grids are defined by specifying the edge points z m 1/2 and z m Nm+1/2 of each flux grid along with the cell widths ∆z m i so that the flux and distribution grid points are distributed according to When evaluating second derivatives it is also useful to introduce the distribution grid spacing Next, the Fokker-Planck and transport equations (1) and (40), are averaged separately over each cell volume, resulting in a set of coupled equations for the average value of X in each cell. By using a central difference approximation for phase-space derivatives, and an Euler backward scheme for time derivatives, the discretized form of equation (1) becomes where l is the time step index and the notation i m ± N indicates that N should be added to/subtracted from the index of the m'th phase space coordinate, z m . The discretized form of equation (40) is almost identical. Since the advection and diffusion fluxes are evaluated on cell faces, the flux into any given cell on the computational grid will also be exactly the flux out of adjacent grid cells. This guarantees exact conservation of the integral of the quantity X within machine precision, in the absence of sources and edge losses.

Discretization of advection terms
The main difficulty of evaluating advection terms comes from the fact that the quantity X must be evaluated on a cell face rather than in the center of the cell, where it is actually computed. To do so we must interpolate in X based on its value in adjacent cells and many possible interpolation schemes could be used to this end. However, choosing the interpolation scheme with care can provide benefits such as improved stability of the numerical scheme as well as the preservation of monotonicity in X.
In Dream, we generally let with δ (im) k denoting a set of interpolation coefficients to be determined. The user can then select among a range of popular schemes, including linear schemes such as the simple centred (δ −1 = δ 0 = 1/2, for uniform grids), first or second-order upwind schemes, the third-order quadratic upwind scheme [60], or nonlinear flux limiter schemes such as SMART [61], MUSCL [62], OSPRE [63] or TCDF [64]. The nonlinear schemes are designed to preserve positivity of the solution. The flux limited schemes are upwind-biased, and for positive flow (A im−1/2 ≥ 0) can be expressed as δ −2 = −kφ(r), .
The flux limiter function φ(r) determines which scheme is used. For negative flows (A im−1/2 < 0), equation (58) is modified by mirroring all indices according to i m +x → i m −1−x and δ k → δ −1−k . The OSPRE and TCDF limiters prescribe continuously differentiable φ, making them robust choices with attractive convergence properties in the Newton solver, whereas the other piecewise linear schemes may sometimes fail to converge, but may provide higher accuracy.

Discretization of diffusion terms
The derivatives appearing in the diffusion terms in equation (1) are conveniently also discretized using a central difference approximation. This causes the diagonal (m = n) terms to take the form which also ensures that monotonicity is preserved for X. Off-diagonal diffusion terms require interpolation, and thus do not preserve monotonicity. In Dream, they are generally discretized as (60) An alternative to the above manner for discretizing off-diagonal diffusion terms was given in [65] where the terms were rewritten as advection terms and combined with a flux limiter scheme to ensure the preservation of monotonicity.

Discretization of avalanche source
The cell average of the bounce-averaged avalanche source (22) is evaluated according tô where the θ integral is taken over all angles for which ξ(ξ j−1/2 ) < ξ (p i ) < ξ(ξ j+1/2 ), with ξ(ξ 0 ) defined as in equation (3), ξ as defined in equation (22), and it has been assumed that n re /B is constant on flux surfaces (consistent with the assumption that runaways have ξ = 1 in the derivation of the local source function). Equation (61) is equivalent to equation (20) in Ref. [66]. We introduce a cutoff γ cut such that, if γ i−1/2 < γ cut < γ i+1/2 , we replace γ i−1/2 = γ cut , and we set the source to 0 for γ i+1/2 < γ cut . Defined this way, the numerical integral of the source term is with γ max = 1 + p 2 max , and p max denoting the maximum momentum resolved on the grid. This expression agrees with an exact integration of (22), independently of grid resolution.

Time evolution
As described in section 4.1, Dream uses an Euler backward (implicit) time discretization. All quantities on the right hand side of (56) should therefore be evaluated at time t = t l+1 , the same time for which a solution to the equation system is sought, thus requiring the nonlinear system to be solved iteratively. In Dream, a standard Newton's method is used. By moving all terms in (56) to the same side of the equality we obtain a system of coupled equations in homogeneous form which we denote with the operator F . By also denoting the vector of unknowns at time t l+1 by x (l+1) , we can write the system of equations compactly as In Newton's method we then linearize F around the true solution and obtain the iterative scheme where the index i indicates the Newton iteration and J −1 denotes the inverse of the Jacobian of F . The Jacobian is constructed based on analytical expressions for all equations, although some derivatives are approximated or even neglected altogether for simplicity and to ensure the sparsity of the Jacobian matrix. Notably, since the kinetic equation is relatively insensitive to the ion charge state distribution, a net performance gain is typically obtained by neglecting the ion Jacobian in the kinetic equation. Also, the advection interpolation coefficients (57) can be set to a linear upwind scheme in the Jacobian, even if a flux limiter method is used in the evaluation of the residual F ; the computational gain by reducing the number of non-zeros in the Jacobian can sometimes offset the additional iterations needed by the Newton solver to reach convergence.
An alternative approach to solving (63), suitable for linear systems and nonlinear systems which vary slowly in time, is obtained with a so-called linearly implicit or semi-implicit time discretization. With this scheme, the equation system (56) is linearized in time so that it can be written where M(x (l) ) is a matrix operator representing the differential operators and coefficients in (56), and S(x (l) ) is a vector representing the sources S in the same equation. In contrast to the Newton's method described above, M and S are here to be evaluated at the current time t l instead of at the time for which the solution x (l+1) is sought. As a result, a solution for equation (63) can be obtained using only a single matrix inversion, instead of a series of repeated inversions as required by the Newton method.
Both the linear and nonlinear solver options are available for all combinations of equations in Dream, allowing the user to easily switch between the two to compare performance and accuracy. To avoid duplicating implementations, the matrix M and the residual F are both built using the same routines, which receive a function pointer that is used to set a single element of M, or add to a single row of the residual F .

Convergence condition
The Newton iteration (64) should be terminated once the solution x (l+1) i+1 is sufficiently close to the true solution x (l+1) . This is typically done by evaluating an appropriate norm of the Newton step ∆x . This method is applied in Dream as well, but since the unknown vector x is a combination of all the unknowns of the equation system, we calculate the norm separately for each unknown quantity. To each unknown X n (r, p, ξ 0 ) we assign a pair of absolute and relative tolerances, abs n and rel n respectively, and demand that be satisfied for every unknown of the equation system X n in order for the solution to be accepted. The norm is taken to be the 2-norm of the solution vector corresponding to X n . By default, the relative tolerance is set to rel n = 10 −6 for all unknowns X n . The absolute tolerance, on the other hand, is disabled for most quantities (i.e. abs n = 0), with the exception of the runaway density n RE and current j RE , for which the absence of an absolute tolerance can sometimes cause the Newton solver to diverge 2 . We set the absolute tolerance for n RE to abs nRE = 10 −10 by default, and for j RE to abs jRE = ec abs nRE .

Treatment of fluid-kinetic boundary conditions
One of the novel features of Dream is the ability to seamlessly evolve different energy regions of the electron population using either fluid or kinetic models. The connection between the different energy regions is handled using specialized boundary conditions on the kinetic grids, and source terms on the fluid grids. Two particularly interesting use cases are when a runaway electron distribution function is included in an otherwise pure fluid simulation, as well as the case where two separate distribution functions are used to model the hot and runaway regions of the electron population.

Fluid runaway sources on kinetic grids
The use of a separate distribution function f RE for the runaway region in Dream makes it easy to resolve the high-energy runaway electron distribution function in simulations which otherwise only evolve fluid quantities. In such simulations, the electron population is split into a (free) thermal electron density n cold and a runaway density n RE . The production of runaway electrons is modelled using fluid generation rates, according to equation (36), with F hot = 0, which move particles from the thermal population to the runaway population.
When the runaway electron distribution function is included in such a simulation, the particles that are moved to the runaway density n RE must also be introduced to the runaway distribution function f RE . In reality, the details of how the particles should be introduced depend on the physics of the runaway production mechanism, but for the purpose of obtaining a lightweight fluid-kinetic model we introduce particles in the cell corresponding to (p, ξ 0 ) = (p min , sign(E)). This choice of source term is motivated by the fact that the newly introduced runaway electrons will generally be rapidly accelerated to even higher energies in the runaway electron distribution function, and thus move close to the ξ 0 = ±1 boundary anyway. The final form of the distribution function is generally not set by the kinetic details of the runaway source term, but rather by effects such as pitch-angle scattering, electric field acceleration, synchrotron emission etc. which are all included fully kinetically in this simulation mode. It is only the kinetic details of how particles cross into the runaway region which are not resolved in this mode.

Flow between hot and runaway distributions
The separation of the electron distribution function into a hot region and a runaway region can have significant computational advantages. The hot electron distribution function is often close to isotropic while varying rapidly with momentum, thus requiring only a few grid points in pitch ξ 0 , but potentially a large number of grid points in momentum p. The runaway electron distribution function, on the other hand, is typically aligned close to ξ 0 = sign(E) with a long, slowly varying tail in momentum. It thus often requires careful placement of pitch grid points, while only a handful of grid points in momentum may be necessary.
To handle the flow of particles between the hot and runaway electron distribution functions, a boundary condition connecting grid cells on both sides of the interface between them is introduced, as illustrated in figure 3. Since the two distribution functions may be defined on grids with widely different resolutions, care must be taken to ensure that the number of particles are conserved when electrons enter the runaway grid and vice versa. The boundary condition connecting the two grids is therefore based on the local density conservation equation where Φ (p),RE and Φ (p),hot denote the particle fluxes into the runaway grid and out of the hot grid, respectively, and the extent ∆ξ jJ by which the cells overlap in ξ 0 is By requiring the flux of particles to be locally conserved as in (67), one obtains for the advective and diffusive fluxes on the hot electron grid Φ (p),hot adv p hot Np+1/2 , ξ hot with the averaged runaway distribution function where J j denotes the smallest integer such that ξ RE Jj ≥ ξ hot j . The interpolation coefficients δ (1) j should be chosen as to minimize the risk of spurious oscillations on the grid boundary, and are therefore determined using an upwind scheme The fluxes on the runaway grid are determined by combining (67) and (69), yielding Note that the advection and diffusion coefficients F hot Np+1/2,j and D hot Np+1/2,j appear in the expressions for the fluxes on both the hot and runaway electron grids. To conserve particles it is necessary to use exactly the same coefficients for both grids.

A number of tests have been implemented for
Dream in order to verify the correctness of both the individual modules in the code, and the overall physics modelled. For the latter, benchmarks against results in the published literature have been performed and in this section we present the results of three such benchmarks. The first two tests verify that the plasma conductivity and Dreicer runaway rates are accurately computed by comparing Dream simulations to simulations with the 2D  Fokker-Planck solver Code [27,28], and primarily validate the Fokker-Planck collision operator used. Code only simulates homogeneous plasmas, but uses the same test-particle collision operator as Dream, albeit with a finite difference discretization in momentum and a Legendre polynomial decomposition in pitch. The third test is to reproduce the tokamak disruption simulations in [13], which were carried out with the 1D fluid code Go [67,49,68].

Conductivity
A typical method for validating Fokker-Planck collision operators is to solve the Spitzer problem arising in the presence of a parallel electric field E , where C is the collision operator, including collisions between electrons and electrons, as well as electrons and ions. In the weak electric field limit (E E c , with E c the critical electric field for runaway [3]) the current density j carried by the distribution function f will be proportional to E , with constant of proportionality σ, i.e. the conductivity of the plasma at the given temperature and effective charge. By solving equation (73) for f , and by extension the current density j, we can obtain the plasma conductivity from the relation σ = j/E .
In figure 4, the conductivity has been calculated with Dream (crosses) in the fully kinetic mode at a few different temperatures and plasma charges, and is compared to the conductivity as calculated with Code (solid lines). Both codes implement the fully relativistic test-particle collision operator of Ref. [69]. The values calculated with Dream are within less than 0.3% of those calculated with Code, indicating that the collision operator is correctly implemented.

Runaway rate
Another quantity of importance to the physics studied in Dream is the so-called Dreicer runaway electron generation rate obtained in a Z eff = 1 plasma when increasing E in the test above to E > E c [3]. In the fully kinetic mode, the number of runaway electrons n re is then defined as the number of particles with momentum p ≥ p re , where the runaway boundary is chosen as p re = 20 2T /mc 2 with T the electron temperature, and the runaway generation rate is taken as γ = ∂n re /∂t. Figure 5 shows the primary runaway rate as calculated with Dream (crosses) and Code (circles). The solid lines are calculated with the formula given in [3]. The x axis ranges from E = 2E c to E = 0.04E Dwhere E D denotes the Dreicer electric field [2] at which all electrons run away-corresponding to marginal and strong runaway electron generation respectively. The Dream and Code runaway rates match to within 3%.

GO ITER simulations
To validate the coupled physics of Dream we will now present a comparison between Dream and the simulations of ITER-like disruptions conducted in [13]. In ref. [13] the effect of injecting impurities in the plasma on the maximum runaway current in a standard ITER scenario was studied using  the 1D fluid code Go [67,49,68]. The Go code can be considered a predecessor of the fluid mode in Dream and uses similar models for the background plasma evolution. Using the fluid model, as described in section 2.2.3, with a cylindrical radial grid in Dream, the two codes should simulate approximately 3 the same physics. While benchmarking the two codes it was discovered that some of the simulations in [13] were slightly under-resolved with respect to time. While it does not lead to any qualitative differences, for this comparison we have re-run the Go simulations with improved time resolution. Figure 6 shows a comparison between the plasma currents and current densities obtained for cases 2, 3 and 4 in [13]. All three cases have the same temperature, main ion density and pre-disruption current density profiles, and only differ in the composition of the injected material. In all cases a mixture of neutral deuterium and neon is injected and is added instantaneously to the plasma, distributed uniformly across radii. The amount of material injected in the different cases is summarized in table 2.
In all three cases considered, Dream closely   Parameter Value Major radius R m 1.65 m Minor radius a 0.5 m Wall radius b 0.55 m Elongation at edge κ(a) 1.15 Toroidal magnetic field B 0 2.5 T Initial plasma current I p,0 800 kA reproduces both the total and runaway plasma currents in both the thermal and current quench phases of the disruptions, as well as the maximum runaway current density. The small deviations between the simulation results are explained primarily by the use of somewhat improved models for the critical electric field E eff c in Dream.

Comparison of electron models
To demonstrate some of the main features of Dream we will now examine two separate disruptions in a toroidal plasma with parameters representative of ASDEX Upgrade [70,71,72]. In section 6.1 we first describe the general parameters used for all simulations and briefly recall the differences between the electron models of Dream. Sections 6.2 and 6.3 then discuss the results and performance of each of the electron models.

Baseline simulation setup
All simulations of this section are conducted in an elongated ASDEX Upgrade-like plasma with magnetic field and vessel parameters as shown in table 3. Figure 7a shows the corresponding flux surfaces along with the plasma boundary (black) and conducting vessel wall structure (red). The flux surfaces are slightly elongated with a linearly varying elongation profile κ(r) = 1+0.15r/a. In figure 7b-d radial profiles of the initial electron density, temperature and current density are shown. The electron density is nearly uniform, close to n e,0 = 2.6 × 10 19 m −3 , while the electron temperature is peaked at T e,0 = 5.8 keV on the magnetic axis and decreases towards 60 eV near the edge. The plasma current density is j(r) = j 0 [1 − (r/a) 4 ] 3/2 , with j 0 = 1.52 MA/m 2 chosen to give the desired total initial plasma current I p,0 = 800 kA.
In the following sections we will insert a combination of neutral deuterium and neutral argon into the plasma outlined above. The material is assumed to be instantly distributed uniformly across the plasma. Using four different electron models, we follow the evolution of the plasma using as it cools down due to radiation losses as well as a prescribed diffusive heat transport, according to (43), with D W = 4000 m 2 /s in section 6.2 and D W = 1000 m 2 /s in section 6.3. The four models used are the fully kinetic, the superthermal and the isotropic models described in section 2.2.2, as well as a fluid model similar to the one used in the Go code [13,49,50], which we briefly commented on in section 2.2.3. The main differences between these models can be briefly summarized as follows: in the fully kinetic model, both cold and hot electrons are modelled kinetically while runaway electrons are modelled as a fluid; in the superthermal model, only hot electrons are modelled kinetically, while cold and runaway electrons are modelled as fluids; the isotropic model makes the same assumptions as the superthermal model, but uses an angle-averaged kinetic equation and evolves only the energy distribution of the hot electrons; in the fluid model, only thermal bulk and runaway electron populations are followed, both as fluids.
Another important difference separating the fully kinetic/fluid and superthermal/isotropic models is the treatment of the cold electron temperature, T cold , which is used both in the test-particle collision operator and in the ion rate equations (41). In all models the temperature is evolved according to equation (43), and the electron distribution is initialized at equilibrium with the temperature of the hot initial plasma. In the fully kinetic and fluid models-which do not distinguish between cold and warm electrons-the temperature T cold starts at the initial warm plasma temperature, and rapidly falls as cold impurities are inserted in the plasma. In contrast, in the superthermal and isotropic models, T cold starts at almost zero, corresponding to the temperature of the injected impurities. The injected electrons will promptly form a Maxwellian at a significantly lower temperature than that of the pre-disruption plasma, and the hot electrons will predominantly slow down in free-free collisions with this cold Maxwellian. Because of these differences, the fully kinetic and fluid models can be expected to agree relatively well, while the superthermal and isotropic models could be expected to differ from the former two in some cases.
The cooling-down process occurring as impurities are injected into the hot plasma leads to an  inherently non-linear evolution for the distribution function as it transitions from containing two Maxwellian electron populations-the injected electrons at a very low temperature and the initial bulk electrons at a warm, but gradually cooling, temperature-into a state dominated by a single Maxwellian electron population. The superthermal and isotropic models provide a numerically efficient way of capturing some of this dynamic without resorting to a fully non-linear collision operator.

Full-conversion scenario
In this scenario we initiate the plasma as described in section 6.1 and insert neutral deuterium and argon with radially uniform densities n D = n Ar = 2.6 × 10 19 m −3 at t = 0. As shown in figure 8, the resulting temperature and current dynamics are very fast, with much of the thermal quench completing in about one hundred microseconds. Due to the rapid cooling, a significant fraction of electrons remain hot at the onset of the current quench, leading to a seamless conversion of ohmic current into superthermal and then runaway current via the hot-tail mechanism. By the end of the current quench, almost all of the original current has been converted into runaway current in all three models.
Some differences are observed in the evolution of the temperature and total runaway current in the fully kinetic/fluid and superthermal/isotropic models, although the total plasma current reached is almost exactly the same in all cases. Note that the temperature shown for the superthermal and isotropic models is that of the cold injected electrons, while the temperature shown for the fully kinetic and fluid models is that of the (initially warm) bulk electrons. The main reason for the differences in temperature evolution is the electron-ion heat exchange Q ij in equation (44), which explicitly depends on the relative temperature difference between electrons and ions. Since the initial value of T cold differs in the four models, so does the heat transferred from the ions to the electrons, and hence also the detailed evolution of T cold . The fact that the runaway currents still agree well is due to the thermal quench being rapid in all models, which leads to the formation of a significant hot electron population that is eventually accelerated and runs away. Figure 9 illustrates the evolution of the current carried by the distribution functions in the superthermal and fully kinetic models. In the former, the ohmic current is obtained from Ohm's law (as in equation (28)) and the distribution function only carries superthermal current, while in the latter the distribution function also contains the thermal bulk electrons, and thus also carries ohmic current (seen as a sharp peak near p = 0 in figure 9b). During the current quench, the remaining superthermal electrons-which carry a significant fraction of the total current-are accelerated to even higher momenta and run away. The absence of the thermal bulk in the superthermal model permits the momentum grid resolution to be nearly halved, making the superthermal model computationally efficient while accurate in capturing the hot-tail formation.

Slow disruption scenario
In this scenario we initiate the plasma as described in section 6.1 and insert neutral deuterium and argon with radially uniform densities n D = 5.2 × 10 20 m −3 and n Ar = 5.2 × 10 18 m −3 respectively at t = 0. The resulting disruption occurs over a relatively long time, with the current quench lasting for up to 4 ms, depending on the model used, as shown in figure 10b. As in section 6.2, the temperature shown in figure 10a for the superthermal and isotropic models is that of the cold injected electrons, while the temperature shown for the fully kinetic and fluid models is that of the (initially warm)  bulk electrons. In contrast to the scenario of section 6.2, this scenario reveals significant differences between the different models used. The fluid and fully kinetic models mostly agree with each other, as do the isotropic and superthermal models, but when comparing the superthermal and fully kinetic models-the two most advanced models-the final runaway currents are found to deviate by a factor of two. The deviations between the fluid/fully kinetic and isotropic/superthermal models in figure 10 are consequences of the self-consistent plasma evolution, although the origin of the different evolutions can be traced to the disparate definitions of the temperature used. Several terms and coefficients depend explicitly on the temperature T cold , including the ion rate coefficients in equation (41), the collisional energy transfer term (44), and the radial heat diffusion term, and as such, differing dynamics are to be expected in the brief initial phase of the TQ when T cold differs significantly between the two groups of models. However, it is only if one or more of these temperature-dependent terms are dominant during the early TQ phase that the final runaway current should be significantly impacted. In the scenario of figure 10 it turns out that the radial heat diffusion term plays an important role in the fluid and fully kinetic models early during the TQ, causing the T cold profile to be flattened. This in turn alters the behaviour of the electric field, which typically grows rapidly in response to the decreased conductivity when the temperature drops.
With the fluid and fully kinetic models, the relatively strong radial heat diffusion causes the temperature to decrease rapidly in the centre of the plasma, but also to be slightly raised at outer radii. As a result, the remaining thermal energy is radiated away more slowly, allowing the ohmic current to be maintained for a longer time, and thus delaying the increase of the electric field. Figure 11 shows the electric field evolution in the superthermal and fully kinetic models during this early phase of the disruption. The strong electric fields during the early phase of the disruption leads to a significant conversion of hot electrons to runaways. In the fully kinetic case, figure 11b, the slower electric field evolution does not allow for as many hot electrons to be immediately converted into runaways, but partially compensates for this later on during the disruption by driving more production of runaways through the avalanche mechanism. The increased avalanche generation in the fully kinetic model is however not sufficient to fully compensate for the early hot-tail generation in the superthermal model.
The deviations between the fluid and fully kinetic models, as well as the isotropic and superthermal models, stem almost entirely from the differences in how the hot-tail generation is modelled. Both the fluid and isotropic models utilise approximations to the Fokker-Planck treatment of the hottail mechanism used in the superthermal and fully kinetic modes. As a result, the number of runaway electrons generated via the hot-tail mechanism is slightly over-and underestimated, respectively, in the fluid and isotropic modes.
Finally, a comment on the physics fidelity of the four considered models is due. The fluid and isotropic models are direct approximations of the more advanced models and are less reliable, although the results here suggest that their results for the temperature and current evolution are reasonably close to the more advanced kinetic models. As for the superthermal and fully kinetic models, it is difficult to clearly state that one is more reliable than the other. The fully kinetic model uses a linearized collision operator and in the early phase of the thermal quench, the process may be inherently non-linear. When a large amount of impurities are injected, the electron distribution will briefly be constituted by two Maxwellians at different temper- atures, which the fully kinetic model is not equipped to handle. The superthermal model, on the other hand, is derived with this exact situation in mind and therefore provides a better approximation of the processes. The superthermal model is however still an approximation to the full disruption physics, and its assumptions, e.g. that a large number of impurities are present, are not necessarily always well satisfied. To verify the hot-tail models considered here, a relativistic non-linear collision operator, as used in Refs. [29,73], coupled to a self-consistently evolving background plasma, is therefore needed. This will be considered in future work.

Summary
The main purpose of the Dream code is to model the self-consistent plasma evolution and runaway electron generation during a tokamak disruption. The output is the evolution of the temperature, densities of the different particle species and poloidal flux (which sets the evolution of the current density) as well as the electron distribution function. The temperature evolution includes ohmic heating, radiated power using atomic rate coefficients, collisional energy transfer from hot electrons and ions, as well as dilution cooling. The poloidal flux evolution includes the option to model rapid current flattening associated with fast magnetic reconnection events, via a helicity-conserving hyperresistivity term. The fluid quantities are solved on a one-dimensional flux-surface averaged grid, and the kinetic equation for the electron distribution is solved in a three-dimensional (1D-2P) bounceaveraged formulation.
The ability to treat electrons at various degrees of sophistication is one of the more novel contributions of Dream. The physics of tokamak disruptions typically involves multiple temporal and spatial scales, and so far often required comprehensive and computationally expensive simulations involving both fluid and full kinetic physics. By separating the electrons into cold, hot and runaway populations, and evolving the cold and runaway electrons using fluid models, Dream avoids resolving the usually uninteresting-but computationally intensive-kinetic bulk and runaway tail dynamics. Furthermore, the possibility to evolve hot electrons using a pitch angle-averaged kinetic equation allows simulation times to be reduced almost to the level of pure fluid models while the electron hot-tail generation is still accurately captured.
Together with the comprehensive physics model, which reaches beyond previous efforts in kinetic disruption modelling, the code is equipped with several attractive numerical features including fully implicit time stepping of the full system, as well as a flux conservative and positivity preserving discretization, which contributes to the robustness of the tool. Due to its flexibility and numerical efficiency, Dream is suitable for extensive investigations of disruption and runaway physics.
Dream has been verified against both kinetic and fluid codes, and has been found to reproduce their results in the appropriate limits. As an application, two disruption scenarios were investigated in an ASDEX Upgrade like tokamak, where the disruption is triggered by the injection of a combination of neutral deuterium and argon atoms. The full hierarchy of electron models was compared: fully kinetic, superthermal, isotropic and fluid models, and reasonable agreement was found. The difference in simulation time between the fluid and kinetic simulations is more than two orders of magnitude.
evaluated with a Gauss-Legendre quadrature rule where the four independent angle-dependent quantities {η n } 4 n=1 as well as the Jacobian J = 1/|∇ϕ · (∇r × ∇θ)| are precomputed at the nodes θ i , where N = 10 is found to typically be sufficient for relative errors < 10 −4 . For bounce integrals over trapped orbits, the integration limits θ b1 (ξ 0 ) and θ b2 (ξ 0 ) depend on pitch and the metric has an integrable singularity (∼ 1/ θ 2 b1,b2 − θ 2 ) at the boundary, and a Chebyshev-Gauss quadrature is employed instead. In this case, the precomputed values must be separately evaluated on individual poloidal angle meshes for each trapped orbit, in contrast to passing orbits for which the same poloidal angle mesh can be reused.
The bounce-orbit metric V contains a logarithmic singularity on the trapped-passing boundary ξ 0 = ±ξ T . In order to resolve these singular points, instead of using the midpoint rule to estimate the cell average in the finite-volume methods, for cells containing a singular point we instead carry out the pitch average as a double integral using the adaptive QAWS routine of QUADPACK, which is designed for integrals with endpoint singularities. For brevity, we have not written out other arguments on which the integrand may depend. Dream supports the use of an analytic up-down symmetric geometry where the flux surfaces are parametrized according to Here, the Shafranov shift ∆, elongation κ and triangularity δ parametrize the shape of the flux surfaces, and G and the reference poloidal flux gradient ψ ref (that is left independent of the poloidal flux ψ, evolved dynamically in the equation system) determine the strength of the toroidal and poloidal components of the magnetic field, respectively. In this geometry, the Jacobian J and scale factor |∇r| 2 are given by J = κrR cos(δ sin θ) + ∆ cos θ + sin θ sin[θ + δ sin θ]

Appendix B. Kinetic equation
Appendix B.1. Heat transport due to fast-electron transport The heat transport associated with the electron particle transport is obtained by integrating the diffusive transport term of the kinetic equation with the diffusion coefficient (23) over a Maxwell-Jüttner distribution function, where we do not keep the contribution from ∂n cold /∂r since electron density cannot be transported independently of the ions; the electron density profile is set by quasineutrality. As such, we assume that the heat transport acts to flatten the temperature profile. Therefore, the heat diffusion coefficient is given by the energy moment (m e c 2 (γ − 1)) For the Rechester-Rosenbluth model (23) the pitch integral is given by with the last line corresponding to the nonrelativistic limit Θ 1, retaining the leading-order relativistic correction.

Appendix B.2. Reduced kinetic equation
For the isotropic electron mode, described in section 2.2.2, a reduced form of the usual kinetic equation detailed in section 2.1 is used, analogous to Ref. [8]. To derive the reduced equation, we first introduce the ordering parameter δ and assume ν D ∼ δ 0 and E ∼ δ 1 , with all other terms being of order δ 2 . Writing f = f 0 + δf 1 + O(δ 2 ) and grouping the kinetic equation by order in δ, we obtain a set of new equations governing the evolution of the distribution function: (B.6c) where in the last equation,Â p denotes the momentum advection that is not due to electric-field acceleration, and F ξ0 denotes the net pitch flux, which will not affect the final result. Solving the δ 0 equation (B.6a) first yields the leading order solution i.e. f 0 is isotropic. Substituting this into the next order equation (B.6b) and integrating the result where the step function H was defined in (24). A solution for f 1 can then be obtained by isolating ∂f 1 /∂ξ 0 and integrating over ξ 0 once more, giving .

(B.10)
With f 0 and f 1 determined, we may now use the δ 2 equation (B.6c) to obtain the final reduced kinetic equation. We first multiply both sides of equation (B.6c) by V and integrate over all ξ 0 . Because of this, the pitch angle scattering term, as well as the ξ 0 component of the electric field term, vanish due to the factors of 1 − ξ 2 0 which are zero at ξ 0 = ±1. The slowing down and transient terms remain mostly unaffected by the integration, only picking up a factor of (B.11) while the electric field term becomes where the expression (B.10) for f 1 was substituted.
After interchanging the order of integration, the ξ 0 integral ranges over the passing region only due to the step function, while the ξ 0 integral is recognized as half the function h(ξ 0 ). Since we will then only evaluate h(ξ 0 ) in the passing region, we may replace it by its actual value there, h(ξ 0 ∈ passing) = 1 − ξ 2 0 . Furthermore, since the integrand is even in ξ 0 in the passing region, we may replace the integral with two times the integral ranging from ξ T to 1: (B.13) The electric field term, which was purely advective in the original kinetic equation, has now become purely diffusive after the averaging procedure.
For the radial transport terms the ξ 0 integral will only apply to the advection and diffusion coefficients, and so we can make use of the relation dξ X ≡ 4πp 2 X ξ , (B.14) where in the next-to-last step we used that B min ξ dξ = Bξ 0 dξ 0 . In what remains, we will use the notation X ξ to denote the combined flux surface and pitch average of a quantity X.
After the steps above, we obtain the final reduced kinetic equation by dividing all terms by V dξ 0 = 4πp 2 V , yielding (B.15) Here we have introduced D E as the momentum diffusion coefficient representing the electric-field acceleration, and the effective passing fraction, de-noted f p , is defined as Current density. The current density corresponding to the distribution evolved by the reduced kinetic equation can be calculated from equation (32). Since f 0 is isotropic, f 1 is the lowest-order component of f to carry any current, and the current density therefore becomes .
(B.17) Just as for the electric field term, we interchange the order of integration in ξ 0 and ξ 0 and recognize that the ξ 0 integral yields a factor of (1/2)h(ξ 0 ). The current density is hence given by An undesirable property of (B.18) is that it allows for the current density to grow larger than ecn e , which is the current density expected when all electrons travel at the speed of light along magnetic field lines. Exceeding this value may destabilize the solver, and so we adjust for this behaviour by smoothly matching the limit of all electrons travelling parallel to magnetic field lines, in which case the current density is given by j iso,2 B = 4πe sgn( E · B ) dp vp 2 f, (B. 19) where sgn(x) denotes the sign function. The actual current density used in the simulation is then taken as the matched formula j = j iso j iso,2 corresponding approximately to the smallest of the two approximations.
Appendix C. Runaway fluid formulae in general tokamak geometry A widely used formula for the avalanche growth rate in a pure plasma was derived by Rosenbluth & Putvinski [8], accounting for geometric effects in a large-aspect-ratio tokamak. This formula has later been generalized to partially ionized plasmas [74,75], which has also been applied to runaway generation due to Compton scattering and tritium beta decay [48,13]. Here, we present generalized fluid formulae for the runaway generation rate that extends the validity of previous work to axisymmetric tokamak geometry with shaped surfaces of arbitrary aspect ratio. approximated by a step function at the critical momentum p defined by (C.10) The validity of the runaway rates can be extended to the near-threshold regime E ∼ E c and to weak pitch scattering ν D ∼ 0 by defining a matched formula for the critical runaway momentum p c according to [75] p c = ν s (p )ν D (p ) + 4ν s (p ) 2 f p e 2 ( E · B / B 2 − E eff c ) 2 1/4 , (C.11) whereν s = p 3 ν s /γ 2 andν D = p 3 ν D /γ denote normalized collision frequencies that depend only weakly on momentum. In this expression E eff c denotes the effective critical field, described in Appendix C.2, and p c = ∞ when the electric field is smaller than E eff c . Compared to ref. [75], the numerator includes a factor ofν 2 s which increases the accuracy of the formula for weakly ionized low-Z plasmas. In Dream, the runaway rate due to source functions are evaluated using equation (C.7) with the runaway probability h(p) = H(p − p c ), (C. 12) where H denotes the Heaviside step function. The source function S due to Compton and tritium decay are modelled as in ref. [13].
where the distribution-weighted bounce averaged momentum flux-or net acceleration-U (p) is The effective critical electric field E eff c is then defined as the minimum value of the electric field for which there exists a real solution to U (p) = 0, that is

Implementation
The calculation of U (p) typically requires repeated evaluation of three nested integrals, two of which stand inside an exponential function, and is hence not entirely straightforward to implement efficiently. To speed up evaluation we use splines to evaluate the function g(ξ 0 ) in (C.15) as well as bounce averages of the advection coefficient A p . For the function g(ξ 0 ), we construct splines representing the integrand ξ 0 / ξ in (C.16) by evaluating the integrand on a uniform ξ 0 reference grid, which subsequently allows the function H(ξ 1 , ξ 2 ) to be efficiently evaluated using routines for exact integration of splines.
The advection coefficient A p can be factorised into A p = i a p i (p)Â p i (ξ 0 , θ) with the sum i taken over equation terms contributing to the force balance, and where the prefactor depends only on momentum and the remainder only on pitch and poloidal angle. The bounce average of the pitchdependent part of the advection coefficients, {Â i p }, are then spline interpolated onto a uniform pitch grid in the interval ξ 0 ∈ [0, 1], since all advection operators considered are either symmetric or antisymmetric in ξ 0 . The distribution-weighted bounce average of the coefficientsÂ p are then spline interpolated to the uniformly sampled variable X = A 2 /(1 + A) 2 ∈ [0, 1] (in which the functions are smoothly varying all the way up to the limit A = ∞, corresponding to all runaways having ξ = 1), where A is the inverse pitch distribution width parameter given in (C.15), allowing rapid evaluation of the acceleration function U . The root of (C. 19) is then solved for as a nested optimization problem with two layers. In the outer layer, a solution is sought to the one-dimensional root finding problem U e ( E · B ) = 0, (C. 20) where U e is the maximum of U (p) with respect to p at a given electric field E · B , i.e. the strongest acceleration experienced by any particle. The problem is solved using an unbounded secant method, assuming for the initial guess that E eff c /E tot c is constant in time, where E tot c denotes the classical critical electric field given in Ref. [3], evaluated with n e being the density of both free and bound electrons.
In the inner layer, the strongest acceleration at any momentum U e = min is determined. This problem is solved using Brent's method [77] from the GNU Scientific Library [78].
To ensure fast and robust convergence of the method, the algorithm is applied to the interval p opt (1 ± 0.02), where p opt is the minimum from the previous solve. If the interval does not contain the minimum, it is expanded in steps of 20% until a minimum is enclosed. Expansion of the interval is typically needed less than once in a thousand solves.
runaways, resulting in increased accuracy at high