I Introduction

# Numerical calculation of ion runaway distributions

## Abstract

Ions accelerated by electric fields (so-called runaway ions) in plasmas may explain observations in solar flares and fusion experiments, however limitations of previous analytic work have prevented definite conclusions. In this work we describe a numerical solver of the 2D non-relativistic linearized Fokker-Planck equation for ions. It solves the initial value problem in velocity space with a spectral-Eulerian discretization scheme, allowing arbitrary plasma composition and time-varying electric fields and background plasma parameters. The numerical ion distribution function is then used to consider the conditions for runaway ion acceleration in solar flares and tokamak plasmas. Typical time scales and electric fields required for ion acceleration are determined for various plasma compositions, ion species and temperatures, and the potential for excitation of toroidal Alfvén eigenmodes during tokamak disruptions is considered.

Numerical calculation of ion runaway distributions

O. Embréus, S. Newton, A. Stahl, E. Hirvijoki, T. Fülöp

Department of Applied Physics, Chalmers University of Technology, SE-412 96 Göteborg, Sweden

CCFE Fusion Association, Culham Science Centre, Abingdon, Oxon, OX14 3DB, United Kingdom

## I Introduction

The phenomenon of particle runaway in a plasma is well known, occurring in both space and laboratory plasmas (1); (2). It arises because the friction force experienced by a charged particle decreases with particle energy, so that the presence of a sufficiently strong induced electric field can allow the particle to be accelerated - or runaway - to high energy.

Electron runaway (3) has been extensively studied in the context of magnetic confinement fusion in tokamaks, where it can lead to the formation of localized high-energy beams which must be carefully controlled (4). The standard analytic method used to determine the initial evolution of the electron distribution function in a fully ionized plasma was introduced by Kruskal and Bernstein (5) and later generalized by Connor and Hastie (6). It takes the form of an asymptotic expansion of the electron kinetic equation in the electric field strength. Once a fast electron population – known as the primary distribution – is established, it can rapidly produce further fast electrons through large-angle, or knock-on, collisions (7). This avalanche process leads to the so-called secondary runaway-electron generation, which is often dominant.

Ion runaway has long been of interest in the astrophysical community, where it is thought to contribute to the observed abundance of high energy ions in solar flares (1). It has also been used to study the behavior of lightning channels (8) and was observed in laboratory plasmas, e.g. in the Mega Ampere Spherical Tokamak (MAST) (9) and in the Madison Symmetric Torus (10). The detailed mechanism of ion runaway differs from that of electron runaway. Friction with the electrons, which are also drifting in the electric field, acts to cancel a portion of the accelerating field. In the ion rest frame, the electrons have net motion anti-parallel to the electric field, and a test ion will experience two main forces beyond friction against the background of charged particles: acceleration in the electric field, and friction due to the electron drift. These forces scale differently with ion charge, and the dominant force is either electron friction – with the consequence that the ions are dragged along with the electrons – or acceleration by the electric field. In a pure plasma the cancellation of electric field acceleration and electron friction is complete, but the presence of impurities, neutrals or effects such as particle trapping in a non-uniform confining magnetic field allow a finite effective field to remain (11); (12).

Furth and Rutherford (12) generalized earlier treatments by adopting a similar expansion procedure to that used to study electron runaway. They determined the steady state ion distribution function in conditions typical of operational fusion plasmas. Helander et al(9) then considered the initial value problem resulting from the onset of an accelerating electric field, produced for example by a plasma instability. An analytic solution for the initial time evolution of the accelerated ion distribution function was developed, but it was noted that its application was limited due to the low electric fields required for it to be valid. Both of these ion runaway studies considered only the presence of a trace impurity population, consistent with typical operating conditions in fusion plasmas.

Plasmas with significant impurity content are also common, however. Astrophysical plasmas often consist of a mixture of dominant species, as well as containing trace elements (1); (13). Disruptions (2), which are a common cause of electron runaway in tokamaks, are also typically associated with an increase in impurity content, either due to deliberate gas injection for mitigation purposes or due to plasma-wall interaction. Therefore, in Ref. (14), the initial value formulation of the problem posed in Ref. (9) was extended to account for arbitrary plasma composition. The potential for ions accelerated during a disruption to excite low frequency plasma instabilities was considered analytically. The results were inconclusive since the asymptotic expansion used to develop the analytic solution was not strictly valid for disruption-type parameters. The limitations of the analytic solutions available in previous work motivate the development of a numerical tool to allow detailed study of the time evolution of an ion runaway distribution.

Here we describe the formulation and implementation of an efficient finite-difference–spectral-method tool, CODION (COllisional Distribution of IONs), that solves the two-dimensional momentum space ion kinetic equation in a homogeneous plasma. CODION solves the ion Fokker-Planck equation as an initial value problem and allows for time-variation of the electric field and bulk distribution parameters (temperature, density, charge and mass) of each plasma species independently. Due to its speed it is highly suitable for coupling within more expensive calculations, e.g. studies of instabilities driven by the fast ions or comprehensive modeling of ion acceleration with self-consistent coupling to solvers of Maxwell’s equations. Using CODION we obtain illustrative two-dimensional ion velocity space distributions, which demonstrate the typical behaviour of runaway ions in a variety of physical scenarios. We show that during typical tokamak disruptions ions are unlikely to be accelerated to velocities high enough for resonant interaction with toroidal Alfvén eigenmodes (TAE). Therefore the experimentally observed TAE activity cannot be explained by the ion runaway mechanism.

