Assessing energy dependence of the transport of relativistic electrons in perturbed magnetic fields with orbit-following simulations
Konsta Särkimäki, Ola Embreus, Eric Nardon, Tünde Fülöp, JET Contributors
AAssessing energy dependence of the transport ofrelativistic electrons in perturbed magnetic fieldswith orbit-following simulations
Konsta S¨arkim¨aki , Ola Embreus , Eric Nardon , T¨undeF¨ul¨op , and JET Contributors * Department of Physics, Chalmers University of Technology, SE-41296G¨oteborg, Sweden Association EURATOM-CEA, IRFM, CEA Cadarache, 13108St-Paul-lez-Durance, France * See the author list of E. Joffrin et al. 2019 Nucl. Fusion 59 112021E-mail: [email protected]
Abstract.
Experimental observations, as well as theoretical predictions,indicate that the transport of energetic electrons decreases with energy. Thisreduction in transport is attributed to finite orbit width (FOW) effects. Usingorbit-following simulations in perturbed tokamak magnetic fields that have anideal homogeneous stochastic layer at the edge, we quantify the energy dependenceof energetic electrons transport and confirm previous theoretical estimates.However, using magnetic configurations characteristic of JET disruptions, we findno reduction in RE transport at higher energies, which we attribute to the modewidths being comparable to the minor radius, making the FOW effects negligible.Instead, the presence of islands and nonuniform magnetic perturbations are foundto be more important. The diffusive-advective transport coefficients calculated inthis work, based on simulations for electron energies 10 keV – 100 MeV, can beused in reduced kinetic models to account for the transport due to the magneticfield perturbations.
Keywords : runaway electrons, stochastic magnetic field, transport, plasma disruption,orbit-following a r X i v : . [ phy s i c s . p l a s m - ph ] J un ssessing energy dependence of the transport of relativistic electrons
1. Introduction
Runaway electrons (REs) generated in disruptions area major concern for future tokamaks with large plasmacurrents, such as ITER [1]. Predictions show thata large fraction of the initial plasma current can beconverted into a runaway electron beam with energiesof tens of MeV due to the avalanche process [2]. Thesubsequent uncontrolled loss of the RE current coulddamage the plasma facing components and has to beavoided.The growth of the runaway electron populationcan be hindered through mechanisms that lead toradial losses. Substantial radial losses can reduce boththe seed runaway population and the avalanche growthrate significantly [3].Losses due to magnetic perturbations are expectedto be present naturally in the early phase of thedisruption, when the topology of the magnetic fieldconfining the particles undergoes drastic changes,including the formation of stochastic regions leadingto rapid radial transport of the REs. Magneticfield perturbations can also be induced by externalfield coils. Such perturbations have been tested inseveral tokamaks, however, while runaways could besuppressed in medium-size tokamaks [4, 5, 6], theperturbations showed no significant effect on runawayson JET [7]. The main focus of this paper is toquantify the radial losses of energetic electrons causedby magnetic perturbations in disruption scenarios.Charged particle transport in a toroidal magneticgeometry with broken flux surfaces is often assumedto be diffusive with a diffusion coefficient thatscales quadratically with the radial magnetic fieldperturbation amplitude δB [8, 9]: D ≈ √ v (cid:107) λ (cid:107) (cid:18) δBB (cid:19) . (1)This type of diffusion arises when particles are movingalong the stochastic field lines with velocity v (cid:107) andtheir orbits become decorrelated after travelling adistance λ (cid:107) , the so-called parallel correlation length ofthe magnetic field.However, Eq. (1) has been found to overestimatethe transport of REs (see the discussion in Ref. [10]).Two possible explanations have been put forward: theenergy dependence of the transport is different fromthat implied by Eq. (1), such that the transport isreduced for more energetic particles, or, the magneticfield is not completely ergodic and the remainingmagnetic islands hinder the transport [11]. In this workwe focus on the former issue—the energy dependenceof the transport.The non-trivial energy dependency of the diffusioncoefficient is introduced via a scaling factor, Υ, [10, 12] as D ≈ √ v (cid:107) λ (cid:107) (cid:18) δBB (cid:19) Υ . (2)Physically, the scaling factor accounts for the finiteorbit width (FOW) effects that arise when energeticparticles drift from their initial flux surface, whichinterfere with the decorrelation process and lead toreduced transport, Υ ≤
1. To clarify, FOW effectsrefer solely to the processes that are present even in astochastic field with uniform structure; in a nonuniformfield there could be additional mechanisms that areenergy dependent e.g. electrons with high energy beingable to make excursions to a non-ergodic region. Theexisting theoretical work does not factor these othereffects in Υ.The exact form and value of Υ depends onthe magnetic perturbations via the parallel and theperpendicular magnetic field correlation lengths, λ (cid:107) and λ ⊥ , and on the electron energy through thegyroradius and the orbit width: ρ g = γmv ⊥ eB , (3) d orb = q γmveB . (4)Here q is the safety factor and γ is the Lorentz factor.For electrons, FOW effects are expected to becomeimportant only at relativistic energies, γ (cid:29)
1. Theexplicit form of Υ is discussed in section 2 of this paper.Note that the scenarios that were studied inRefs. [10, 13, 14] were turbulent flat-top scenarios,i.e. the level of magnetic perturbations was consider-ably lower than the one typical of the thermal quenchof a disruptive plasma. Although the theoretical esti-mates relevant to flat-top scenarios are not expectedto hold in disruptions, they can be used to benchmarkthe numerical models.The energy dependence of the RE transportcan also be studied numerically using orbit-followingsimulations. The benefit of such simulations is thatthey can accurately model motion of particles and theirtransport due to magnetic perturbations using realisticmagnetic fields, if those can be provided. In this work,we use orbit-following simulations to: • Assess whether significant reduction in REtransport occurs (i.e. Υ <
1) at high energyin disruption magnetic fields, or when externalcoils are used for RE mitigation, for plausible REenergies (section 3). • Verify the theoretical estimate for the energy-dependence of the transport given in Ref. [10](section 4).We quantify the cross-field transport of REs by com-puting the corresponding (energy-dependent) advec-tion and diffusion coefficients. The advantage of this ssessing energy dependence of the transport of relativistic electrons No orbit width effects O r b i t a v e r a g i n g O r b i t d e c o r r e l a t i o n log Energy T r a n s p o r t Figure 1.
Illustration of the expected energy dependence ofthe transport. Here λ (cid:48)⊥ > λ , and the conditions d orb > λ (cid:48)⊥ and τ ⊥ < τ (cid:107) are met simultaneously. The dotted lines include finitegyro-radius effects. method is that the transport coefficients enable one tointroduce 3D magnetic field effects to reduced kineticmodels, which makes it possible to assess the reductionin avalanche growth rate due to the field stochasticity.However, this task is beyond the scope of the presentwork.
2. Transport in the presence of finiteorbit-width effects
Figure 1 illustrates the expected energy dependence ofthe radial transport. After increasing initially due toincreasing velocity, the transport ceases to grow when v ≈ c , where it would remain flat if FOW effects werenot taken into account. At higher energy, the exactbehaviour depends on two mechanisms, orbit-averaging and perpendicular decorrelation , that are responsiblefor FOW effects. These mechanisms are presented indetail in Ref. [10] and here we only review the mainpoints.Particles with finite orbits do not trace the samefield line exactly due to the radial drift. This oscillationin radius is accompanied by toroidal precession, whichis the transit-time averaged drift toroidally. Whilethe particle returns to same radial position after eachtransit, it moves toroidally from the original field lineby a distance equal to x = v prec τ orb , where v prec = γmv ˆ seBR , (5)is the toroidal precession velocity, R the major radius,ˆ s the magnetic shear, and τ orb = 2 πqR v , (6)is the transit time. The particle returns to the zonewhere it still remains correlated with the original perturbation if (cid:15) ≡ v prec τ orb λ ⊥ < . (7)The parameter (cid:15) is referred as the orbit-averagingvalidity parameter .When the particle returns to (or does not leaveat all) the zone of correlation, (cid:15) <
1, the transportdecreases because the particle experiences effectivelysmaller perturbation as the fluctuating magneticperturbation is averaged along the particle trajectory.This effect becomes more prominent as the orbit widthbecomes comparable to the perpendicular correlationlength, eventually leading to a Υ ∝ γ − scaling intransport when d orb (cid:38) λ ⊥ . The explicit formulae forΥ are given in Ref. [10]Υ ≈ (cid:34) − (cid:18) d orb λ ⊥ (cid:19) + 52 (cid:18) d orb λ ⊥ (cid:19) (cid:35) ∀ d orb (cid:46) λ ⊥ , (8)Υ = λ ⊥ d orb ∀ d orb (cid:38) λ ⊥ . (9)If the electron energy is high enough, the orbit-averaging also happens on the gyromotion scale if ρ g (cid:38) . λ ⊥ . This additional scaling is obtained fromEqs. (8) – (9) when d orb is replaced with ρ g . Since theperturbation is averaged both along the poloidal andthe gyro orbit, eventually the transport decreases asΥ ∝ γ − .If the particle does not return to the zone ofcorrelation, (cid:15) >
1, the orbit-averaging mechanism isnot valid. However, the particle becomes decorrelatedduring the radial excursion if the orbit width is wideenough, d orb > λ ⊥ , and the characteristic time-scalefor this perpendicular decorrelation τ ⊥ ≡ τ orb λ ⊥ πd orb , (10)is smaller than that of the parallel decorrelation τ (cid:107) ≡ λ (cid:107) v (cid:107) , (11)i.e. τ ⊥ < τ (cid:107) . The energy scaling for this process isgiven in Ref. [10]Υ = λ ⊥ λ (cid:107) d orb ∀ d orb > λ ⊥ . (12)The case d orb < λ ⊥ is not covered in the literature,but we expect that the transport to remain constantin that case or it might increase to its original level ifit was reduced before orbit-averaging becomes invalid(the dashed line in Fig. 1). The orbit-averaging is ssessing energy dependence of the transport of relativistic electrons ρ g (cid:38) . λ ⊥ and,thus, one again observes Υ ∝ γ − at sufficiently highenergies.These scaling laws are based on the small Kubonumber approximation,
K ≡ ( λ (cid:107) /λ ⊥ )( δB/B ) < K <
3. Energy scan in perturbed JET and ITERfields
Even though the transport is expected to decrease atsufficiently high energy, a more interesting point is toknow whether this happens within plausible range ofRE energies. The answer to this question depends onthe perpendicular correlation length λ ⊥ , which dictatesthe scaling laws and whether the orbit-averaging isvalid. The value of λ ⊥ depends on the magnetic fieldstructure and the nature of the perturbation, makingit impossible to give a universal answer. Instead,here we assess the energy dependence of the transportby investigating three cases which are relevant forRE dynamics and mitigation in disruptions, eachrepresenting a different type of magnetic configuration.The first case represents a magnetic field that hasbecome fully stochastic, as shown in Fig. 2a, duringthe thermal quench phase in a plasma disruption. Thesecond case is during the pre-thermal quench phase,when the core has not yet become stochastic and somemagnetic islands remain at the stochastic edge region(Fig. 2b). These two cases are important because itis during the thermal quench, when the initial seedpopulation for the RE avalanche is born. Duringthe current quench phase following the completestochastization, the flux surfaces begin to heal, whichcan be expected to lead to better RE confinement.The magnetic field data for the previous twocases are obtained from a JOREK 3D non-linearMHD simulation [15, 16]. The simulated case is amassive injection of Argon gas in JET pulse number85943, which led to the formation of a RE currentof 1 MA. Details of the model will be described ina forthcoming publication [17]. This simulation hasbeen chosen in particular because a plasma current ( I p )spike of comparable magnitude to the experimental oneis obtained, which we interpret as a likely sign thatthe level of magnetic stochasticity is well reproduced.Indeed, theoretical work [18] as well as JOREKsimulations indicate a correlation between the I p spikeheight and the level of magnetic stochasticity.The third case (Fig. 2c) corresponds to a proposedscenario for RE mitigation, where the field is made stochastic using external coils. In ITER, the ELMcontrol coils could be used for this purpose, althoughthey are believed to be insufficient for RE mitigationon their own [19, 20]. Nevertheless, we use thecurrent flat-top magnetic field (no post-disruption datais available) of the ITER baseline scenario perturbedwith ELM control coils for the transport study. Thefield is similar and constructed with same means asto what was used in an earlier study [21] that didnot investigate FOW effects. The coils are set to an N = 3 mode, where N is the toroidal periodicity, withcoil current I = 45 kA which is half the maximumvalue (this is the setup foreseen to be used for ELMmitigation [22]). The field is constructed by addingthe perturbation due to the coils, calculated with aBiot-Savart solver [23], on the equilibrium field. Thisresults in a field where a stochastic layer is generatedat the edge but leaves the core intact.In addition to these three realistic cases, wehave constructed an artificial field by imposingarbitrary resonant magnetic perturbations to the ITERequilibrium. These perturbations create a stochasticlayer at the edge, with no islands, as shown in Fig. 2d,for which the quantities of interest, namely δB/B , λ (cid:107) , and λ ⊥ , are adjustable. The perturbations arestationary and have the form δ B = δ ∇ × ( α B D ) (13)where B D is the (unperturbed) axisymmetric field, δ a perturbation scaling factor, and α ( ρ, θ, ζ ) = (cid:88) n,m α nm ( ρ ) cos( nζ − mθ − φ ,nm ) . (14)Here ( ρ, θ, ζ ) are the radial, poloidal, and toroidalstraight-field-line coordinates, respectively. Each modeis defined by its radial eigenfunction α nm , poloidal, m , and toroidal, n , mode numbers, and phase φ ,nm ,which is chosen to be random.The eigenfunctions are Gaussian functions thatpeaks at the corresponding resonant surface, α nm ( ρ ) =exp( − ( ρ − ρ nm ) / σ ) where q ( ρ nm ) = m/n , andthe width σ is left as a free parameter. A total of25 modes were used to create the stochastic layerat the edge, and the mode numbers n and m werechosen so that the resonant surfaces were distributedat approximately equal intervals in ρ (see Fig. 3). Thescaling factor δ is chosen so that δB/B ≈ × − .Note that Ref. [12], where the orbit-averaging processwas originally presented, used similar perturbations tostudy the transport in perturbed magnetic fields. It was found in Refs. [21, 24], that modellingthe transport as a purely diffusive process can be ssessing energy dependence of the transport of relativistic electrons a) JET fully stoc. b) JET edge stoc. c) ITER coil d) ITER RMP Figure 2.
Magnetic field Poincar´e-plot at φ = 0 ◦ Rz -plane for the investigated cases. a) Fully stochastic JET magnetic field duringthe thermal quench. b) Partially stochastic JET magnetic field during the pre-thermal quench phase. c) ITER current flat-topequilibrium perturbed with ELM control coils. d) ITER current flat-top equilibrium perturbed with artificial resonant magneticperturbations. The red and black orbits illustrate the trajectories of 100 MeV and 10 keV electrons which are the extreme valuesused in the simulations. On each plot, the trajectories begin at the outer mid-plane at the same radial position where the markerswhere initialized in the actual simulations.
Normalized minor radius01 M od e a m p lit ud e ( a r b . un it s ) Figure 3.
Radial distribution of the mode eigenfunctionsin ITER RMP simulations. Each mode has an eigenfunctionwhich is a Gaussian whose peak is located on the mode’s rationalsurface. The modes (dashed lines) are distributed radially almostevenly. The solid line shows the total perturbation. inaccurate. Better correspondence between the actualand the modelled transport is obtained by alsoincluding an advection component to the transportmodel [21] (that is in addition to the effectiveadvection caused by the gradient in the diffusioncoefficient). However, we emphasize that the choiceof modelling the transport as an advection-diffusionprocess was motivated by its numerical practicality forimplementation in reduced models and the conveniencein evaluating the transport coefficients.The transport coefficients are evaluated with theMonte Carlo method from several markers, represent-ing electrons, which are traced in a perturbed mag-netic field using the orbit-following code ASCOT5 [25]. The transport coefficients are evaluated by simulatinga population of markers that were initially located onthe same radial position in the stochastic region. Theinitial position in each case is displayed in Fig. 2. Themarkers were simulated until the time distribution ofthe cumulative losses saturated, and the distributionwas subsequently used to extract the transport coeffi-cients. In effect, the coefficients obtained via this pro-cess do not describe transport locally, but they canbe considered to describe the average transport withinthe entire stochastic region. The details of obtainingtransport coefficients from the orbit-following simula-tion results are explained in Appendix A.In these simulations, all transport will be dueto magnetic field perturbations, as the electrons areassumed to be collisionless and the radiation reactionforce is ignored. In addition, also the accelerationdue to electric field is omitted so the electron energyremains constant, allowing us to scan the transport fordifferent values of electron energy.Since we expect to see a reduction in transportnot only due to finite poloidal orbit width effects butalso due to finite Larmor radius (FLR), each setting issimulated twice: once solving for the full gyro-motion,and another time solving for the guiding center motion.The FLR effects can be isolated by comparing theresults.The simulations are done for different energiesranging from 10 keV to 100 MeV. The purpose isto observe whether the expected energy dependence(recall Fig. 1) can be recovered in magnetic fieldsrelevant for RE mitigation. Furthermore, we carry outadditional simulations where the field lines themselvesare traced instead of electrons, in order to obtain thezero-orbit-width results for comparison. ssessing energy dependence of the transport of relativistic electrons J ET f u ll s t o c . Field line transport
Advection [m/s]
FLR/GC : v / v = 0.99FLR/GC : v / v = 0.90
01 1e5 Diffusion [m /s]01 J ET e dg e s t o c . k e V k e V M e V M e V M e V Electron energy01 I TE R c o il k e V k e V M e V M e V M e V
03 1e2
Figure 4.
Transport coefficients evaluated with orbit-following simulations. Each row corresponds to a different case. The thinbright lines correspond to gyro-orbit simulations, i.e., the ones including finite Larmor radius effects. The thick faded lines correspondto the guiding-center results. Blue color corresponds to particles with pitch v (cid:107) /v = 0 .
99 and red color to v (cid:107) /v = 0 .
90. Each lineis drawn from a simulation with 32 distinct energies, and the Monte Carlo noise was reduced by taking a moving average beforeplotting. The dashed black line is the magnitude of the field-line transport coefficient, i.e., transport of massless particles travellingat the speed of light along the magnetic field.
The transport coefficients evaluated using the orbit-following simulations are shown in Fig. 4 for the threecases with realistic magnetic field. The transport inthe ITER field with artificial perturbations will bediscussed separately in the next section. The resultsshow good agreement between the guiding-center andgyro-orbit simulations, indicating that the FLR effectsplay only a minor role in the studied cases. Thetransport does depend on pitch as the results with ξ ≡ v (cid:107) /v = 0 .
90 are generally lower than the ones with ξ = 0 .
99. However, the difference is in line with thetrivial pitch dependence, v (cid:107) ∝ ξ which is present evenin the Rechester-Rosenbluth diffusion coefficient givenin Eq. (1). The energy dependence of the transportshows the expected behavior for low energies ( E < v (cid:107) . When v ≈ c , the transport either peaks or flattens at thelevel of the field-line transport. For higher energies,the behavior differs in each case. The first row:
The fully stochastic magnetic fieldcase has both advection and diffusion rising above the field-line transport level at
E >
10 MeV. This resultis inconsistent with with the hypothesis that Υ ≤ > The second row:
The edge stochasticity case hasa drop in both advection and diffusion above a fewMeV. There is a simultaneous decrease in the numberof markers that are lost within the simulation time.For energies above 20 MeV, the few markers whichare lost, do so in groups: the transport has becomecompletely laminar and the coefficients cannot besensibly evaluated, hence the abruptly ending lines.These effects are caused by the magnetic field topology:the electrons with high energy cross the region wherethe flux surfaces are intact, and no stochastic fieldline transport can occur. For example, 100 MeVelectrons spend a major part of their orbit within the ssessing energy dependence of the transport of relativistic electrons
The third row:
The ITER coil case has markersinitialized at the very edge where there are no islandsand even the E = 100 MeV electrons cannot reachthe intact flux surface region initially. Hence, thereis no similar cutoff in transport at higher energies asobserved in the previous case. However, the transportdecreases with increasing energy, which could beattributed to FOW effects. On the other hand, theslope indicates a reduction weaker than the predictedΥ ∝ γ − .To summarize, none of the considered JETcases exhibit a reduction of transport with energycompatible with orbit averaging or orbit decorrelationeffects. In the fully stochastic case the transport increases as a function of energy for reasons to beclarified. In the case with stochastic edge, thetransport decreases but this is due to particles spendingmore time in their excursion in the healed flux surfaceregion. Only the ITER coil case shows a reductionin transport with increased energy that could beattributable to FOW effects. However, the reduction ofthe transport at high energy is weaker than expected.
4. Comparison to analytical results
In order to better understand the results in the realisticcases, we now study the transport in the presenceof artificial magnetic perturbations superimposed onthe ITER equilibrium, allowing us to adjust theperturbation parameters at will.Based on the discussion in section 2, FOW effectsbecome apparent when the perpendicular correlationlength, λ ⊥ , is comparable to the electron orbit width.The orbit width for a 100 MeV electron in ITER isapproximately 0.1 m. To observe FOW effects over theenergy range of 10 keV – 100 MeV, we therefore seekto have λ ⊥ ≈ d orb ∝ γ when γ (cid:29) λ ⊥ ∼ σ , when the perturbations have Gaussian radialprofiles. For the parallel correlation length, we mayuse the estimate λ (cid:107) ≈ πqR/n .We choose the toroidal and poloidal modenumbers as n (cid:46) and m (cid:46)
20, so that the mode resonant surfaces are evenly spaced radially ‡ . In total there are25 adjacent modes, and the smallest possible widththat still allows the modes to overlap and keeps theedge region stochastic is σ = 0 .
01 m. We repeat thesimulations for three cases where the mode width isvaried from this smallest possible value to the width ofthe entire stochastic layer.With these choices, we estimate that λ (cid:107) ≈ λ ⊥ ≈ .
01 – 0.2 m. Appendix B presents howthe correlation lengths are computed numerically. Inthe following discussion, we rely on the numericallyevaluated values.The magnetic field parameters for each simulatedcase are listed in Table 1. The first and second columnsshow the parallel and perpendicular correlation length,respectively, and the third column is the flux surfaceaverage of the perturbation magnitude on the initialradial position. The fourth column is the Kubo numberevaluated from the aforementioned parameters. In allcases
K (cid:28) D num , and the analytical estimate,Eq. (1). This ratio is tabulated in the fifth column,and it is close to unity in almost all cases, with theexception of the stochastic edge JET case, possiblybecause there are major magnetic islands present.Table 1 also lists threshold energies , correspond-ing to the following situations: the orbit-averaging be-comes invalid, the orbit width exceeds the perpendic-ular correlation length, and when the orbit-averagingoccurs on the gyromotion scale. The final column liststhe threshold energy for τ ⊥ < τ (cid:107) , the third conditionfor the perpendicular orbit decorrelation mechanism(the two others were (cid:15) > d orb > λ ⊥ ).We first compare these threshold energies to thecalculated transport coefficients in the artificial field(shown in Fig. 5) before returning to the discussion ofthe realistic cases studied earlier. Each case in Fig. 5exhibits an increasing transport as v → c , until itsaturates and begins to decrease. In each case Υ ≤ The first row:
Orbit averaging becomes invalidright when v ≈ c is reached and no transport reductionis observed until d orb > λ ⊥ . Afterwards there isa slight reduction (when the transport should beconstant), which could be caused by a transition inthe mechanism driving the transport, from paralleldecorrelation to perpendicular decorrelation as τ (cid:107) → τ ⊥ . The reduction in transport steepens as theperpendicular decorrelation and the gyro-orbit effects ‡ The dominant modes in the realistic cases were n ≤
3, butwith n = 3 there were not enough rational surfaces to enableseveral adjacent modes. ssessing energy dependence of the transport of relativistic electrons Table 1.
Magnetic field parameters and threshold energies.
Magnetic field parameters Threshold energies [MeV] λ (cid:107) [m] λ ⊥ [m] δB/B K D/D num (cid:15) > d orb > λ ⊥ ρ g > . λ ⊥ τ ⊥ < τ (cid:107) n = σ = 0 .
01 m 2 0.006 9 × − σ = 0 .
03 m 7 0.04 7 × − σ = 0 . × − × − × − × − × × (ITER) = . m Field line transport
Advection [m/s] 06 1e2 Diffusion [m /s]04 = . m k e V k e V M e V M e V M e V Electron energy02 = . m FLR/GC : v / v = 0.99FLR/GC : v / v = 0.90 k e V k e V M e V M e V M e V
05 1e2
Figure 5.
Transport coefficients evaluated with the orbit-following simulations for the ITER RMP case with n = 10 and fordifferent values of mode width σ . The various grey regions in the background indicate when the orbit averaging has become invalid, (cid:15) >
1, (light grey with ” \ ”-stripes), and when τ ⊥ < τ (cid:107) (dark grey with ”/”-stripes). The plain white background is where theorbit-averaging is valid. The black vertical lines indicate when the threshold energies are reached: thin lines indicate the d orb > λ ⊥ threshold whereas the thick lines are for the ρ g > . λ ⊥ threshold (for the v (cid:107) /v = 0 .
90 case). become relevant—by a coincidence, almost at the sameenergy. Here we observe the γ − relation betweenthe transport and the electron energy. However, eventhough the FLR effects should be present, only a smalldifference is observed when those are accounted for. The second row:
The orbit-averaging is valid for higher energies than in the previous case. Equation (8)predicts a maximum of 10% reduction until the orbit-averaging becomes invalid. However, no reduction isobserved until the point where the perpendicular orbit-decorrelation becomes valid. At the very end of theenergy range the FLR effects become valid, and there ssessing energy dependence of the transport of relativistic electrons
The third row:
In this case, the orbit averaging isvalid for even higher energies, but the notable featureis the valley near E = 50 MeV. If the orbit-averagingwould be valid at that point, the reduction in transportwould be 15%, which is close to the observed value. Wesuggest that the increase beyond this point is due tothe orbit-averaging transitionally becoming invalid andthe transport returning to Υ = 1. In other words, whatis seen here corresponds to the dashed line in Fig. 1.At the end of the energy range, the transport peaksand begins to decrease again. This is at the thresholdwhen the perpendicular decorrelation becomes valid.In summary, we observe a good, though not anexact, agreement between the theoretical estimates andthe simulations. The differences are mostly related tothe orbit-averaging process.
5. Realistic cases revisited
Having verified the theory through the simulations,and vice versa, we can now compare the transport inthe realistic cases to the theoretical predictions.In the JET cases, recall Fig. 4, the orbit averagingis valid until E = 50 MeV. At this energy, where theelectron orbit width is approximately 10 cm, Eq. (8)predicts a reduction of 35% in transport. However,this reduction is not visible in the results due tothe transport barrier effect being dominant in thefully stochastic case and the more energetic electronsbecoming confined in the case with the stochastic edge.The perpendicular decorrelation mechanism becomesactive at 100 MeV, when d orb > λ ⊥ , and one canobserve a drop in transport in the first row in Fig. 4that could be attributed to this.We can confirm the presence of the transportbarrier at the edge in the fully stochastic JET caseby moving inwards the loss surface, i.e., the boundarybeyond which markers are considered lost. Figure 6shows the energy dependence of the transport whenthe loss surface has been moved over the barrier. Nowthe reduction in transport is similar to what the theorypredicts.The transport barrier appears because δB/B isnot uniform radially, see Fig. 7a, and the perturbationis significantly weaker at the edge. However, thisbarrier is weaker for more energetic electrons becausetheir orbit drift takes them closer to the core where theperturbation is stronger. Therefore, the orbit-averagedvalue of δB/B increases with energy for electrons at theedge. This increases transport the effect being largerthan the reduction due to FOW effects. From Fig. 7b,we can observe that δB/B in the fully stochastic casepeaks where the dominant (1 ,
2) mode has a peak. The k e V k e V M e V M e V M e V Electron energy01 J ET f u ll s t o c . Field line transport
Diffusion [m /s] FLR/GC : v / v = 0.99FLR/GC : v / v = 0.90 Figure 6.
Diffusion coefficient in the fully stochastic JETcase when the loss surface is moved inwards over the transportbarrier.
Normalized poloidal flux ( )01 M od e a m p lit ud e ( a r b . un it s ) ( , ) ( , ) ( , ) ( , ) ( , ) P e r t u r b a ti on o r b it - a v e r a g e B / B
1e 2
Orbit widthat 100 MeV E = M e V E = k e V a)b) Figure 7. a) Mean perturbation amplitude along theparticle trajectory for 10 keV and 100 MeV electrons in thefully stochastic JET case. Horizontal axis is the marker initialposition at the outer mid-plane. Dashed vertical line denotes theposition from where the markers were launched when evaluatingthe transport coefficients. The solid vertical line is the newloss surface that has been moved inwards from the edge. Thehorizontal grey bar illustrates the orbit width of an 100 MeVelectron. b) Radial profiles of the dominant modes in the fullystochastic JET case. The brackets ( n, m ) are the toroidal andpoloidal mode numbers of the mode. modes (1 ,
3) and (1 ,
4) that have resonant surface at theedge are comparably lower which explains the drasticreduction in δB/B .As for the ITER case, the orbit averaging isvalid for the entire energy range. At 100 MeV, thepredicted decrease in transport is 10% which is asignificantly smaller reduction than what was observedin the simulations. One should note that this casediffers from all the others in that the perturbation is ssessing energy dependence of the transport of relativistic electrons δB/B is strongest near the coilsat the low-field side. The theory assumes uniformperturbation, which might explain the difference seenhere.As the final point, we note that the numericallyevaluated parallel correlation length in the artificialITER cases is close to the estimate based onthe toroidal mode number, and the perpendicularcorrelation lengths are comparable to the radial widthof the modes. In the JET cases, the dominantmodes (recall Fig. 7b) have n = 1, which gives anestimate λ (cid:107) ≈
50 m that is an order of magnitudehigher than the numerically evaluated value. On theother hand, the perpendicular correlation length iscomparable to the mode width which is roughly halfof the minor radius, i.e. 0.6 m. Therefore, one canuse the estimate λ ⊥ ∼ mode width, in order to makean initial assessment whether FOW effects should beconsidered.
6. Summary and conclusions
While the orbit-following simulations broadly agreewith the theoretical results in Ref. [10], it was foundthat in disruptions magnetic field structures andnon-uniform perturbation may dominate the energydependence of the transport. In extreme cases, thetransport of energetic particles ceases completely whentheir orbit crosses confined field line regions, or,somewhat unexpectedly, the transport increases abovethe transport level that would be observed if theparticles would follow the field lines exactly.The evaluation of the radial advection anddiffusion coefficients was motivated by the prospect ofincluding them in a reduced kinetic model, in orderto capture the effect of the 3D magnetic field on therunaway electron dynamics. Based on our results, ifthe magnetic field is spatially inhomogeneous, it ismore appropriate to evaluate the transport coefficientsnumerically, instead of relying on the analytical model.Simulation results differ from theory mostly inthe orbit-averaging regime. This was especially truein the case where the magnetic field was perturbedwith external coils, which could be attributable tothe perturbation being poloidally localised in thiscase. Simulations done within the guiding centerapproximation yielded results that were practicallyidentical to those obtained by tracing the completegyromotion. Only for
E >
50 MeV in some cases, theguiding center results showed higher transport. Wenote that Ref. [26] reported a difference between gyro-orbit and guiding center results, but this was due toparticles being born inside the islands which was notconsidered here.One of the motivations for investigating the energy dependence of the transport was the discrepancy wherethe RE transport observed in experiments was orders ofmagnitude lower than that predicted by the Rechester-Rosenbluth diffusion coefficient. However, here FOWeffects were found to reduce the transport by an orderof magnitude only at E = 100 MeV, which is onthe high end of plausible runaway electron energyspectra in ITER-like disruptions and is well beyond theenergy of seed populations in any runaway scenario.The energy reduction depends on the perpendicularcorrelation length of the magnetic perturbation, whichwas found to be of the same order as the radial widthof the mode. Importantly, when the stochastic field iscaused by MHD activity, such that the mode widths arecomparable to minor radius, and not due to microscaleturbulence, the transport of REs is not reduced in theenergy range of interest.In our simulations the Rechester-Rosenbluthestimate fared well in the zero orbit-width limit exceptfor one case, the one with major islands present,where it over-estimated the transport by two ordersof magnitude. We emphasize that this reduction wasnot because particles are born inside islands where theyremain trapped, as was the case in Ref. [26], as in thiscase the markers where initialized within the stochasticregion; instead the mere presence of the islands reducedthe transport. Acknowledgments
We acknowledge the CINECA award under the ISCRAinitiative, for the availability of high performancecomputing resources and support. The authors aregrateful to M. Hoppe, P. Svensson and I. Pusztaifor fruitful discussions. This work was supported bythe European Research Council (ERC-2014-CoG grant647121), the Swedish Research Council (Dnr. 2018-03911), and the EUROfusion - Theory and AdvancedSimulation Coordination (E-TASC). The work hasbeen carried out within the framework of theEUROfusion Consortium and has received fundingfrom the Euratom research and training programme2014-2018 and 2019-2020 under grant agreement No633053. The views and opinions expressed hereindo not necessarily reflect those of the EuropeanCommission.
Appendix A. Transport coefficient evaluation
We assume the radial transport of particles is given bythe following advection-diffusion equation ∂f∂t = − ∂∂ρ Kf + ∂ ∂ρ Df, (A.1)where f ( ρ, t ; µ, E ) is the distribution function, and K ( ρ ; µ, E ) and D ( ρ ; µ, E ) are the advection and ssessing energy dependence of the transport of relativistic electrons Time [s]01 L o ss fr ac ti on L a m i n a r Figure A1.
Illustration of the transport coefficient evaluationand loss times. The cumulative first-passage-time distribution,Eq. (A.2), is fitted (dotted black lines) to the cumulative particleloss-time distribution (solid blue lines), and the coefficients areobtained from the fit parameters. Two cases are shown: theone on the left has also laminar losses which appear as “steps”in the loss time distribution. The orange line is the loss-timedistribution where the coefficients obtained from the fit are usedin a transport model, but employing a reflective inner boundary. diffusion coefficients, respectively. The energy, E , andmagnetic moment, µ , are kept as parameters since theyremain invariant in the transport process we consider.The transport coefficients are likely to have aradial dependence because the magnitude of δB/B may vary and any islands that are present affectthe transport. However, the radial dependence isneglected here and, instead, the transport coefficientsare evaluated only at a single radial position. Atthis position, we initialize a population of markersrepresenting electrons. The markers are initialized atthe outer mid-plane and their toroidal coordinate issampled from a uniform distribution.The markers are traced with ASCOT5 and thetransport coefficients are extracted from the results.There are several ways this can be done [21], and herewe choose to evaluate the coefficients based on thedistribution of loss time, i.e. the time it takes for amarker from the beginning of the simulation to crossthe separatrix for the first time. Some markers mightbe radially confined, e.g. inside islands, and, therefore,the marker population is simulated until this loss-timedistribution becomes saturated.The transport coefficients are extracted from theloss-time distribution by fitting an analytical modelto the results; however, the analytical model isbased on several simplifications. When one assumesthat the transport is uniform radially (i.e. K and D do not depend on position), the advection termis positive (i.e. the flow is towards the edge), andthe inner boundary extends to minus infinity § ,then the cumulative loss-time distribution obeys the § In reality, the inner boundary is reflective and it is located atthe magnetic axis or at the boundary between the confined andstochastic field-line regimes. so-called first-passage time distribution [27] which,in probability theory is known as the inverseGaussian distribution. The corresponding cumulativedistribution is T ( t ) =Φ (cid:18) √ Dt ( Kt − ∆ ρ ) (cid:19) + exp (cid:18) K ∆ ρD (cid:19) Φ (cid:18) − √ Dt ( Kt + ∆ ρ ) (cid:19) , (A.2)where Φ is the cumulative normal distribution. Thefirst-passage time refers to the time instant when arandom walker launched from a given position firstpasses through a pre-defined coordinate (separatedby a distance ∆ ρ ), which exactly is the case herewhen the simulations measure the time when a markercrosses the separatrix for the first time. In additionto fitting, the moments of the loss-time distributioncan be used to estimate the transport coefficients as K = ∆ ρ/ Mean[ t ] and D = (∆ ρ ) Var [ t ] / t ] ,although the results are more susceptible to noise.Even though the approximations behind Eq. (A.2)are drastic, the formula fits the data well: at any timethe difference in losses deviate less than 10 % fromthe simulation results. Two example fits are shown inFig. A1: the one on the left (faster losses) is for 100MeV electrons ( ξ = 0 .
9) and the one on the right is for10 keV in the fully stochastic JET case.The data on the left has a step-like structure atthe beginning indicating that multiple markers arelost at the same time. In this laminar transport particles are lost before they have decorrelated, andthe advection-diffusion model does not describe thisearly-time behaviour. Still, the trend is captured bythe fit even though the details are not.For the second case the shape of the fit is slightlydifferent than the data, which is due to the actualtransport not being uniform in ρ . We can assesswhat effect omitting the reflective inner boundary inEq. (A.2) has, by including one in Eq. (A.1) and solvingthe equation numerically. We use the coefficientsobtained from the fit (second case), an initial radialdistribution that is identical to the one used inthe orbit-following simulations, and set the reflectingboundary right next to the marker initial position. Thecumulative loss-time distribution obtained this way(the orange line) deviates from the earlier results asparticles are lost faster due to the reflection. The take-away message here is that since a reflective boundaryis inherently present in the data to which Eq. (A.2)is fitted, we can expect that the coefficients we obtainoverestimate the actual transport.Finally, the advection is negative only in thepresence of islands but islands are not present in mostcases studied in this work. ssessing energy dependence of the transport of relativistic electrons Appendix B. Correlation length evaluation
The (auto)correlation length is defined as λ ≡ (cid:90) ∞ C ( l ) dl, (B.1)where C ( l ) is the (auto)correlation function of themagnetic field perturbation: C ( l ) ≡ (cid:68)(cid:82) ∞−∞ ˜ B ( r ) ˜ B ( r + l ) dr (cid:69)(cid:68)(cid:82) ∞−∞ ˜ B ( r ) dr (cid:69) , (B.2)where ˜ B is the perturbation component that isperpendicular to the unperturbed field. The bracketsare averages over all realizations of ˜ B .In most cases the perturbation is not isotropic;in this paper, for example, the perturbation potentialoscillates along the field lines. Therefore, it iscustomary to consider the parallel and perpendicularcorrelation lengths and functions separately. For theparallel correlation length (function), the integrals inEqs. (B.1) – (B.2), are performed along the field lines.The correlation function peaks at l = 0 due tothe convolution in the numerator (the denominator isjust a normalization factor). It is generally assumedthat, in both directions, C ( l ) has the form of Euleriancorrelation function: C ( l ) ∼ exp (cid:32) − l λ (cid:107) (cid:33) , (B.3)but this is valid only if ˜ B is a Gaussian random noise.We evaluate the parallel correlation functionnumerically by choosing a radial position, and thenfinding the field-line ( R, φ, z ) coordinates in theaxisymmetric field. The integrals in Eq. (B.2)are evaluated in the perturbed field along thesecoordinates, after which we move the toroidalcoordinate by a random angle, and repeat theevaluation. The process is repeated until theensemble average has converged and we obtain thecorrelation function. Figure B1 illustrates how thenumerically evaluated correlation function compares tothe estimate in Eq. (B.3).
References [1] Lehnen M, Aleynikova K, Aleynikov P, Campbell D,Drewelow P, Eidietis N, Gasparyan Y, Granetz R, GribovY, Hartmann N, Hollmann E, Izzo V, Jachmich S, KimS H, Koˇcan M, Koslowski H, Kovalenko D, Kruezi U,Loarte A, Maruyama S, Matthews G, Parks P, PautassoG, Pitts R, Reux C, Riccardo V, Roccella R, Snipes J,Thornton A and de Vries P 2015 J. Nucl. Mater.
39– 48 doi: 10.1016/j.jnucmat.2014.10.075
Distance along the orbit [m]101 N o r m a li ze d c o rr e l a ti on f un c ti on Figure B1.
Numerically evaluated magnetic field (parallel)correlation function. Horizontal axis is the distance in meters.The blue lines are the correlation functions calculated forindividual field lines, and the red line is the ensemble averageof these, i.e. the actual correlation function in Eq. (B.2). Theblack line is the Eulerian correlation function, Eq. (B.3), where λ (cid:107) was obtained by integrating the actual correlation function(the red line) according to Eq. (B.1).[2] ITER Physics Expert Group on Disruptions, PlasmaControl, and MHD and ITER Physics Basis Editors1999 Nucl. Fusion (25)255003 doi: 10.1103/PhysRevLett.100.255003[6] Lin Z F, Chen Z Y, Huang D W, Huang J, Tong R, WeiY N, Yan W, Li D, Hu Q M, Huang Y, Yang H Y,Li Y, Zhang X Q, Rao B, Yang Z J, Gao L, DingY H, Wang Z J, Zhang M, Liang Y, Pan Y and andZ H J 2019 Plasma Phys. Control. Fusion https://doi.org/10.1088/1361-6587/aaf691 [7] Riccardo V, Arnoux G, Cahyna P, Hender T C, Huber A,Jachmich S, Kiptily V, Koslowski R, Krlin L, Lehnen M,Loarte A, Nardon E, Paprok R and and D T 2010 PlasmaPhys. Control. Fusion ssessing energy dependence of the transport of relativistic electrons [15] Huysmans G and Czarny O 2007 Nucl. Fusion TBD
To be submitted to Phys.Plasmas[18] Boozer A 2019 Plasma Phys. Control. Fusion Preprint )URL https://arxiv.org/abs/1511.01629 [24] Papp G, Drevlak M, Pokol G I and F¨ul¨op T 2015 J. PlasmaPhys. doi: 10.1017/s0022377815000537[25] Varje J, S¨arkim¨aki K, Kontula J, Ollus P, Kurki-Suonio T,Snicker A, Hirvijoki E and ¨Ak¨aslompolo S 2019 High-performance orbit-following code ASCOT5 for MonteCarlo simulations in fusion plasmas Submitted to Comp.Phys. Comm. ( Preprint ) URL https://arxiv.org/abs/1908.02482 [26] Carbajal L, del Castillo-Negrete D and Martinell J J 2020Phys. Plasmas33