J. Plasma Fusion Res. SERIES, Vol. 9 (2010)

Particle Simulation Study of SOL Heat Transport to theDivertor Plate
Aaron FROESE1) , Tomonori TAKIZUKA2) and Masatoshi YAGI3)
1)

IGSES, Kyushu University, Kasuga 816-8580, Japan
Japan Atomic Energy Agency, Naka 311-0193, Japan
3)
RIAM, Kyushu University, Kasuga 816-8580, Japan

2)(Received: 18 November 2009 / Accepted: 25 February 2010)

Parallel heat ﬂux in a steady-state SOL is determined via fully kinetic simulations using the 1d3v
particle-in-cell code PARASOL. In order to better evaluate the heat load on the divertor plates, theion and electron velocity distributions functions (VDF), heat ﬂux, and heat ﬂux limiters αi and αeare calculated. The eﬀects of collisionality, recycling rate, ion-electron mass ratio, and source modelare examined. Results for the VDF are explained qualitatively and the average values of both ﬂux
limiters are found to be approximately 0.1, but with a large spread from 1.4 × 10−3 ∼ 2 with a complex
dependence on the simulation parameters.
Keywords: tokamak scrape-oﬀ layer, parallel heat transport, heat ﬂux limiter, velocity distributionfunction, particle-in-cell

1. Introduction

collisional plasma.

In a tokamak, hot plasma is lost from the core byanomalous transport and ELMs and carried by parallel transport along the open ﬁeld lines in the scrapeoﬀ layer (SOL) to the divertor plates. The SOL andplasma-surface interaction are highly complex featuresthat need to be modeled numerically in order to simulate realistic fusion devices. Heat ﬂux through the
SOL must be regulated so as to limit heat loads onthe divertor plates high enough to cause physical damage, necessitating lengthy and expensive maintenance.
The most common approach is to model a 2d3v plasma
with computationally fast ﬂuid codes, but this requiresthat many kinetic factors be implemented manually,such as boundary conditions at the wall, heat conductivity, and plasma viscosity. Valid quantities must be
acquired from fully kinetic simulations.
To accurately calculate energy-dependent sputtering rates, both the heat ﬂux and velocity distribution function (VDF) at the sheath edge must be
known. The parallel heat ﬂux can be divided into
the ion and electron conductive and convective heat
ﬂuxes,
��
(1)
qconv,σ = nσ mσ V||,σ v 2 /2,
�
�
2
(2)
qcond,σ = nσ mσ (v|| − V||,σ )v /2,

qSH,σ = −κSHσ ∇|| Tσ= 3.9ni vti λii
κSH
i

(3)

κSH
e = 3.2ne vte λeewhere κSHσ is the heat conductivity, nσ is the density,
vtσ = (Tσ /mσ )1/2 is the thermal speed, λσσ = vtσ τσσ3is the thermal mean free path (MFP), and τσσ ∝ vtσ
is the same-species collision time. This expression isvalid when the actual heat ﬂux is less than the orderof the Maxwellian one-way free-streaming heat ﬂux
�
�1/2
,
(4)
qFS,σ = nσ Tσ T||,σ /mσ�where T||,σ = n−1
dv · v||2 fσ (x, v) is the temperatureσparallel to the magnetic ﬁeld. The collision-dependent
heat ﬂux is approximated via the harmonic mean [3],
�−1
�
1
1
+
,
(5)
qeﬀ,σ =qSH,σασ qFS,σwhere the ﬂux-limiting coeﬃcient ασ≡
limλmfp →∞ qσ /qFS is deﬁned as the ratio between
the actual heat ﬂux and the one-way free-streaming
heat ﬂux in the collisionless limit. The heat ﬂux
limiters are usually assumed to be ασ ≈ 0.1, althoughprevious kinetic simulations have given values in therange αi = 0.1 ∼ 2 and αe = 0.03 ∼ 3 [1].
We have studied the heat transport in SOL plasmas with the PARASOL code and shown that in acollisionless plasma where the ion-electron mass ratio is mi /me = 400, the electron heat ﬂux limiter is
αe = 0.75 in a plasma at equilibrium [4]. Further exploration showed that for a more realistic mass ratio ofmi /me = 1800, αe varies smoothly from a very smallvalue when the electron radiation energy-loss rate is
low to approximately unity when it is high [5]. In addition, examination of the eﬀect of collisionality has