The rest of the paper is organized as follows. In Section II we describe the ion Fokker-Planck equation, and in Section III we outline its numerical implementation in CODION. In Section IV we explore the numerical solution, including the effect of various approximations to the collision operator. In Section V the model is applied to a variety of physical scenarios, illustrating typical acceleration time scales in laboratory and space plasmas. Finally we close with concluding remarks in Section VI.

## Ii Runaway Ion Distribution

We consider the problem of ion acceleration by induced electric fields with a component parallel to the background magnetic field in a plasma. We restrict ourselves to straight field line geometry, and assume a homogeneous background plasma. The time evolution of the ion distribution is then given by the Fokker-Planck equation, and particle acceleration will be opposed by various friction forces. The effect of friction with neutral particles can be significant in various physical situations, for example giving rise to charge exchange losses, which was studied in the context of lightning discharges in Ref. (8). Here we will focus on fully ionized quasi-neutral plasmas, in which case the friction is the result of inter-species Coulomb collisions only.

We are interested in the initial value problem where an electric field appears in what was previously an equilibrium state. Therefore we assume that the initial particle distribution functions are stationary Maxwellians – possibly at different temperatures – and consider their distortion from this state by the electric field. We linearize the collision operator about this background Maxwellian, and neglect the non-linear contribution to the evolution. This restricts the study to situations where only a small fraction of the ion population is accelerated, or to the initial stages of ion runaway. A reason for why this is sufficient in the present context is that, once a high energy population forms, the runaway ions have the potential to excite instabilities (14), which will have a strong impact on the further evolution of the distribution. Note that the non-linear terms of the kinetic equation are sometimes required to account for the transfer of energy from the electric field into heating the distribution. These aspects of the longer term evolution of the distribution are beyond the scope of the work presented here.

The linearized collision operator consists of the test and field-particle terms. Momentum and energy conservation in ion self-collisions is retained by approximating the field-particle piece with restoring terms, as described in Ref. (15), without the need to include the full field-particle operator which is numerically more expensive. Collisions with the other ion species could be treated similarly, however this would require the simultaneous evolution of the distribution functions of multiple species. Therefore, we consider only the evolution of the ion species with the highest runaway rate, so that the other ion species remain approximately stationary and only the test-particle piece of the unlike ion collision operator needs to be retained. While it is difficult to verify a priori that this condition is satisfied, it can be determined by numerical simulation of each ion species individually, assuming the others to remain stationary. Because of the sensitivity of acceleration rate to ion charge and mass (as demonstrated in Section V), the condition can typically be well satisfied as the acceleration rate of different species is often separated by several orders of magnitude.

The velocity-dependent friction on a test particle resulting from collisions with a Maxwellian background species has a peak near the thermal velocity of the background due to the form of the Coulomb interaction between charged particles. In the case of electrons, the friction force will be a monotonically decreasing function for velocities above the electron thermal velocity, allowing an electron to run away to large energy (where relativistic (6) and synchrotron effects (16) become dominant). We focus on situations where the ion, , and electron, , thermal velocities satisfy , meaning that their temperatures are sufficiently similar that . Then, if an ion is accelerated away from the bulk, friction against the electron population will increase with velocity. This has the consequence that the ion acceleration will be naturally balanced by the electron friction for some , if the electric field is below a threshold value similar to the Dreicer field (17) for electrons.

Since an electron reacts to the electric field on a time-scale times shorter than the ions, we assume the electron population to be in a quasi-steady state on all time-scales of interest for ion acceleration. Parallel force balance for the electron distribution then requires that the total electron-ion friction cancels the acceleration by the electric field, , where the sum is taken over all ion species in the plasma. Due to the small mass ratio, the electron-ion interaction is dominated by pitch-angle scattering, so that , and we can solve for the friction between electrons and the ion species of interest, , where the effective charge is .

The electron distribution can be written as , with a Maxwellian distribution drifting with the bulk ion velocity , and a correction varying over velocities of order , accounting for the electron drift behavior in the electric field. Then the linearized ion-electron collision operator, neglecting terms quadratic in , can be simplified (18); (1), noting that momentum conservation in binary collisions requires that ,

 Cie{fi,fe}=Reimini⋅∂fi∂v+Cie{fi,fMe(v−Vi)}. (1)

The first term, which describes the friction arising from the drifting electron distribution and was calculated above, readily combines with that describing acceleration by an electric field, giving rise to the effective electric field . Thus, as noted in the Introduction, in a pure plasma where , net ion acceleration will not occur. Light ions with can be accelerated in the direction of the electric field. Heavy impurities with will be accelerated in the opposite direction, that is in the direction of electron runaway (the latter case was studied by Gurevich (11)). The second term in Eq. (1) describes the slowing down of the fast ions on the electrons, as well as the slow energy exchange between the bulk species. Note that the ion flow velocity correction is time dependent and formally small in the runaway density. This term will become more significant in the ion-electron momentum exchange as the runaway distribution builds up.

Finally then, the kinetic equation we consider for the evolution of the ion distribution function in the presence of an accelerating electric field and arbitrary plasma composition is the following:

 ∂fi∂t+ZiemiE∗(ξ∂∂v+1−ξ2v∂∂ξ)fi=∑sCis{fi}, (2)

where , and the effect of collisions with the background Maxwellian populations is described by the sum over all particle species in the plasma,

 ∑sCis{fi} = 1τie∑snsZ2sne{ϕ(xs)−G(xs)2x3i∂∂ξ[(1−ξ2)∂fi∂ξ] (3) +1x2i∂∂xi[2TiTsx2iG(xs)fi+xiG(xs)∂fi∂xi]} +1τii[2νs(v)xiξu∥+νE(v)x2iQ]fMi,

