Electron heat diffusivity in radially-bounded ergodic region of toroidal plasma

Drift-kinetic δf simulations are performed to investigate effect of ergodic field lines caused by resonant magnetic perturbations (RMPs) on radial heat diffusivity of electrons in the edge region of toroidal plasma of collisionality ν∗∼0.1. The following is assumed in the simulations. The ergodic region is bounded radially on both sides by closed magnetic surfaces. The pressure gradient remains nonzero in the ergodic region because of an incomplete flattening of the pressure profile, and the characteristic scale length of the pressure gradient is of larger order than the overlapping width of the neighbouring magnetic islands. It is found in the quasi-steady state of δf that the electron heat diffusivity is of smaller order than the theoretical estimate derived by the Rechester–Rosenbluth model (Rechester and Rosenbluth 1978 Phys. Rev. Lett. 40 38). The radial heat conduction is dominated not only by parallel motions along the ergodic field lines, but also by trapped particle motions generated by the RMP field. The contribution of the trapped particles reduces the radial heat conduction enhanced by the parallel motions.


Introduction
In recent tokamak experiments, resonant magnetic perturba tions (RMPs) have been used to suppress/mitigate edge local ized modes (ELMs) in the plasma edge region of collisionality ν * ∼ 0.1 [1,2]. The definition of ν * is given in section 2. In these ELM suppression experiments [3,4], it was found that the experimental estimate of the radial heat diffusivity of electrons affected by the RMPs is fairly small compared to the theor etical estimate derived by the Rechester-Rosenbluth (hereafter, R-R) model [5,6]. On the other hand, if magn etic field lines affected by the RMPs become ergodic (or stochastic), and if the width of the ergodic region is approx imately the same as the width of the edge region, the R-R model is considered to be applicable for estimating the heat diffusivity in the edge region. The applicability of the R-R model in the estimation of the diffusion coefficient has been verified by a Monte Carlo simulation of guiding centre motions [7]. Therefore, there is an open question: why does there occur a difference between the theoretical results and the experimental results?
This question is composed of two parts. The first issue to be discussed is whether the RMP field penetrates the plasma in the ELM suppression experiments. This is not yet com pletely resolved at present. If the RMP field is screened by the plasma rotation [8], then field lines in the perturbed region do not become ergodic. As a result, the fieldline diffusion given by the R-R model is not realised. However, the RMP penetra tion in the case of ELM mitigation is reported in magneto hydrodynamic (MHD) simulation studies [8][9][10]. We assume in this paper that an ergodic region exists in the plasma edge region affected by the RMP field, and is bounded radially on both sides by closed magnetic surfaces, and that good confine ment in the plasma edge region is maintained, in which the ELMs are suppressed. Note that in general, the heat confine ment deteriorates if the ergodic region is not bounded radially, and field lines in the ergodic region connect to the divertor target plates [11,12]. The second issue is whether the field line diffusion given by the R-R model is realised in such a radiallybounded ergodic region. The phenomenon in which ergodic field lines are confined in a finite domain plays an important role for the magnetic confinement of the toroidal plasma. That phenomenon is discussed, showing the recent results in detail, in [13]. For example, it is well known that guiding centre motions in the radiallybounded ergodic region differ from random motions given in the R-R model [14]. Therefore, if ergodic field lines are confined in the radially bounded region, and if the pressure gradient remains nonzero in the ergodic region because of an incomplete flattening of the pressure profile, there is a possibility that the radial heat diffusivity is different from the estimate of the R-R model.
Other possibilities may explain the difference between the theoretical results and the experimental results. Recently, the kinetic ion-electron-neutral guiding centre code XGC0, which is a kinetic simulation code for treating the socalled full distribution function (i.e. a fullf simulation code), has been applied to plasma transport in an ergodic region [15], where the screening of RMPs is not investigated. It is reported that the electron heat transport is significantly reduced from the R-R model. In this simulation, many physical effects are considered. These include selfconsistent radial electric field, Coulomb/neutral collisions, particle sourcing from wallrecy cling and neutral transport, electron cooling from impurity radiation, and others. Thus, it is not clear whether the combi nation of many physical effects is essential for the explanation of the difference, or the difference is obtained in a simpler situation.
While driftkinetic δf simulation [16,17] is known to be a powerful tool for estimation of neoclassical transport in toroidal plasmas, including threedimensional field geometry, the simulation code has rarely been used in studies of col lisional transport in the ergodic region. When applying the driftkinetic δf simulation code to the transport in the ergodic region, care should be taken that the presupposition |δf /f (0) | 1 is required to be satisfied because the quadratic collision term C(δf , δf ) is ignored in the simulation code, where f (0) is the background distribution function given in advance in this paper. Here, the roles of f (0) and δf are explained in section 2. In the previous work [18], we developed a driftkinetic δf simulation code that includes the effect of the radial electric field in order to estimate the electron heat diffusivity in the ergodic region. However, the approach to the open question was insufficient because it was not clear which assumption is necessary for In this paper, this point will be discussed. Then, in order to answer the second question, which is the key point in the open question raised in this paper, we will examine the effect of ergodic field lines on the radial heat diffusivity of electrons in the radiallybounded region using the driftkinetic δf simulation under an assumption of zero electric field, which is the same assumption as provided in the R-R model. For the sake of a simpler comparison between the simulation and the R-R model, it is also assumed that effects of MHD activities, turbulence, mean flow, impurities and neu trals are ignored. Because many physical effects are neglected, it is difficult to directly compare our simulation and the fullf simulation. However, our simulation will show that the elec tron heat diffusivity is significantly reduced, when compared to the R-R model, in the simpler situation than in the fullf simulation.
This paper is organised as follows. In section 2, the ergodic region used in the simulations is briefly explained.
The assumption for f (0) satisfying |δf /f (0) | 1 is discussed. The effect of ergodic field lines on the electron heat diffu sivity is examined. The summary and discussions are given in section 3.