where v is the particle velocity,� V = �v� is the ﬂuid
velocity, and brackets �a� ≡ afσ (x, v)dv indicate
an average of the variable a over the phase-space distribution function fσ of species σ ∈ e, i. Unlike theelectron heat ﬂux that is divided between convective
and conductive, the ion heat ﬂux is expected to bemainly convective [1]. The conductive heat ﬂux is approximated by the Spitzer-Härm expression [2] in a
author’s e-mail: aaron@riam.kyushu-u.ac.jp557©2010 by The Japan Society of Plasma
Science and Nuclear Fusion Research

A. Froese et al., Particle Simulation Study of SOL Heat Transport to the Divertor Plateshown that when there is no electron radiation, theheat ﬂux does not follow a harmonic average transition from the collisional to the collisionless state, but
peaks at moderate collisionalities when the mean freepath is of the order of the system length and then decreases to the collisionless limit. The reason for thedecrease was found to be the development of a large
high-energy tail in the electron energy distribution as
the plasma becomes more collisionless [6]. The tail effect is overpowered when electron radiation creates anasymmetry between the incoming and outgoing electron distributions, thereby raising the heat ﬂux in thecollisionless limit and eliminating the peak. Addinga Langevin heat bath model to the source region was
also found to suppress the tail asymmetry, preventing
formation of the peak.
These previous studies have concentrated on the
response of the electron heat ﬂux to collisionality, electron radiation rate, and the type of source model.
In contrast, this work investigates the behaviour of
the ion heat ﬂux in response to diﬀerent mass ratios
and recycling rates, as well as collisionality and source
model eﬀects.

There is no diﬀerence between inboard and outboard
divertor plates to avoid masking the heat ﬂux physics
we wish to observe. The half-domain shown in Fig.
1 is divided into three major regions: source, intermediate, and radiation/recycling. The source region
lies in the center, where hot particles are generated
and a Langevin heat bath may also be employed.
In the Langevin model, each particle j lying within
the source region experiences heating given by theLangevin equation,∆vj= −νvj + A,
∆t

� 2�
TL0 ν
,
A =2
mj ∆t

(6)

where ν is the Langevin relaxation constant, A is arandom variable of uniform distribution with range set
by the given temperature TL0 , and mj is the particle
mass. The heating parameters are ν∆t = 10−3 andTLi0 = TLe0 = Te0 . The Langevin model ensures thatthe velocity distribution tends towards Maxwellianeven in a collisionless plasma with poor mixing [6].
The radiation/recycling regions lie near the divertor
plates, where an electron energy sink and and an ambipolar cold particle source model radiative coolingand recycling, respectively. Explanation of the radiation model is given in [7]. In this case, the radiationenergy-loss fraction is set to frad = Qrad /Qsrc = 0.6.
Between the source and radiation regions lie the intermediate regions, where sampling occurs. As reported
by others [8], the values of ασ vary by orders of magnitude along the ﬁeld line, so the poloidal average overthe intermediate region is taken.
Ions lost to the divertor plates are replaced at
the hot particle and cold recycling sources. Particle
generation is ambipolar, so electron and ion source
ﬂuxes are equal, Γsrc,e = Γsrc,i , Γrec,e = Γrec,i . Since
the source particle ﬂux is equal to the ion ﬂux tothe divertor, a sheath potential forms and a positive
bias develops in the plasma bulk. The ratio of recycling to hot source particle generation is an input
parameter R = Γrec,σ /(Γsrc,σ + Γrec,σ ). Hot particles are generated uniformly within the source region Lsrc,a = 0.4L < x < L/2, with an isotropic
thermal distribution of temperature Te0 = Ti0 = 1.
Cold particles are generated uniformly within the recycling region Lrec,a = 0.01L < x < 0.21L = Lrec,b ,with an isotropic thermal distribution of temperature
Trec,e /Te0 = 5 × 10−3 , Trec,i /Te0 = 1 × 10−3 .