where , and the usual error function and Chandrasekhar function appear. The dimensionless moments and of the distribution function appearing in the momentum and energy restoring terms of the self-collision operator are

 Missing dimension or its units for \hskip (4)

with the scattering frequencies given by

 νs(v)=4G(xi)xi,νE(v)=2(4G(xi)xi−ϕ(xi)x3i). (5)

The runaway behavior of interest can be demonstrated by considering the simpler dynamics of an ion test particle (8); (19) in the presence of the electric field. The ion equation of motion takes the form , where the test particle collisional friction is given by

 Ftesti(v)=−Z2ieED∑snsZ2sneTeTs(1+msmi)G(xs), (6)

and is the Dreicer field. Thus ion acceleration can occur when . The required electric field values are illustrated in Fig. 1, for low and high effective charge. The figures illustrate the non-monotonic ion friction force, with one maximum near the thermal velocities of the ion species, and another near the electron thermal velocity.

Therefore, as first described by Furth and Rutherford (12), for a sufficiently strong electric field we may expect ions to be accelerated from their initial velocity to a higher velocity at which friction on the electrons dominates, giving rise to a suprathermal population in the plasma. We will compare the characteristics of the test-particle behavior, governed by the friction force illustrated in Fig. 1, to the numerical solution of Eq. (2) in the next section.

Figure 1 also illustrates that the electric fields needed to accelerate ions are highly sensitive to ion charge and plasma composition, due to their effect on the effective electric field . Note that the electric fields needed to accelerate ions beyond the electron thermal velocity are significantly larger than the minimum electric field necessary for acceleration. Such strong fields will not be considered here, since the validity of the linearization typically breaks down before the ions reach a significant fraction of the electron thermal velocity.

## Iii Codion

In this section, we outline the implementation of Eq. (2) in the numerical tool CODION, which solves the ion Fokker-Planck equation numerically as an initial value problem to give the evolution of the ion distribution function in the presence of an accelerating electric field. It uses a continuum-spectral discretization scheme based on that used in CODE (20). Illustrative solutions for a tokamak-like plasma are presented, and a comparison of the obtained distribution function is made with the behavior predicted by test-particle equations, demonstrating the importance of collisional diffusion for the runaway of ions. In addition we investigate the effect of different choices for the self-collision field-particle operator.

The pitch-angle dependence of the distribution function is represented by a truncated Legendre polynomial expansion, while velocity is discretized on a uniform grid , :

 fi(vn,ξ,t)=lmax∑l=0fl(vn,t)Pl(ξ), (7)

where the Legendre polynomials obey the orthogonality relation , and

 fl(vn,t)=2l+12∫1−1dξPl(ξ)fi(vn,ξ,t). (8)

The integral operation is applied to the kinetic equation in Eq. (2) for each , producing a linear set of equations for the quantities , using the boundary condition for all . Well-known recurrence relations for the Legendre polynomials are used to obtain analytic expressions for all the terms appearing in the equation. For example, the Legendre polynomials are eigenfunctions of the linearized collision operator, while the electric field-term will produce a coupling between and modes. This procedure exactly captures number conservation for any choice of . The velocity derivatives appearing in the kinetic equation are represented with five-point stencils

 ∂fl∂v∣∣vn =112ΔvN−1∑m=0(−δn,m−2+8δn,m−1−8δn,m+1+δn,m+2)fl(vm), (9) ∂2fl∂v2∣∣vn =112Δv2N−1∑m=0(−δn,m−2+16δn,m−1−30δn,m+16δn,m+1−δn,m+2)fl(vm), (10)

formally introducing an error of order . The integral moments of the ion distribution appearing in the self-collision model operator are discretized with a quadrature of the form , where the quadrature weights are chosen according to Simpson’s rule, also with error of order .

For the distribution function to be single-valued and smooth at the origin, we enforce the boundary condition for all and . Since we restrict ourselves to cases where electron friction will dominate the electric field at sufficiently high velocities, the maximum resolved velocity can always be chosen so that only insignificant numbers of particles are near the edge of the grid, minimizing the effect of the choice of boundary condition. We use the Dirichlet boundary condition for all at the maximum velocity. A detailed investigation of the convergence properties of the solutions with respect to discretization parameters is described in Ref. (21).

The discretization procedure outlined above casts the kinetic equation, Eq. (2), into the form

 ∂F∂t+MF=0, (11)

where is a sparse matrix and is a vector representing the discretized distribution function. Time integration is performed with the first order implicit scheme

 F(tn+1)=[I+ΔtM(tn+1)]−1F(tn), (12)

where any time-dependence of the operator is due to time-variation of electric fields and background plasma parameters. An arbitrary plasma composition is determined by a set of input vectors containing particle masses , corresponding charge states , charge densities and temperatures .

Figure 2(a) shows a typical example of the evolution of the ion distribution for a case where the electric field is above the minimum required for runaway acceleration. The plasma parameters used are characteristic for tokamaks with a hot bulk deuterium plasma at keV and fully ionized native carbon impurities. Note that the loop voltage is typically  V in normal tokamak operation, corresponding to  V/m (2); (9). A contour plot of the distribution in velocity space when steady state is reached is shown in Fig. 2(b). For the parameters used here, approximately of the ion population has been accelerated and the linearization used to derive Eq. (2) is well satisfied.