Simulation model
In the δf Monte Carlo simulation code KEATS [18,19], we introduce the following ordering into the representa tion of a distribution function f in the driftkinetic equation: Here, the zerothorder distribution function f (0) is set to a Maxwellian f (0) = f M = n(m/2πT) 3/2 exp{−mv 2 /2T} with fixed pro files of temperature T and number density n, where m is the par ticle mass, v is the velocity, and the mean velocity is ignored. The firstorder distribution function f (1) is the difference from f M , i.e. f (1) = δf = f − f M . The temperature and the number density are assumed to be given functions of a label of refer ence surfaces, where the reference surfaces are defined by the unperturbed magnetic flux surfaces. Due to the presupposi tion | f (1) /f (0) | 1, the driftkinetic δf simulation is required to ensure that orders of characteristic scale lengths L (0) and L (1) of the zerothorder and firstorder distribution functions are well separated from each other, i.e. L (1) /L (0) 1. In this paper, the characteristic scale lengths L (0) and L (1) where T e is the electron temperature and W island is the overlap ping width of neighbouring magnetic islands. Note that the zerothorder distribution function f M is characterised by the scale length of the temperature gradient because the gradient of the number density is set to be zero in the estimation of the electron heat diffusivity. On the other hand, the firstorder distribution function δf is strongly affected by the scale length of the ergodic region, i.e. W island , as discussed below. Here, the ratio of the electron banana width δr e b to the overlapping width W island is set to satisfy δr e b /W island 1, where δr e b is estimated as δr e b ∼ ρ e p √ t , ρ e p is the electron poloidal gyroradius, and −1 t is the aspect ratio.
Applying the decomposition f = f M + δf to the drift kinetic equation in 5dimensional phase space (x, v , v), we have the following equation with respect to electrons confined in static magnetic field B: [20,21]. Here, the electric field E is ignored, and is set to E = 0. The parallel and drift velocities are given as elementary charge, and the electron mass is written as m. The energy is given as E = (1/2)mv 2 = constant in the case of E = 0. The Coulomb collision operator is described by C T and C F , where C T represents the electron-electron and electronion pitchangle scattering operators [21], and C F is the fieldparticle collision operator which ensures the local momentum conservation property for the electron-electron collision, given in section 3.8 of [21]. The boundary condition of the simulation is that the test particles (i.e. the sample guidingcentres having the two weights in the Monte Carlo method [18]) are absorbed when arriving at the wall of the calculation box. In the elec tron transport simulations, the ion distribution function-with the species set to be proton-is assumed to be Maxwellian, with the same profiles of temper ature and number density as the electron profiles because the condition of | f (1) /f (0) | 1 is also satisfied in the case of ions [19,22]. Here, when seen by the electrons, the Maxwellian of the ions is approximately a delta function without the case that the ion temperature is significantly larger than the electron temperature. Thus, the ion temperature is assumed to be equal to the electron temperature for the purpose of simplifying the simulation.
The magnetic configuration is set to a simple tokamak field affected by RMPs. The major radius of the magnetic axis is R ax = 3 m, the minor radius of the toroidal plasma is a = 1 m, and the strength of the magnetic field on the magn etic axis is B ax = 2.5 T. The temperature profile is given as T e = T i = T(r) = T ax − (T ax − T edge ) r/a and the number density n e = n i = n is assumed to be constant (i.e. ∇ ln n = 0), where r = r(R, Z) = (R − R ax ) 2 + Z 2 is the label of the reference surfaces in the cylindrical coordinates (R, ϕ, Z). The reference surfaces are given by magnetic flux surfaces of a circular tokamak field B 0 . The unperturbed magnetic field [23], where q is the safety factor given as q = 1.8 + 3.1(r/a) 4 , and R, ϕ, and Z are the unit vectors in the R, ϕ, and Z direc tions, respectively. The RMP field, which causes resonance with rational surfaces satisfying 3 < q = k/ < 4 and = 2, 3, 4, 5, 6, 7, is given by a perturbation field δB = ∇ × (αB 0 ), which is the same type of RMP used in [7,14]. Here, the function α is used to represent the structure of the RMP field; α(R, ϕ, Z) = k, a k (r) sin{kθ − ϕ}, where tan θ = tan θ(R, Z) = Z/(R − R ax ), a k = A RMP exp{−r− r k ) 2 /∆r 2 } with A RMP /a = 2 × 10 −3 , r = r k is the rational surface of q = k/ , and ∆r is a small parameter controlling the width of the perturbation, which is set to ∆r/a = 5 × 10 −2 .
The total magnetic field B is B = B 0 + δB. The Poincaré plots of magnetic field lines affected by the RMP field and the reference surfaces in the calculation box are shown in figure 1. The Chirikov parameter, which is the indicator of the overlapping of two magnetic islands [13], is estimated as σ Chir = (1/2){w (k/ ) + w (k / ) }/∆r (k/ , k / ) 1 around the centre of the ergodic region. Here, ∆r (k/ , k / ) is the width of the neighbouring rational surfaces of q = k/ , k / , and w (k/ ) is the width of the (k, ) component of the RMPs, where W island /a w (k/ ) /a ∼ 0.02 because several magnetic islands overlap each other. Note that all of the parameters mentioned above are in an artificial device, but are set for the sake of simplicity in comparing the simulation with the R-R model.
The radial heat diffusivity of electrons, χ e r , is effectively estimated by under the conditions of n e = constant, zero electric field, and zero mean flow, where Q e r (r) is the radial energy flux of electrons in a quasisteady state of δf . In the simulations, the radial energy flux crossing a reference surface labeled by r is given as Here, · means the timeaverage, and the averaging time is of the same order as the typical time scale of δf , i.e. the electronelectron and electron-ion collision times, which are τ ee and τ ei respectively. The timeaveraging is carried out after sufficient exposure to the Coulomb collision. The average · is defined as · = (1/δV) δV · d 3 x, where δV is a small volume, and lies between two neighbouring reference surfaces with vol umes V(r) and V(r) + δV . The effective heat diffusivity in equation (2) is interpreted as the radial energy diffusivity. The statistical properties of guiding centre motions in the ergodic region bounded radially by closed magnetic surfaces have been discussed in [14]. The relation between the socalled subdiffusive behaviour of the guiding centres and the diffu sivity given from the driftkinetic equation has been studied in detail in [19]. The verification between the δf simulation and the R-R model has also been performed in [22]. Those discus sions are not repeated here. Note that the Monte Carlo simula tion of guiding centre motions in the previous works [7,14] corresponds to solving only the lefthand side of equation (1), i.e. Dδf /Dt = 0 and δf = j δ(x − x j (t)) δ(v − v j (t)) in consideration of the Maxwellian background, where j means the number of a sample guidingcentre.