2. Simulation model
The PARASOL code simulates a 1d3v selfconsistent electrostatic particle-in-cell model with abinary collision model [4, 7]. A slab geometry is usedwith motion on the separatrix parallel to the magnetic ﬁeld. Position in the toroidal direction is given
by coordinate x (Fig. 1). Ions are fully traced (1d3v),
electrons follow their guiding centers (1d2v), and forthis study neutrals are ignored. The ratio of thepoloidal magnetic ﬁeld to the total ﬁeld strength isset to Θ = Bx /B = 0.2. The ion gyro-radius is set toρi /L = 5×10−3 , the number of spatial cells to 800, andthe number of ions to N0 = 106 . The time step length
in the runs presented is the inverse of the plasma fre−1. Each run requires approximately
quency, ∆t = ωpe
5
K = 1 × 10 time steps to reach equilibrium, whichtakes approximately 10 hours on 12 modern CPUs.
Equilibrium data is acquired by averaging over timefrom K = 1.8 × 105 to 2.0 × 105 .

3. Results and discussion
The eﬀects of four diﬀerent properties are explored: the mean free path is given a value of λmf p,i ∼L for a moderately collisionless SOL and λmf p,i ∼103 L for a very collisionless SOL, the recycling rate isgiven a value of 0% or 99%, the Langevin heat bathcan be absent or present, and the ion-electron mass
ratio is varied between mi /mp = 1 to mimic a hydro-

Fig. 1 Half the simulation domain showing source/sink regions, magnetic ﬁeld, and ﬂux values.

The spatial domain is symmetric across the midpoint, with divertor plates located at x = 0 and x = L.

558

A. Froese et al., Particle Simulation Study of SOL Heat Transport to the Divertor Plate

genic plasma and 4 to mimic an impure plasma, where
mp /me = 1800 is the proton mass. The electron and
ion VDF are averaged over the range x = [0 : 0.09]
which includes the sheath edge. We conﬁrm that theinclusion of the entire sheath region in the averaging
does not change the VDF signiﬁcantly, because theranges x = [0 : 0.09] and x = [0.09 : 0.18] give similar results. The eﬀects of the former three properties
are shown in Fig. 2 for the case with mi /mp = 1.
At ﬁrst glance, we see that the electrons are generally composed of two Maxwellian populations [5]: thebulk that are trapped by the sheath potential and ahigh-energy tail of electrons that travel directly fromthe hot source to the divertor plate. Despite not being trapped, the ions also have a bulk population andhigh-energy tail, which we have diﬃculty explaining,
especially in the case of a collisionless plasma.
By comparing the moderately collisional VDF in
Figs. 2 (a,c,e,g) to the corresponding collisionlessplasma in Figs. 2 (b,d,f,h), we see that the temperature of the bulk populations of both electrons andions is greatly reduced as mixing is reduced, since thetendency to return to a thermal distribution is eliminated. Less heat is transferred from the hot source
particles to the bulk as they travel to the divertor.
This leaves more and more electrons trapped in theradiation regions, reducing the electron heat ﬂux tothe divertor plate and hence the ﬂux at the source (asconﬁrmed in Fig. 3), decreasing the plasma potential.
At equilibrium, electron and ion loss is ambipolar, so
the ion ﬂow rate to the divertor also decreases, despite
reduced backscatter by the ions.
Fig. 2 (a,b,e,f) and (c,d,g,h) show a comparison of
the ion and electron VDF for the case of 0% recycling
versus 99% recycling. The addition of recycling reduces the ion temperature and ﬂow rate at the sheathedge, while concurrently increasing the electron temperature. The reduced ion temperature and ﬂow ratenaturally occurs because of the cold recycled ions entering the plasma in the radiation/recycling region.
The eﬀect is ampliﬁed as the collisionality increases
and the recycled ions cool more of those directly fromthe source. Just as electron radiation creates a lowenergy population of electrons and reduces the electron heat ﬂux, recycling creates a distinct populationand signiﬁcantly reduces the ion heat ﬂux to the divertor. However, if the origin of the bulk distribution iscool, recycled ions, it is unusual that the bulk temperature is many times larger than that of the recycled
ions, even without mixing. This must be investigatedfurther. Electron temperature increases because therecycled electron temperature is higher relative to therecycled ion temperature (Trec,e = 5Trec,i ) and theelectron radiation rate is very high (60%). Once electrons from the source reach the sheath region, they are
already cooler than the recycled electron temperature.