The steady state distribution is typically established in 10-20 s at this temperature and density, and the time-scale varies with plasma parameters as the collision time defined in connection with Eq. (3), . For stronger electric fields, the initial evolution of the distribution can be followed, but the linearization breaks down before the steady state can be reached. The entire ion distribution will eventually run away when for the case considered here, and the linearization breaks down within ms with such an electric field.

## Iv Results

Expansions of the collision operators appearing in Eq. (2) have been used previously in the literature to consider ion runaway analytically. Here we compare the effect of various approximations to the collision operator on the numerical solution of the distribution function. The plasma parameters of Fig. 2 are used as a basis in the comparisons, but the behavior illustrated is characteristic of a wide range of parameters. We first consider the effect of neglecting the conservative field-particle terms in the self-collision operator. Figure 3 shows the cut of the distribution of the bulk deuterium population for two cases, with effective charge and , respectively. It can be seen that, using only the test-particle operator, the dominant behavior of the fast-ion population as given by the fully conserving case is reproduced, with the main difference being the rate at which it builds up. This is expected, since the conserving field-particle terms are proportional to , and therefore act mainly on the thermal bulk of the distribution.

Figure 3(a) shows a case with . With a lower amount of impurities to which momentum and energy can be transferred, the fully conserving linearized collision operator for self-collisions will exhibit unphysical behavior before a significant runaway population has time to form, which is clearly illustrated by the distortion near . The reason is that when the conserving terms are kept in the kinetic equation, the distribution is heated by the electric field, causing the linearized equation to break down after some time. This is also observed in the solution obtained using only the energy conserving term, albeit less pronounced. The distribution functions obtained with the momentum conserving self-collision operator tends to stay regular for longer. Figure 3(b) shows a similar case but with higher impurity content, , for which all operators yield well-behaving results. The main difference between the conserving operators and only using the test-particle operator is that the former typically lead to a runaway rate which is at least twice as large. An additional consequence of the low effective charge was demonstrated in Fig. 1(a), where impurity ions were shown to be more easily accelerated by an electric field than the bulk species, implying that for low the assumption of stationary impurity ions may be violated.

In conclusion, the high-energy part of the ion distribution obtained using only the test-particle operator is in qualitative agreement with the result obtained with conservative operators, but the runaway rate is expected to be lower in the test-particle case. A quantitative investigation of runaway rates for impurity species is presented in Section V.

It is instructive to compare the behavior of the numerical solution to the characteristic behavior indicated by the test particle friction given in Eq. (6). Noting that the runaway ion velocity satisfies , we can expand the contributions to Eq. (6) using the known low and high-velocity limits of the Chandrasekhar function. The test particle friction in this limit reduces to

 Ftesti≈−mivTiτie[Zeff+¯¯¯nxi+43√π(TiTe)3/2√memixi], (13)

where allows for arbitrary impurity content. Consider first the minimum value of the magnitude of the collisional friction force; this will determine the minimum electric field which can accelerate a fast test ion. Minimizing Eq. (13) yields

 vmin =vTe[3√π2memi(Z% eff+¯n)]1/3, (14) Ftesti(vmin) =−2mivTiτieTiTe[32πmemi(Zeff+¯n)]1/3. (15)

From this it follows that the minimum, ”critical”, value of the electric field above which a test ion can be accelerated is given by

 Missing or unrecognized delimiter for \right (16)

By taking , we can find the range of test ion velocities, , for which acceleration in a given electric field occurs, as discussed in Refs. (1); (12). Using the expression for the friction given in Eq. (13) results in a third order equation, however simpler approximate formulae can be obtained by noting that will fall near to the region dominated by ion friction, and in the region dominated by electron friction. Retaining only the corresponding terms in Eq. (13) yields, for arbitrary impurity content,

 vc1 = vTi√ZiTe2Ti(EED)−1Zeff+¯¯¯n|1−Zi/Zeff|, (17) vc2 = vTe3√π2EED|1−Zi/Zeff|Zi. (18)

These equations generalize the corresponding expressions in Ref. (1) to arbitrary plasma composition. Note that these formulae are only valid when is sufficiently large compared to , since at ion and electron friction contribute equally. We may expect that in steady-state, the position of the high-velocity maximum of the distribution function, denoted , is close to , which scales linearly with in the approximate form given by Eq. (18). This is confirmed numerically and illustrated in Fig. 4, where we show the variation with electric field of , obtained from steady-state CODION solutions of Eq. (2). Also shown are the boundary velocities of the acceleration region, resulting from numerical solution of the force balance using the full test particle friction, Eq. (6), as well as their approximate forms Eqs. (17-18). The values converge when the system is strongly driven by a large . The linear dependence of is clearly visible at large , where it approaches the value given by the test particle approximation. The analytic approximation for , Eq. (16), is only accurate to , however, indicating that collisional diffusion contributes significantly to the evolution at lower electric fields. Since the linearization breaks down more rapidly at larger electric fields, it is mainly at fields near the threshold for runaway generation that the model can consistently be applied to study the long-term evolution of the ion runaway tail, making a full kinetic simulation essential for capturing the important physics. For the more massive impurities, the features of the test-particle model become increasingly accurate since the runaway ion population is further separated in velocity space from the thermal bulk.

It is important to point out that neither the diffusion terms nor the field-particle self-collisions have been accounted for in the derivations of the above estimates, which are meant to give simple expressions that show how the essential quantities scale with the plasma parameters, and to provide a useful physical picture for illustrating the ion runaway phenomenon. A complete description will be provided only by numerical solution of the kinetic equation.

## V Applications