Simulation results
Using the δf simulation code based on equation (1), we esti mate the effective heat diffusivity of electrons in the toroidal plasma affected by the RMPs. The effective heat diffusivity χ e r in the perturbed region in the case of T ax = 1.137 keV, T edge = 0.8T ax , and n = constant = 1 × 10 19 m −3 is shown in figure 2, where the collisionality ν * in the perturbed region is estimated as ν * ≈ 0.1. In general, the collisionality ν * is defined as ν * = νqR/ 3/2 t v th . In the simulations, the collision frequency ν is estimated to be ν ∼ ν ei ∼ τ ei −1 , the thermal speed v th is given as v th = 2T/m, and R ∼ R ax , where ν ei is the electron-ion collision frequency. The effective heat diffusivity χ e r evaluated in the simulation is fairly small com pared to χ RR r = k, πqR ax v th |δB k, | 2 /|B ax | 2 given by the R-R model [5,6], i.e. χ e r /χ RR r 10 −6 at the centre of the ergodic region r/a ≈ 0.8, where δB k, is the (k, ) comp onent of ∇r · δB. Note that in this case the perturbed region is ergodic within 0.6 < r/a < 1.0, but a closed magnetic surface remains, and divides the ergodic region, as shown in figure 1, and the simulation satisfies L (1) /L (0) 1/10.
We confirm that | d 3 v δf /n| 10 −4 and | d 3 v (1/2)mv 2 δf / (3nT/2)| 10 −3 for all r after several collision times in the simulation. Here, the number of the test particles which are absorbed at the wall of the calculation box is less than 5 × 10 −3 % of the total amount (≈2.58 × 10 8 ), where the test particles are initially distributed to be randomly uniform in the region of r/a 1 and 0 ϕ < 2π. In par ticular, the effect of the RMPs on the temperature profile is estimated by The modified temperature profile T (0+1) = T (0) + δT agrees well with the original profile T (0) = T(r) = T ax − (T ax − T edge ) r/a, as shown in figure 3(a). Therefore, the presup position of the δf method, i.e. |δf /f M | 1, is appropri ately fulfilled in this simulation. When the condition of L (1) /L (0) 1/10 is satisfied, the effective heat diffusivity χ e r given by equation (2) is independent from the temperature gradient, as shown in figure 4. It is also confirmed that the extent to which the simulation satisfies the presupposition of the δf method gradually decreases if L (1) /L (0) → 1, e.g. L (0) changes to a smaller scale, which means a steeper temper ature gradient, or L (1) changes to a larger scale, which means a broader overlapping width of magn etic islands. For example, if L (1) changes to be larger than in the case of figure 3(a), T (0+1) becomes clearly different from T (0) and flattened, as shown in figure 3(b). However, T (0) is replaced by a more gentle profile, such that the condition of L (1) /L (0) 1 is sat isfied, and T (0+1) ≈ T (0) is realised. See figure 3(c).
In the ELM suppression experiments [3,4], the col lisionality regime in the plasma edge region is estimated as ν * ≈ 0.05-1. As shown in figure 5, the result that χ e r /χ RR r 1 is not sensitive to the collisionality ν * , where ν * is set by changing the value of the constant number density. It is also shown that the electron heat diffusivity, which is effectively evaluated from the radial energy flux in equation (2), is almost independent of the collisionality ν * ≈ 0.05-1.
While the effective heat diffusivity in equation (2), i.e. the radial energy diffusivity, is independent of the collisionality ν * , the ratio of the radial heat flux q r to the radial energy flux Q r increases as the collisionality increases, ν * → 1, as shown in figure 6, where Q r = q r + (5/2)TΓ r , Γ r is the radial particle flux, and q r and Γ r are estimated in the same manner as in equation (3). This result means that the radial heat conduction estimated by q r /(n|∂T/∂r|) is weakened and, inversely, the radial particle diffusion estimated by Γ r /(n|∂ ln T/∂r|) is strengthened, when the collisionality χ χ Figure 2. Effective heat diffusivity of electrons χ e r is affected by the RMPs in 0.6 < r/a < 1.0, as shown by the solid line. The value of the diffusivity affected by the RMPs at the centre of the ergodic region r/a ≈ 0.8 is a hundred times larger than the neoclassical diffusivity illustrated by the dashed line, i.e. χ e r without RMP. Here, the temperature is set to T(r) = T ax − (T ax − T edge ) r/a with T ax = 1.137 keV and T edge = 0.8T ax . The number density is n = constant = 1 × 10 19 m −3 . The collisionality ν * is estimated as ν * ≈ 0.1 at r/a = 0.8.
is lower. Here, in all cases in figure 6, the gradient of the number density is set to be zero, and the temperature pro file is T(r) = T ax − (T ax − T edge ) r/a with T ax = 1.137 keV and T edge = 0.8T ax . Then, there is a new question: why does the radial heat conduction become weakened in the case of lower collisionality? It is naturally expected that the parallel motions of guiding centres along ergodic field lines enhance the radial heat conduction in the case of lower collision ality. In the simulation results, this tendency is confirmed, as shown by the open squares in figure 7. The contribution of the parallel motions is estimated as the contribution of untrapped particles to q r , which is evaluated as the radial heat flux satisfying the condition of κ 2 > 1, where κ 2 is the pitch angle parameter κ 2 = {(1/2)mv 2 − µB ax (1 − t )}/2µB ax t [21] and q r = q r (κ 2 > 1) + q r (κ 2 1) > 0 in all cases. The contrib utions of untrapped and trapped particles are repre sented as q r (κ 2 > 1) and q r (κ 2 1), respectively. On the other hand, the contribution of trapped particles satisfying the condition of κ 2 1 negatively enhances the radial heat conduction against the positive enhancement of the par allel motions, as shown by the closed squares in figure 7. The contribution of the trapped particles reduces the radial heat conduction enhanced by the parallel motions. Here, it is confirmed in the simulation that q r = q r (κ 2 1) > 0 and  Effective heat diffusivity of electrons around the centre of the perturbed region r/a ≈ 0.8 does not substantially depend on the temperature gradient when the condition of L (1) /L (0) 1/10 is satisfied. Here, the temperature profile is given as T(r) = T ax − (T ax − T edge ) r/a and the number density is fixed to be n = constant = 1 × 10 19 m −3 . The simulation result in the case of T ax = 1.137 keV and T edge = 0.8T ax is shown by the solid line, which is the same as in figure 1. In this case, the temperature gradient is estimated as ∆ T := |a ∂ ln T/∂r| ≈ 0.238 at r/a = 0.8. The simulation result in the case of T ax = 1.336 keV and T edge = 0.644T ax is shown by the dashed line, where the temperature gradient is 2∆ T at r/a = 0.8. Note that in both cases, the value of temperature at r/a = 0.8 is set to be approximately 0.955 keV. Figure 5. Dependence of the effective heat diffusivity of electrons χ e r on the collisionality ν * is investigated in the radiallybounded ergodic region, where χ e r is evaluated at r/a = 0.8. The diffusivity affected by the RMPs is almost independent from the collisionality ν * when ν * ≈ 0.05-1, as shown by the closed squares. In all cases, the temperature profile is set to T(r) = T ax − (T ax − T edge ) r/a with T ax = 1.137 keV and T edge = 0.8T ax , and the gradient of the number density is zero. q r /Q r ≈ 1 in the unperturbed magn etic field under the condi tions of ∇ ln n = 0 and ν * ≈ 0.1, which is the same result as in the neoclassical transport theory [21].