If the radiation model was deactivated, both ions andelectrons would experience an decreased temperatureand ﬂow rate due to recycling.
To see the eﬀect of the heat source model, we
compare heating using only the hot particle source in
Fig. 2 (a-d) with runs that also include Langevin heating in Fig. 2 (e-h). The Langevin heat bath increases
overall heating and forces the particles towards a thermal distribution. Since it aﬀects all particles within
the source region, mixing behavior occurs even in acollisionless plasma. In all cases, the ﬂuid velocity
and heat ﬂux towards the divertor plate increase, as
seen in Figs. 3 and 4. The density of particles in thehigh-energy tail also increases relative the total density. Due to the bath, the temperature in the source
region stabilizes at a high value, which increases thetotal heat ﬂux and reduces its susceptibility to otherfactors.
Heat ﬂuxes corresponding to each VDF in Fig. 2
with and without Langevin heating are shown in Figs.
3 and 4, respectively. These ﬁgures also show dependence of the heat ﬂuxes on the ion-electron mass ratio.
The mass ratio is found to introduce a small, smoothly
varying eﬀect on the heat ﬂux that typically varies
over only a factor of 2, equivalent to the change in(mi /me )1/2 . The only exception is the ion conductive
heat ﬂux, which increases by up to 1 order of magnitude as the ion mass grows from mi /mp = 1 to 4. A
greater ion mass should raise the sheath potential according to the relation eφe /Te ≃ 1/2 ln mi /me , whichtraps a larger fraction of the electrons and raises thebulk electron temperature. This further increases thesheath potential, which helps to accelerate the ions to
the divertor. However, we can see that ion conduction
is not always dominated by convection.
Other eﬀects on the heat ﬂux follow directly fromthe VDF. As the collisionality is reduced, when thermalizing behaviour is not alternatively provided byLangevin heating, the ion conductive heat ﬂux increases and electron heat ﬂux is reduced. This corresponds to the higher ion ﬂow velocity and reduced
electron temperature. Recycling reduces the ion temperature and increases the electron temperature, so
the heat ﬂuxes likewise decrease and increase, respectively. Langevin heating increases the particle temperature across the board, and so the heat ﬂuxes also
increase prodigiously.
The value of αi and αe for each set of conditions
is shown in Table 1. The left half shows the result ofusing only the hot particle source. In this case, thestrongest eﬀect is caused by collisionality. There is anincrease of 1 order of magnitude in αe and 2 orders in
αi when going from a moderately collisional plasma toa low collisional plasma. Without recycling, both limiters have a weak dependence on the mass ratio closeto ασ ∝ (me /mi )1/2 . However, with 99% recycling,559

A. Froese et al., Particle Simulation Study of SOL Heat Transport to the Divertor Plate(a) particle source only, λmfp = L, recycling=0%(b) particle source only, λmfp = 103 L, recycling=0%(c) particle source only, λmfp = L, recycling=99%(d) particle source only, λmfp = 103 L, recycling=99%(e) Langevin heat bath, λmfp = L, recycling=0%(f) Langevin heat bath, λmfp = 103 L, recycling=0%(g) Langevin heat bath, λmfp = L, recycling=99%(h) Langevin heat bath, λmfp = 103 L, recycling=99%

Fig. 2 Ion and electron distribution functions of random energy (ﬂuid velocity removed) averaged over the region
x/L=[0:0.09] for diﬀerent collisionalities, recycling rates, and either a hot particle or Langevin heat sources. Averagetemperatures (Ti , Te ), tail temperatures (Tt,i , Tt,e ), and ion ﬂuid velocities (vi ) are labelled.

560