In this section, CODION is applied to calculate runaway ion distributions for typical solar flare and fusion plasmas. The rate at which a fast ion population forms due to the runaway mechanism is determined, and it is investigated whether the difference in acceleration rate between different ion species can explain the enhanced abundance of heavy ions in the solar wind. We also consider the possibility of Alfvénic instabilities being driven by runaway ions during tokamak disruptions.

Throughout this section, time-scales are chosen so that significant fast ion populations have time to form, which typically takes a few hundred ion-electron collision times. We define the runaway density, , as the fraction of ions with velocity larger than the low-velocity numerical solution of , denoted , which if is taken as the velocity minimizing the friction .

### v.1 Ion acceleration in solar flares

Solar flares are thought to be initiated by reconnection in the corona (22), but the origin of the observed fast ion populations in flares is still not completely understood (23). Both stochastic acceleration by waves and the direct acceleration of the particles by the electric field have been considered, and it appears likely that a combination of the two can be at work (1); (24); (25).

The effective accelerating field experienced by a given species varies with its charge and the effective charge of the plasma, as discussed in Section II. This can give rise to preferential acceleration of heavier elements under certain circumstances, and this effect was considered in (1), where estimates of the runaway rate were given based on the approximate formula of (11). With CODION we can determine the time evolution of the ion distribution function numerically, and evaluate the dependence of the acceleration rate on various ion parameters.

The composition of the solar plasma, particularly the metallic elements, has been studied extensively in recent years, however much uncertainty remains. We will choose parameters consistent with the choices made by Holman (1). We use the plasma temperature eV for all particle species, and hydrogen density m. Elements with atomic number can readily be assumed to be fully ionized at this temperature. The plasma composition is based on the ion abundance recommended by Schmeltz et al. (13). We use a helium population of density , and represent all heavier impurities by a carbon population of density , corresponding to an effective charge . Electric field strengths in solar flares are not well constrained by experimental observation, and we will investigate the rate of acceleration at a range of values. The Dreicer field is mV/m for this set of plasma parameters.

The critical electric fields for the ion species in such a plasma are  mV/m for hydrogen,  mV/m for helium and mV/m for carbon. Note that the acceleration rate depends strongly not only on , but also on . Figure 5 shows the cut of the distribution functions of hydrogen (H), helium (He) and carbon (C) after 30 s of acceleration from initial Maxwellians, with the plasma parameters specified above and an electric field mV/m. This is significantly below the hydrogen critical field, and no hydrogen runaway population forms. Runaway ion populations of both helium and carbon do form however, with runaway densities and , respectively.

Positive values of represent the direction of the electric field. Therefore, Fig. 5 also illustrates how heavier ions (charge ) are accelerated in the direction opposite to the electric field, dragged by electron friction. The corresponding 2D carbon distribution function is shown in Fig. 6, displaying a strong directional anisotropy (compare to Fig. 2(b)). This can be understood by the observation that the accumulation velocity is located at a higher value of for heavier ions. Since pitch angle scattering of the energetic heavy ions scales with velocity like , the mechanism will be less effective than for light bulk ion species in increasing the perpendicular energy of the distribution.

We will now investigate how the rate of acceleration varies with electric field strength for different ion species. To determine the runaway density of heavier ions, we have introduced trace amounts of each ion species with charge between 2 (helium) and 18 (argon), assumed to be fully ionized. In practice, ions of charge will typically not be fully ionized at the temperature considered due to their high ionization energy, meaning that the results shown here will overestimate the acceleration rate of the heavier ions. The ion masses have been set to that of the most common isotope, i.e. Li, Be, Ne etc. Both He and He are shown, with He showing a significantly enhanced runaway rate as compared to that of He, for all values of .

Figure 7(a) shows how the runaway density depends on electric field after 1 s of acceleration for various ion species present in the solar flare plasma. The figure illustrates how the runaway rate is sensitive to ion parameters; at low electric fields the heavier ions tend to be accelerated at a lower rate than light ions, while at higher fields they are the most readily accelerated. Note that above the critical electric field, mV/m for helium, the runaway rate increases significantly faster than for . The runaway rate of He is seen to be orders of magnitude higher than the runaway rate of He for all electric fields considered.

Finally we illustrate the dependence of the acceleration rate on ion charge and mass. Figure 7(b) shows the runaway density after 1 s of acceleration as a function of ion charge for various electric fields. Ions of charge between (Be) and (O) are seen to be preferentially accelerated over lighter or heavier elements for low electric fields. For , the trend depends on electric field. For low electric fields, the runaway rate decreases with charge, while for larger electric fields (in this case above approximately 50 mV/m), heavier ion species may be accelerated at a higher rate than the lighter species. Further studies are required to determine the full effects of variation in background composition, temperature, density, and charge states on the relative rates of acceleration between different ion species. The results presented here demonstrate the utility of CODION for the problem.

As previously noted, acceleration by quasi-static electric fields is not the only mechanism for ion acceleration in a flare plasma. Interaction with Alfvén waves can accelerate ions which have velocities above the Alfvén velocity. As this is usually well above the thermal ion velocity, an initial acceleration by electric fields may be required before the process becomes significant (26); (1). CODION provides the means for more accurate modeling of the effects of such interactions.

### v.2 Tokamak disruptions

During tokamak disruptions the plasma temperature drops from the typical operating regime of around several keV to a few eV in a couple of milliseconds. A large electric field is initially induced parallel to the magnetic field to maintain the plasma current of several MA, potentially leading to the formation of a beam of energetic electrons through the runaway mechanism. The potential for damage by such a focused high-energy beam on contact with the vessel wall is large, and runaway generation must as far as possible be suppressed. To study the physics and mitigation of runaway electrons, disruptions can be induced by the injection of large quantities of noble gas, often in amounts comparable to the initial plasma inventory or larger (2).