Summary and discussion
Using the driftkinetic δf simulation code KEATS [18], we investigate the effect of ergodic field lines caused by resonant magnetic perturbations (RMPs) on radial heat diffusivity of electrons in the cases of the collisionality ν * ≈ 0.05-1. The simulations are performed under conditions in which the ergodic region is bounded radially on both sides by closed magnetic surfaces, and the characteristic scale length of the pressure gradient L (0) is of larger order than the overlapping width of the neighbouring magnetic islands L (1) . Note that the characteristic scales of the zerothorder and the firstorder distribution functions are well separated from each other in the simulations.
It is found in the quasisteady state of δf that the elec tron heat diffusivity is of smaller order than the theoretical estimate derived by the Rechester-Rosenbluth (R-R) model [5,6]. This result indicates that the electron heat diffusivity is significantly reduced when compared to the R-R model, even in a simpler situation than the fullf simulation of [15]. The contribution of the trapped particle motions generated by the RMP field reduces the radial heat conduction enhanced by the parallel motions along ergodic field lines. The radial heat conduction is dominated not only by the parallel motions, but also by the trapped particle motions. When the collisionality is lower under the condition of ∇ ln n = 0, it is also seen that the radial heat conduction estimated by q r /(n|∂T/∂r|) becomes weakened, and that the radial particle diffusion estimated by Γ r /(n|∂ ln T/∂r|) becomes strengthened. The electron trans port properties mentioned above differ from the properties derived by the R-R model and its successor models, which are developed from the R-R model under the consideration of guiding centre motions in ergodic field lines [13], and are included in the modelling of Dδf /Dt = 0. In conclusion, the radial heat diffusion given by the R-R model is not realised in the radiallybounded ergodic region in the cases of the col lisionality ν * ≈ 0.05-1.
In general, detailed properties of radial electron and ion transport in threedimensional geometry are determined under the ambipolar condition. In addition, dependence of the elec tron transport on the detailed structure of field lines in the per turbed region is not investigated in this paper. The mechanism of the negative enhancement of the trapped particles in the radiallybounded ergodic region is not yet clear. These topics are left as future issues. Radial energy transport becomes conductive as the collisionality increases, ν * → 1, where Q r = q r + (5/2)TΓ r and Γ r is the radial particle flux. Note that in all cases, Q r is positive, and the temperature is set to be 0.955 keV. Contributions of untrapped and trapped particles to the radial heat flux, q r (κ 2 > 1)/q r and q r (κ 2 1)/q r , depend on the collisionality ν * in the radiallybounded ergodic region, as shown by open and closed squares, respectively, where the values are evaluated at r/a = 0.8 and κ 2 is defined as κ 2 = {(1/2)mv 2 − µB ax (1 − t )}/2µB ax t . Note that q r = q r (κ 2 > 1) + q r (κ 2 1) is positive in all cases. Regression curves are given by the dashed lines.