A. Froese et al., Particle Simulation Study of SOL Heat Transport to the Divertor Plate

Fig. 3 Hot particle source
only case:massdependence of ion
conduction(+),
ion convection (×),
electron conduction(∗),
and electron
convection (⊡) heat
ﬂuxes measured at
the sheath edge.

Fig. 4 Langevinheating
case:mass dependenceofion
conduction(+),
ion convection (×),
electron conduction(∗),
and electron
convection (⊡) heat
ﬂuxes measured at
the sheathedge.
Negative values (circled) are plotted bymagnitude.

parametersRmi /mpλmfp /L

hot particle source only
αe
αi

1
1
1
1
1
1
1
1
103
103
103
103
103
103
103
103

1.0 × 10−1
7.0 × 10−2
5.5 × 10−2
4.6 × 10−2
1.4 × 10−2
1.0 × 10−2
8.7 × 10−3
7.3 × 10−3
7.4 × 10−1
5.1 × 10−1
4.0 × 10−1
3.5 × 10−1
7.8 × 10−2
3.9 × 10−2
2.4 × 10−2
3.1 × 10−2

0.00
0.00
0.00
0.00
0.99
0.99
0.99
0.99
0.00
0.00
0.00
0.00
0.99
0.99
0.99
0.99

1
2
3
4
1
2
3
4
1
2
3
4
1
2
3
4

4.5 × 10−3
2.2 × 10−3
1.8 × 10−3
1.4 × 10−3
3.2 × 10−3
3.9 × 10−3
4.0 × 10−3
4.2 × 10−3
2.4 × 10−1
1.8 × 10−1
1.4 × 10−1
1.1 × 10−1
2.2 × 10−1
4.2 × 10−1
5.3 × 10−1
5.5 × 10−1

Langevin heatingαeαi
3.9 × 10−1
3.8 × 10−1
3.5 × 10−1
2.7 × 10−1
2.8 × 10−1
1.7 × 10−1
1.3 × 10−1
1.1 × 10−1
7.2 × 10−1
5.4 × 10−1
3.7 × 10−1
3.1 × 10−1
1.1
6.7 × 10−1
4.7 × 10−1
3.6 × 10−1

3.8 × 10−2
4.4 × 10−2
1.6 × 10−2
1.4 × 10−2
7.0 × 10−2
3.5 × 10−2
2.4 × 10−2
1.8 × 10−2
4.8 × 10−2
4.3 × 10−2
5.0 × 10−2
6.5 × 10−2
1.8
1.6
1.4
1.3

Table 1 Heat ﬂux limiters ασ measured at x/L = [0.27 : 0.36] for diﬀerent collisionality, recycling rate, ion-electron mass
ratio, and source model.

561

A. Froese et al., Particle Simulation Study of SOL Heat Transport to the Divertor Plate

4. Summary

this trend reverses itself for αi , switching instead to asmall increase with mass ratio. The inclusion of recycling decreases αe by about 1 order of magnitude.
The right half of Table 1 shows the result ofLangevin heating. The overall eﬀect is to increase andstabilize the heat ﬂux limiters. Without the Langevinheating, αi lies in the range [1.4×10−3 : 0.55], but withheating, the range becomes [1.4 × 10−2 : 1.8]. Likewise for the electron heat ﬂux limiter, αe goes from therange [7.3 × 10−3 : 0.35] to [0.11 : 1.1]. In both cases,
the maximum and minimum increase by more thanan order of magnitude and their ratio becomes lessthan half as much as for the hot particle source only
case. Because of the additional heating, collisionality becomes less important, with the collisionless SOL
giving a slightly higher value for both limiters by asmall factor ∼ 1. Going from small to large mass ratio (mi /mp = 1 to 4) causes αe to decrease, and causes
αi to decrease in two cases: when there is neither recycling nor Langevin heating and when there is bothrecycling and Langevin heating. This is because thepresence of either recycling or Langevin heating creates a strong temperature gradient, making the ion
conductive heat ﬂux very sensitive to the mass ratio.
Since the ﬂux limiters ασ are deﬁned as the ratioof the actual heat ﬂux to the one-way free-streaming
heat ﬂux according to Eq. 5, they are most aﬀected bychanges to the shape of the actual velocity distribution. Looking at the VDF in Fig. 2, we see the distributions typically become more symmetric by decreasing the recycling rate, increasing collisionality, or using Langevin heating. The electron distribution reacts
somewhat more strongly to the changes than the ion
distribution. We can deduce that ασ is increased byconditions that cause an irregularly-shaped VDF: increasing the recycling rate, decreasing collisionality, or
deactivating Langevin heating. Table 1 conﬁrms that
this assumption is correct for the collisionality condition and for the heating method condition, except forthe collisionless case without recycling where both αeand αi become much larger when Langevin heating is
activated. In fact, it is seen that adjusting the plasma
collisionality is the most eﬃcient way to increase αiand especially αe . The recycling condition does notlend itself to such clear analysis. (Note that one cannot technically make a direct comparison between Fig.
2 and the Table 1, because the sampling regions arenext to the divertor (x/L = [0 : 0.09]) and inside theintermediate region (x/L = [0.27 : 0.36]), respectively.
However, only the distribution of electrons proceeding
away from the divertor is signiﬁcantly diﬀerent, andusually only at low velocity.)