The large induced electric field will usually decay rapidly on a timescale of a few ms in response to the formation of a narrow runaway electron (RE) beam. With runaway electrons reaching energies of order tens of MeV, they can carry a significant fraction of the pre-disruption plasma current and can drive high frequency electromagnetic instabilities through resonant interactions (27); (28); (29); (30). Recently, low-frequency magnetic fluctuations in the range kHz, have been observed in the TEXTOR tokamak during induced-disruption studies with argon massive gas injection (MGI). These fluctuations take the form of either a strong signal at a distinct frequency (31), or accompanied by broadband activity (32). The fluctuations appear to limit the RE beam formation in these cases, as the magnetic perturbations may scatter the runaway electrons and provide passive mitigation. Aside from the potential consequences for mitigation, observed instabilities offer a non-intrusive diagnostic for both bulk plasma and fast-particle properties, through the extensively applied technique of MHD spectroscopy (33).

Fast ions resulting, for example, from certain heating schemes are well known to resonantly drive low frequency Alfvénic instabilities in typical operational scenarios (34); (35). Runaway ions may thus also provide a potential drive for the fluctuations observed. Interestingly, toroidal Alfvén eigenmodes (TAEs) (36) can have frequencies and mode numbers in the same range as the post-disruption magnetic activity. Therefore, the excitation of TAEs by runaway ions was recently considered in Ref. (14) using an analytical approximation for the runaway distribution function. As the validity of the approximate distribution was limited, definite conclusions could not be drawn. With CODION we can extend the study using the numerically calculated ion distribution. The cold post-disruption plasma is highly collisional, motivating the use of a homogeneous background plasma and the neglect of magnetic trapping when evaluating the effective electric field.

The TAE perturbation is typically dominated by two neighboring toroidally coupled harmonics at large aspect ratio, with poloidal mode numbers and . The mode is localized about the minor radius at which the magnetic safety factor is , where is the radius of the magnetic axis and is the toroidal mode number. The TAE frequency is , where is the Alfvén speed and the mass density. The two component harmonics allow resonant interaction with particles whose parallel velocity satisfies or . It was argued in Ref. (14) that as the runaway ions accelerate, the inverted region of their energy distribution, where is the particle energy, can reach the lower Alfvén resonance, and may drive the TAE. If the radial runaway ion profile peaks on axis, the spatial gradient will give an additional positive contribution to the growth rate. Taking parameters characteristic of argon MGI-induced disruptions, m, , , and anticipating a native background carbon impurity with and so that , a background ion temperature of , toroidal magnetic field T and major radius m means that the resonance condition requires deuterium ions with velocities .

The electric field required to accelerate bulk ions to the resonant velocity at these low temperatures is substantial, varying in response to changes in but is typically V/m. Such field strengths are unlikely to occur during a disruption, and they would be short-lived if they did (37). Therefore, we conclude that whilst ion runaway may be of interest in hot fusion plasmas, runaway ions are unlikely to provide the drive for the observed fluctuations during disruptions.

To quantify the electric field needed for significant ion runaway, we show in Fig. 8 how the deuterium distribution evaluated after ms of acceleration from an initial Maxwellian – a typical time scale for the induced electric field – varies with (a constant) electric field. The parameters are m,  eV and the same plasma composition as before with . It is seen that for electric fields below  V/m, no runaway tail tends to form, and even with  V/m the fast ions are far from the resonant velocity near . The behavior is sensitive to which temperature is chosen for the plasma: increased temperature decreases the electric fields needed to accelerate ions, but makes the acceleration timescale longer.

Note that a higher amount of assimilated argon, or a weaker magnetic field, would lead to a lower Alfvén velocity and TAE frequency, and therefore the runaway ions would reach the resonance condition more easily. The higher electron density leads to an increased collision frequency, however, so that the runaway ion distribution requires a longer time to form. The level of argon absorption during mitigation varies between machines (2) and the absorption appears to decrease for bigger machines. Therefore, although we may expect a higher electric field in large tokamaks such as ITER, the higher bulk plasma density and likely lower absorption suggests no ion runaway will occur. Note that as the ITER magnetic field is so high, even if the assimilated argon was equal to the initial deuterium inventory, the ions would need to reach (at 10 eV) to reach the typical resonance, assuming T, m and m.

Finally we note that the model presented in this paper assumes the initial ion distribution to be a stationary Maxwellian. Fast ion populations present due to heating schemes in use before the disruption may not be completely expelled, and it is not certain that the initial distribution is accurately described by a Maxwellian. The dynamics of ion acceleration starting from these non-Maxwellian distributions may yield a fast-ion population more readily. Furthermore, since plasma conditions change on a time scale of a few collision times during the sudden cooling associated with instabilities or disruptions, so-called “hot-tail” runaway generation may occur. This has been shown to be an important effect for runaway of electrons, where a seed of runaway electrons are provided by fast electrons present before cooling (38). These fast electrons are cooled at a slower rate than the low-energy electrons, and may find themselves in the runaway region when the plasma has reached its final temperature. However, using CODION to investigate the effect of hot-tail ions, it has been concluded that the effect is small for realistic fields. The reason is that, if the electric field is not high enough for the ions to overcome the friction and become runaways in the first place, all hot-tail ions will necessarily also be slowed down. A low-velocity inverted ion population may however form as a result of the cooling process even for such low electric fields, with its peak near the velocity which minimizes the collisional friction at the final temperature (typically around 6-8 thermal velocities). This is still significantly lower than that needed for resonant interaction with Alfvén waves.