To better evaluate the heat load on the divertor plates, the ion and electron VDF, as well as theion and electron heat ﬂux and heat ﬂux limiters αiand αe were examined. Both ion and electron VDF
were found to have high-energy tails separate fromlow-energy bulk populations, but for diﬀerent reasons.
The tails derive from the hot particle source, while theelectron bulk is due to trapping by the sheath. The ion
bulk is strengthened by recycling, but also occurs in
the absence of recycling. Collisionality and Langevinheating are found to force the plasma towards a thermal distribution, but the former reduces ﬂow rate andheat ﬂux, while the latter increases them. Recycling
reduces the ion temperature, ﬂow rate, and heat ﬂux.
Finally, the eﬀects of the ion-electron mass ratio aresmall, except when a strong temperature gradient canbe induced with recycling or Langevin heating.
The average values of the ﬂux limiters ασ arefound to be approximately 0.1, as currently employed
in ﬂuid codes. However, they show a wide deviationfrom 1.4 × 10−3 to 2 for only the tested simulation parameters, so the actual variation is likely even greater.
Therefore, using a constant value is not recommended.
At the very least, some linear correction should be
employed to account for changes to the ion-electronmass ratio, recycling rate, collisionality, and heating
model. As already pointed out [8], ασ even varies
by orders of magnitude along the ﬁeld line, a detailwhich we corroborate. Ultimately, limiting the heatﬂux may provide no substantial beneﬁts. Therefore,
future work should be directed to developing a moreaccurate model for incorporation in ﬂuid codes. In
addition, ﬂuid codes require a number of other kinetic factors, including the heat transmission coeﬃcients, ion viscosity, sheath potential, and polytropic
constant. Further inquiry will be directed towards theeﬀect of the plasma parameters on these factors.
[1] W. Fundamenski. Plasma Phys. Control. Fusion,
47(11):R163–R208, 2005.
[2] S. I. Braginskii. Transport processes in a plasma, volume 1. Consultant Bureau, New York, 1965.
[3] Peter C. Stangeby. The Plasma Boundary of Magnetic
Fusion Devices. Taylor & Francis, New York, 2000.
[4] T. Takizuka, M. Hosokawa, K. Shimizu. Trans. Fusion
Technol., 39:111–118, 2001.
[5] A. Froese, T. Takizuka, M. Yagi. Electron parallel heattransport in the scrape-oﬀ layer using a particle-in-cellcode. to be published in Plasma Fusion Res.
[6] A. Froese, T. Takizuka, M. Yagi. Eﬀect of source andsink on heat transport in the SOL. to be published inContrib. Plasma Phys.
[7] T. Takizuka, K. Tani, M. Azumi, K. Shimizu. J. Nucl.
Mater., 128:104–110, 1984.
[8] D. Tskhakaya et al. Contrib. Plasma Phys., 48(13):89–93, 2008.

562