One refinement to the model would be the inclusion of “knock-ons”, i.e. large-angle collisions, which have been neglected in the Fokker-Planck equation. It is well known that single collisions can change the momenta of the interacting particles significantly, and a runaway ion interacting with a bulk ion could cause both to end up in the runaway region. In a situation where the electric field is low enough that runaway ions are produced at a low rate through the standard acceleration mechanism, knock-on collisions could possibly contribute significantly to the runaway generation rate. This has been demonstrated to be the case for electron runaway, where this effect drastically affects the rate at which runaways are produced. A simplified runaway ion knock-on operator could potentially be constructed from the Boltzmann collision integral under the assumption that fast ions accumulate near and collide mainly with the bulk distribution, since the fast ion distribution is assumed to be a small perturbation in our linearized model. However, there are differences between ion runaway and electron runaway that suggest that knock-on runaway generation is a less significant effect for ions than for electrons. Since our linearized model restricts the study to cases where (which is also the regime where knock-on generation would be expected to be significant), the accumulation velocity near will not be significantly larger than the runaway velocity . Because of this, collision events where both particles end up in the runaway region will be less frequent. This is in contrast to electron runaway, where the electrons have unbounded energy (neglecting radiation effects).

There are applications for knock-on operators other than avalanche generation. It has been suggested that fast ion populations due to other sources – for example hot alpha particles created in fusion reactions or ions heated by external sources such as neutral beam injection (NBI) or radio frequency (RF)-heating – could accelerate bulk impurity ions, which could in turn be used for diagnostics (39); (40). The suggested collision operator could be implemented in CODION, and the time-evolution of bulk impurities solved for in the vicinity of an assumed or numerically obtained background of fast ions, however this is outside the scope of the present paper.

## Vi Conclusions

Electron runaway resulting from the occurrence of a strong electric field in a plasma has been the subject of extensive study, and numerical tools exist to simulate the electron dynamics. The analytic description of the associated ion acceleration was developed at the same time, but its application has been much more limited and is restricted by the various approximations to the collision operator which were required.

We have developed an efficient open-source numerical tool, CODION (41), which solves the ion Fokker-Planck equation as an initial value problem in a fully ionized plasma. A uniform background magnetic field is assumed, along with initially stationary Maxwellian distributions, however arbitrary impurity densities and temperatures may be specified. A model operator for ion self-collisions based on that used in the gyrokinetic code GS2 (42) has been employed, satisfying momentum and energy conservation, non-negative entropy production and self-adjointness. A simplified analytical model based on the large mass ratio is used for ion-electron collisions, allowing a description of ion-electron friction caused by the perturbation of the electron distribution due to the electric field. However, we wish to note that our model will break down for strong electric fields, as the electrons – which are assumed to be in force balance with the electric field and ion friction – will be rapidly accelerated by electric fields approaching the Dreicer field. A full description of such scenarios would require the simultaneous evolution of the ion and electron distributions, for example by coupling the CODION and CODE (20) codes.

The effect of various approximations to the collision operator commonly used in the literature has been studied numerically. It has been demonstrated that the addition of momentum and energy conservation in self-collisions mainly acts to increase the rate at which the fast ion population builds up, while the qualitative behavior is largely unaffected. For strong electric fields, the test particle description is seen to reproduce the characteristic velocity achieved by the runaway population well, and the predicted critical electric field for ion acceleration is accurate to . Using the test particle approximation, we derived concise analytic expressions for the critical electric field for ion runaway, Eq. (16), and the typical runaway energy, Eqs. (17-18), and tested these expressions against direct numerical simulation with CODION.

The output of CODION is the evolution of the 2D velocity space ion distribution. The utility of this has been demonstrated for calculating acceleration rates of ions in solar flare plasmas. The rate at which ions are accelerated has been evaluated for a range of ion masses and charges for a solar flare scenario based on that considered by Holman (1), and an exponential dependence of acceleration rate with charge for has been illustrated for this scenario.

Low-frequency instabilities, in the range characteristic of TAE modes, have been observed in post-disruption tokamak plasmas, where the disruption was induced by massive gas injection. Using CODION we have considered the potential for ions accelerated in the disruption-induced electric field to drive such modes resonantly. The post-disruption discharge parameters are not well constrained experimentally and simulations of a range of values indicate that ion acceleration is possible. The typical maximum ion velocity achieved is too low for resonant interaction to occur, however, and the rate of runaway generation is too slow for a significant runaway density to be reached in the short-lived electric fields of a typical disruption.

## Acknowledgements

The authors are grateful to Gergely Papp, Matt Landreman, István Pusztai and Joan Decker for fruitful discussions. This project has received funding from the Knut and Alice Wallenberg Foundation and the RCUK Energy Programme [under grant EP/I501045].

### References

1. G. D. Holman, Astrophys. J. 452, 451 (1995).
2. E. M. Hollmann, M. E. Austin, J. A. Boedo, N. H. Brooks, N. Commaux, N. W. Eidietis, D. A. Humphreys, V. A. Izzo, A. N. James, T. C. Jernigan, A. Loarte, J. Martin-Solis, R. A. Moyer, J. M. Muñoz-Burgos, P. B. Parks, D. L. Rudakov, E. J. Strait, C. Tsui, M. A. Van Zeeland, J. C. Wesley and J. H. Yu, Nucl. Fusion 53, 083004 (2013).
3. R. M. Kulsrud, Y.-C. Sun, N. K. Winsor and H. A. Fallon, Phys. Rev. Lett. 31, 690 (1973).
4. M. Lehnen, A. Alonso, G. Arnoux, N. Baumgarten, S. A. Bozhenkov, S. Brezinsek, M. Brix, T. Eich, S. N. Gerasimov, A. Huber, S. Jachmich, U. Kruezi, P. D. Morgan, V. V. Plyusnin, C. Reux, V. Riccardo, G. Sergienko, M. F. Stamp and JET EFDA Contributors, Nucl. Fusion 51, 123010 (2011).
5. M. D. Kruskal and I. B. Bernstein, Princeton Plasma Physics Lab. Report No. MATT-Q-20 (1962) (unpublished).
6. J. W. Connor and R. J. Hastie, Nucl. Fusion 15, 415 (1975).
7. M. N. Rosenbluth and S. V. Putvinskii, Nucl. Fusion 37, 1355 (1997).
8. T. Fülöp and M. Landreman, Phys. Rev. Lett. 111, 015006 (2013).
9. P. Helander, L.-G. Eriksson, R. J. Akers, C. Byrom, C. G. Gimblett and M. R. Tournianski, Phys. Rev. Lett. 89, 235002 (2002).
10. S. Eilerman, J. K. Anderson, J. S. Sarff, C. B. Forest, J. A. Reusch, M. D. Nornberg and J. Kim, Phys. Plasmas 22 020702 (2015).
11. A. V. Gurevich, Sov. Phys. JETP 13, 1282 (1961).
12. H. P. Furth and P. H. Rutherford, Phys. Rev. Lett. 28, 545 (1972).
13. J. T. Schmelz, D. V. Reames, R. von Steiger and S. Basu, Astrophys. J. 755, 33 (2012).
14. T. Fülöp and S. Newton, Phys. Plasmas 21, 080702 (2014).
15. I. G. Abel, M. Barnes, S. C. Cowley, W. Dorland and A. A. Schekochihin, Phys. Plasmas 15, 122509 (2008).
16. A. Stahl, E. Hirvijoki, J. Decker, O. Embréus and T. Fülöp, to appear in Phys. Rev. Lett. (2015).
17. H. Dreicer, Phys. Rev. 115, 2 (1959).
18. P. Helander and D. J. Sigmar, Collisional Transport in Magnetized Plasmas (Cambridge University Press, Cambridge, 2002).
19. F. L. Hinton, Collisional Transport in Plasma in Handbook of Plasma Physics Vol. 1 (North-Holland Publishing Company, 1983).
20. M. Landreman, A. Stahl and T. Fülöp, Comp. Phys. Comm. 185, 847 (2014).
21. O. Embréus, Master’s thesis, Chalmers University of Technology (2014) http://publications.lib.chalmers.se/records/fulltext/210276/210276.pdf.
22. C. Liu and H. Wang, Astrophys. J. 696, 121 (2009).
23. T. L. Garrard and E. C. Stone, Adv. Space Res. 14, 589 (1994).
24. D. Braddy, J. Chan and P. B. Price, Phys. Rev. Lett. 30, 669 (1973).
25. B. Sylwester, J. Sylwester, K. J. H. Phillips, A. Kȩpa and T. Mrozek, Astrophys. J. 787, 122 (2014).
26. E. R. Harrison, J. Nucl. Energy Part C: Plasma Physics 1, 105 (1960).
27. T. Fülöp, G. Pokol, P. Helander and M. Lisak, Phys. Plasmas 13, 062506 (2006).
28. G. Pokol, T. Fülöp and M. Lisak, Plasma Phys. Control. Fusion 50, 045003 (2008).
29. A. Kómár, G. I. Pokol and T. Fülöp, Phys. Plasmas 20, 012117 (2013).
30. G. I. Pokol, A. Kómár, A. Budai, A. Stahl and T. Fülöp, Phys. Plasmas 21, 102503 (2014).
31. G. Papp, Ph. W. Lauber, M. Schneller, S. Braun, H. R. Koslowski, TEXTOR Team, ECA, Vol. 38F (2014) http://ocs.ciemat.es/EPS2014PAP/pdf/P2.032.pdf
32. L. Zeng et al., Phys. Rev. Lett. 110, 235003 (2013).
33. H. A. Holties, A. Fasoli, J. P. Goedbloed, G. T. A. Huysmans and W. Kerner, Phys. Plasmas 4, 709 (1997).
34. W. W. Heidbrink, Phys. Plasmas 15, 055501 (2008).
35. F. Zonca et al., Nucl. Fusion 47, 1588 (2007).
36. C. Z. Cheng and M. S. Chance, Phys. Fluids 29, 3695 (1986).
37. H. Smith, P. Helander, L.-G. Eriksson, D. Anderson, M. Lisak and F. Andersson, Phys. Plasmas 13, 102502, (2006).
38. P. Helander, H. Smith, T. Fülöp and L.-G. Eriksson, Phys. Plasmas 11, 5704 (2004).
39. V. G. Nesenevich et al., Plasma Phys. Control. Fusion 56, 12 (2014).
40. P. Helander, M. Lisak and D. D. Ryutov, Plasma Phys. Control. Fusion, 35, 3 (1993).
41. CODION is available at https://github.com/Embreus/CODION.
42. M. Kotschenreuther, G. Rewoldt and W. M. Tang, Comp. Phys. Comm. 88, 128 (1995).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters