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

Modeling of neutral transport for self-consistent transport
simulations in tokamaks
Mitsuru HONDA, Tomonori TAKIZUKA, Atsushi FUKUYAMA1) and Katsuhiro SHIMIZU
Naka Fusion Institute, Japan Atomic Energy Agency, Ibaraki 311-0193, Japan1)
Graduate School of Engineering, Kyoto University, Kyoto 606-8501, Japan(Received: 26 October 2009 / Accepted: 29 January 2010)

A model of the neutral transport in tokamaks has been developed suitable for a one-dimensionaltransport code, TASK/TX. The model consists of the one-dimensional diﬀusion equations for slow, thermal and halo neutrals, which are simultaneously solved together with two-ﬂuid equations and Maxwell’sequation. The behavior of the slow neutrals is easily modelled by a one-dimensional diﬀusion equation
owing to short mean free path, while that of the thermal neutrals with long mean free path is notstraightforwardly. The development of the way to evaluate an eﬀective mean free path reﬂecting their
two-dimensional motion and velocity realized the reasonable estimation of the thermal neutral diﬀusivity for a one-dimensional diﬀusion modeling. The validity of the model has been conﬁrmed by a Monte
Carlo code in the cases with various plasma parameters.
Keywords: neutral transport, mean free path, diﬀusion modeling, transport simulation, tokamak
dynamics occur near the wall, as a surface-plasma interaction. For instance, the dissociation process ofhydrocarbons has been included in detail in a moresophisticated SOL-divertor integrated code [6]. However, we now develop a neutral transport model so as
to precisely study the behavior of the neutrals as aparticle source, and in that sense a detailed consideration of the surface-plasma interaction is outside of
the scope of our study. We just assume that cold neutral atoms with 3 ∼ 5 eV are born from the wall after
the dissociation processes of molecules.
A Monte Carlo approach provides a ﬂexible
and convenient framework for developing numerous
physics calculations in terms of neutrals [7]. It can easily include the countless energy groups of neutrals andthe geometrical eﬀect and accurately evaluate a neutral density proﬁle within a statistical error. Monte
Carlo codes have been widely used to analyze the neutral population in experiments and to predict the neutral transport in simulations, over the entire plasma
including the core region [6, 7, 8, 9, 10, 11]. The Monte
Carlo approach is, however, not always suitable fora coupling with a time-dependent transport code because it requires much computation cost and its solution inevitably includes a statistical error. From the
aspect of the compatibility of a ﬂuid-type transport
code and the long-time simulation, a ﬂuid modelingof the neutral transport is still considered to be important.
There have been many neutral transport modelings except a Monte Carlo approach: analytic models for the estimation of the neutral density proﬁle[12, 13, 14] and computational diﬀusion models in a
two-dimensional geometry [15, 16, 17]. These modelings are essentially based on the assumption of short

1. Introduction
A prediction of density proﬁles in coming burning plasmas is one of the most important issues because the fusion output increases with the peakedness of the density proﬁle. The density proﬁle is determined by the particle transport and the particle
source. While the particle transport is usually governed by the turbulence-induced diﬀusion and inward
pinch [1], the source for charged particles in a plasma
has strong ties with the behavior of neutral hydrogen atoms and plays an important role especially inthe formation of the density proﬁle in the edge region.
Since the neutrals cannot be conﬁned by the magnetic ﬁeld, for a transport modeling they should betreated in a manner diﬀerent from charged particles.
In this paper, we model the behavior of the neutrals
as a set of equations suitable for a one-dimensional
multi-ﬂuid transport code, TASK/TX [2]. TASK/TX
has been used to study the ripple-induced toroidal rotation [3, 4] and to analyze the characteristic of the� torque induced by charge separation of fast par�j × B
ticles generated by a neutral beam (NB) injection [5].
The detail of the code will be described in the next
section.
In an actual plasma, there are countless energy
groups of neutrals: they gradually penetrate into
the core plasma, suﬀering from multi-step charge exchange. Actually from the aspect of ﬂuid modelingfor neutrals, however, only ﬁnite energy groups can behandled. In our modeling, we consider three groups ofneutrals and a ﬁrst-step charge exchange. We need to
conﬁrm the validity of the assumptions made in ourmodeling by using a Monte Carlo calculation.
Several kinds of atomic processes and molecularauthor’s e-mail: honda.mitsuru@jaea.go.jp529©2010 by The Japan Society of Plasma
Science and Nuclear Fusion ResearchM. Honda et al., Modeling of Neutral Transport for Self-Consistent Transport Simulations in Tokamaks

mean free path of neutrals, i.e. frequent charge exchange. In the region inside the edge region in presentday tokamak discharges, however, thermal neutrals
could have long mean free path. It is then diﬃcult
to directly apply these modelings to a wide rangeof plasma parameters, with which transport codes
must deal. We therefore would like to develop a onedimensional neutral transport model, compatible withTASK/TX and applicable to a wide range of plasma
parameters.
In this paper, a one-dimensional neutral diﬀusion
model in cylindrical coordinates is proposed for thethree-group neutrals, on the basis of the diﬀusion approximation idea [15]. The rest of this paper is asfollows. The next section describes a brief summaryof TASK/TX, and a neutral transport model consistent with TASK/TX is shown in section 3. Simulation
results for the validation of our model are given in section 4. Finally we summarize our results in section 5.

2. Multi-ﬂuid
TASK/TX

transport

There are several kinds of particle sources in the
TASK/TX system. Thermalization of fast beam ionsis the major particle source of the bulk ions. Furthermore, ionization of neutrals is also the major source,
usually localized near the edge region. Neutrals in atokamak plasma have many diﬀerent origins, such as
recycling, gas puﬃng, NB injection and so on. An accurate estimation of the particle source is required to
study the behavior of the particle transport and theresultant density proﬁles; hence these processes should
be carefully taken into account in the code.
To that end, TASK/TX solves the slowing down
equations for fast ions and the neutral diﬀusion
equations together with the two-ﬂuid equations andMaxwell’s equations. We should deal with the neutraltransport equations in a manner diﬀerent from otherequations. This is because the motion of neutrals is
not restricted by the magnetic ﬁeld unlike charged
particles. The usual diﬀusive transport equations arethen not directly applied for the neutral transport
analysis.

code,

We brieﬂy summarize the main characteristicsof the one-dimensional multi-ﬂuid transport code
TASK/TX with an emphasis on the diﬀerence fromconventional transport codes. More details about thecode and actual expression of the equations are foundin [2].

3. Modeling of neutral transport
In order to develop a one-dimensional neutral diffusive transport model, we should reduce an essentially three-dimensional eﬀect to a one-dimensional
one. As with most of the neutral transport modelings, we assume that neutrals move around in a twodimensional plasma cross-section, without considering
their movement in the toroidal direction.
In tokamak plasmas, there are many kinds of neutrals from the aspect of their energy and their origin.
Neutral atoms with cold temperature of typically 3
to 5 eV, after dissociation of hydrogen molecules, aremainly provided by the wall through recycling and/or
gas puﬃng. We call them slow neutrals hereafter.
They are randomly moving around, and subsequently
some of them are ionized and the other turn into
thermal neutrals with temperature equivalent to thatof the thermal ions with which the neutrals collide,
due to charge exchange. The so-called ﬁrst-generation
thermal neutrals have long mean free paths and areable to reach the core region, leading to the possible
particle source there. In a NB-heated plasma, the situation surrounding the neutrals is somewhat diﬀerent.
For example, a positive-ion based NB unit in JT-60U
provides high energetic neutrals with around 80 keV,much higher than the temperature of thermal ions.
Due to charge exchange between these neutrals andthe thermal ions, the neutrals turn into the fast ionswhile the thermal ions become the thermal neutrals,
which are called halo neutrals. The halo neutrals are
similar to the above-mentioned ﬁrst-generation thermal neutrals in that they have thermal velocities comparable to the background ion temperature. On the• The code solves the continuity equations, thethermal transport equations and the two-ﬂuid
equations of motion for electrons and ions coupledwith Maxwell’s equations as well as the equations
for neutrals and fast beam ions in the cylindrical coordinates (r, θ, φ), where r, θ and φ denote
the radial, poloidal and toroidal directions respectively.
• Since the multiple continuity equations for all
charged particle species are solved as well as Poisson’s equation, an explicit quasi-neutrality condition need not be imposed.
• Neoclassical eﬀects such as the bootstrap current
and the Ware pinch appear through the neoclassical viscous force term in the poloidal equation of
motion. The neoclassical viscosities are evaluatedby the NCLASS module [18].
• The calculation domain is extended to the scrapeoﬀ layer across the plasma surface (separatrix)
and is bounded by the wall: It enables us not to
impose boundary conditions on the plasma surface.
• The formation of the radial electric ﬁeld is accurately evaluated through Poisson’s equation andthe radial force balance in a manner consistent
with the ﬂows and background plasma proﬁles.

530

M. Honda et al., Modeling of Neutral Transport for Self-Consistent Transport Simulations in Tokamaks

other hand, they deﬁnitively diﬀer in where they areborn. The former can be born anywhere in the pathof NB, even in the core region, while the latter can begenerated only in the edge region where the cold neutrals can exist. In order to develop a neutral transportmodel suitable for TASK/TX, we separately treat the
behavior of these three groups of neutrals.
A time-dependent basic diﬀusion equation forneutrals is expressed by
�
�
∂n0j
1 ∂
∂n0j
rD0j
=
∂t
r ∂r
∂r
1
− ne �σiz v� n0j + S0i ,
(1)
Zi

2D

where the subscript 0j denotes the neutral species j:
j = 1 for the slow neutrals, j = 2 for the thermal
neutrals and j = 3 for the halo neutrals. n denotes thedensity, Zi the charge number of bulk ions, �σiz v� theMaxwellian averaged ionization cross-section and S0j
the source and sink. The diﬀusivity D0j for neutrals
is deﬁned by [16]
D0j ≡

2
v0j
λ0j v0j=,33ni (�σiz v� + �σcx v�)

r

O

1D
O

λ

r

Fig. 1 Illustration of a wrong modeling of ﬁrst-generation
neutral trajectories, denoted by slim arrows, bornat the magnetic surface r. The magnetic surface isdenoted by the broken line and the magnetic axis
by O.

(2)

where λ0j denotes the mean free path, v0j theMaxwellian averaged thermal velocity and σcx thecharge exchange cross-section.
The mean free path of the slow neutrals is muchsmaller than those of other kinds of neutrals becauseof their slow thermal velocity. In the edge region withlow electron temperature, they quickly suﬀer from ionization, not long after their birth at the wall. Therefore the population of the slow neutrals is highly localized in the edge region. Their density signiﬁcantly
decreases with distance from the wall and almost vanishes around typically ρ ≡ r/a ≈ 0.9, although theboundary varies depending on the background electron temperature and other conditions. Here r is the
local minor radius and a the minor radius. Owing to
the short mean free path, i.e. frequent collisions, thebehavior of the slow neutrals can be properly modelled by a one-dimensional diﬀusion equation. For theslow neutrals, the source and sink term consists of thecharge exchange loss, recycling and gas puﬃng as follows:of the thermal neutrals if we directly utilized eqs. (1)and (2) without any special treatment in terms of thelong mean free path. Since the source of the thermal
neutrals are localized in the edge region as explained
above, the thermal neutral density decreases with distance from there. In a diﬀusive manner, the neutrals
tend to move towards region of the lower density, thatis, the magnetic axis from the high density region. If
we model their movement in a two-dimensional space
as a one-dimensional diﬀusion equation (1) without
careful consideration, the movement towards the magnetic axis in the radial direction is regarded as if theneutrals intentionally funneled towards the magnetic
axis in two-dimensional space, as seen in ﬁg. 1. In
this case eq. (2) yields the excess estimation of D02 ,predicting the excess neutral density in the core. Of
course this cannot happen in an actual plasma. In
fact the neutrals born at a certain ﬂux surface overall move inwards as mentioned above, while they canmove any direction. This eﬀect should be appropriately included in the one-dimensional modeling.
Here we consider a thermal neutral born at theﬂux surface n corresponding to r, as shown in ﬁg. 2, inorder to estimate the more realistic diﬀusivity D02 (r),which can be applied to eq. (1). We assume thatit moves its mean free path length, labelled by λn ,and then arrives at the surface i. Since in this case we
consider only the neutral that moves inside the surface
n, a maximum angle is deﬁned as θn . The projection

rcyl
Puﬀ
+ S01
.
S01 = −ni �σcx v� n01 + S01

On the contrary, the mean free path of the ﬁrstgeneration thermal neutrals is rather long because of
their high temperature due to charge exchange withthermal ions. It could be half of the minor radiusa. The thermal neutrals are born in the region where
the slow neutrals exist, and they can penetrate deep
into the core due to their long mean free path. We
would fail to reasonably estimate the density proﬁle531

M. Honda et al., Modeling of Neutral Transport for Self-Consistent Transport Simulations in Tokamaksnweighted volume-averaged temperature is a representative of the various thermal neutral temperatures, asdeﬁned below:�bn02 (r)Ti (r) rdr,
T̄02 ≡ r01� bn (r) rdr
r01 02magnetic surface
i
i-1
λn

m
O

θi
λi

prj

θn

where r01 is the innermost boundary of the region
where the slow neutrals exist and b denotes the radiusof the wall. �By using the Maxwellian averaged velocity v02 = 8k T̄02 /(πmi ) and the locally evaluatedcross-sections, we obtain λ02 appearing in eq. (3) andthen D02 by eq. (2). Here k denotes the Boltzmann
constant and mi the ion mass. The source term canbe expressed asrλn

S02 = ni �σcx v� n01 ,Fig. 2 Illustration of the computation of an eﬀective mean
free path of a ﬁrst-generation neutral at a positionr, corresponding to the magnetic surface denoted
by n. A mean free path is denoted by λ and anangle between the intersection points of a magnetic
surface and a neutral trajectory by θ. The subscripts correspond to the labels of the magnetic surface.

which is the particle supply from the slow neutrals due
to charge exchange.
Unlike the ﬁrst-generation thermal neutrals, thehalo neutrals are usually produced in the core region
near the magnetic axis because a source proﬁle for thehalo neutrals coincides with a birth proﬁle due to NB
heating, i.e.
NB.
S03 = S03of λi (= λn ) on the midplane is λprj
= λn cos θi . If
i
as
the
mean
free
path
of the neutrals
we regard λprjiwhich are born at n and fall into the area between i
and i − 1, the ratio of them to all the neutrals bornat n is (θi − θi−1 )/θn . Some of these neutrals may
arrive at the surface n after moving their mean freepath, meaning that they eventually remain at the born
surface n. Based on this concept, an eﬀective mean
free path used in the one-dimensional modeling must
be shorter than λn (= λ02 ). Since this argument canbe applied from the surface m to n, where the surface
m denotes the farthest surface at a distance of λn fromthe surface n, we ﬁnally obtain the eﬀective mean freepath λeﬀ
02 (r) evaluated at r in the form:
λeﬀ
02 (r) =
=

n
�

i=m+1
n
�
i=m+1

λprj
i

The typical density proﬁle of the halo neutrals is centrally peaked and gradually decreases towards theedge region, although it depends on the birth proﬁle. Then they mainly tend to diﬀuse from the coreto the edge, with the velocity faster than that of thethermal neutrals. It means that their mean free path
is longer than that of the thermal neutrals. Owing
to the core-localized source proﬁle and the fairly longmean free path, the halo neutrals can be regarded asa simple outward diﬀusion in cylindrical coordinates,
and in addition many of them may be ionized in thecore due to the high plasma density. Therefore we canreasonably use eqs. (1) and (2) for them.

θi − θi−1
θnλ02 cos θiθi − θi−1.
θn

4. Validation of the neutral transportmodeling
In the preceding section, we developed a set ofdiﬀusion equations for three-group neutrals based on
the physical consideration. We should then conﬁrm
the validity of our modeling by comparing a neutral
density proﬁle with that by a Monte Carlo code TOPICS/NT [9]. In this section, we focus our attention
on the behavior of the slow and thermal neutrals onlyin an ohmic plasma without NB heating. This is because we devoted our eﬀort to the development of thetransport model for the thermal neutrals.
We use JT-60U-like plasma parameters: the major radius R0 = 3.2 m, a = 0.8 m, the toroidal magnetic ﬁeld Bφ = 2.68 T, the plasma current Ip =(3)

We note that we assume no poloidal variation of theneutral density on a ﬂux surface in this estimation.
All the ﬁrst-generation neutrals are assumed tohave mono-energy in our modeling. When we estimate λ02 for D02 , we must properly determine their
thermal velocity v02 , i.e. their temperature. The
source region of the thermal neutrals exactly coincides
with the region where the slow neutrals exist. Evenin this narrow edge region, the ion temperature Ti
signiﬁcantly varies. Then we assume that the density-532

M. Honda et al., Modeling of Neutral Transport for Self-Consistent Transport Simulations in Tokamaks
1e+16

1e+16

1e+15

1e+15

T e=Ti [keV]

ne [1020 m-3]

1e+14

1e+14

1e+13

1e+12

T e=Ti [keV]

ne [1020 m-3]

TASK/TX w/o

Total neutral density
0

0.2

0.4

r/a

0.6

0.8

TOPICS/NTTASK/TX

1e+13

TOPICS/NTTASK/TX

1e+12

1Fig. 3 Comparison of the neutral density proﬁle calculated
by a Monte Carlo code TOPICS/NT with that by
TASK/TX in a transient phase at t = 0.1 s, in thelower plasma density case.

Total neutral density
0

0.2

0.4

r/a

0.6

0.8

1

Fig. 4 Comparison of the neutral density proﬁles calculated by TOPICS/NT, TASK/TX and that without any special treatment for the thermal neutral transport modeling in a quasi-steady state at
t = 0.5 s, in the lower plasma density case.
1e+16

1.0 MA, the gas puﬀ rate Γ0 = 1.0 × 1020 m−2 s−1 andthe recycling rate γ0 = 0.8. For simplicity, in the following simulations we assume that temperature proﬁles for electrons and ions are the same and are ﬁxed atinitial proﬁles without loss of generality. TASK/TX
calculates the particle conﬁnement time τp at every
time step, and then for comparison it provides TOPICS/NT with the calculated τp as a constraint condition.
First of all, we study the case in a lower density plasma, as observed in typical JT-60U ohmic discharges. In a transient phase at t = 0.1 s, we comparea neutral density proﬁle calculated by TOPICS/NT
with that by TASK/TX, as shown in ﬁg. 3. At
this time the volume averaged electron density n̄e is
9.86 × 1018 m−3 , and the background proﬁles of theelectron density and the temperatures are also seenin the ﬁgure. The neutral transport model developed well reproduces the neutral density proﬁle withboth the slow and thermal neutrals mixed predicted
by TOPICS/NT. Roughly speaking, the slow neutrals
dominates the neutral density proﬁle outside ρ = 0.9
and the thermal neutrals does inside ρ = 0.9. In bothregions, the slope of the proﬁle agrees well between twocases. The maximum deviation between both proﬁles
is just by a factor of ∼ 1.3. Fig. 4 shows the comparison of the proﬁles in a quasi-steady state at t = 0.5 s.
We again have a very good agreement between two
predictions even in the steady state. The tendencyof the proﬁle is almost similar to the transient case.
Just for reference, the ﬁgure also includes the result by
TASK/TX without any special treatment for the thermal neutral transport modeling. More speciﬁcally, inthis case λ02 is adopted instead of λeﬀ
02 , and the local thermal neutral temperature Ti (r) is used instead

1e+15
T e=Ti [keV]

ne [1020 m-3]

1e+14

1e+13

1e+12

1e+11

TOPICS/NTTASK/TX

Total neutral density
0

0.2

0.4

r/a

0.6

0.8

1

Fig. 5 Comparison of the neutral density proﬁle calculated
by TOPICS/NT with that by TASK/TX in a transient phase at t = 0.5 s, in the higher plasma density
case.

of T̄02 . This model clearly overestimates the thermal neutral density in the core, indicating that themodel developed in this study plays a crucial role inthe correct estimation of the neutral density proﬁle.
We conclude that the neutral density prediction by
TASK/TX agrees well with that by TOPICS/NT, regardless of whether the plasma is in a transient phaseor a steady state.
It is meaningful to check the model validity ina diﬀerent operation regime. We then compare neutral density proﬁles in a higher density plasma witha little high temperature compared with the previous
cases. A comparison result with background density
and temperature proﬁles is shown in ﬁg. 5. In thiscase the ohmic plasma with n̄e ≈ 1.73 × 1019 m−3 is533

M. Honda et al., Modeling of Neutral Transport for Self-Consistent Transport Simulations in Tokamaks

still in a transient phase at t = 0.5 s. The density proﬁles are in good agreement with each other again. In
comparison with the previous cases, a slight disagreement in the shape of the proﬁles is found, while themaximum deviation is very small by a factor of ∼ 1.4.
These results show that the model works well in plasmas with both low and high densities, indicating that
it is applicable to the wide range of plasma parameters, which is a very important feature for the model
cooperated with a transport code.

[1] C. Bourdelle, Plasma Phys. Control. Fusion 47, A317
(2005).
[2] M. Honda and A. Fukuyama, J. Comput. Phys. 227,
2808 (2008).
[3] M. Honda, et al., Nucl. Fusion 48, 085003 (2008).
[4] M. Honda, et al., J. Plasma Fusion Res. SERIES 8,
316 (2009).
[5] M. Honda, et al., Nucl. Fusion 49, 035009 (2009).
[6] K. Shimizu, et al., Nucl. Fusion 49, 065028 (2009).
[7] M.H. Hughes and D.E. Post, J. Comput. Phys. 28, 43(1978).
[8] K. Shimizu, et al., J. Nucl. Materials 313-316, 1277
(2003).
[9] N. Hayashi, et al., J. Phys.: Conf. Ser. 123, 012025(2008).
[10] D. Heifetz, et al., J. Comput. Phys. 46, 309 (1982);
http://w3.pppl.gov/degas2/.
[11] D. Reiter, et al., Plasma Phys. Control. Fusion 33,
1579 (1991); D. Reiter, J. Nucl. Materials 196-198,
80 (1992); http://www.eirene.de/.
[12] R.J. Goldston, Plasma Phys. 20, 1199 (1978).
[13] P. Bachmann, Physica Scripta 38, 835 (1988).
[14] D.F. Düchs, D.E. Post and P.H. Rutherford, Nucl.
Fusion 17, 565 (1977).
[15] E.L. Vold, A.K. Prinja, F. Najmabadi and R.W.
Conn, J. Nucl. Materials 176 & 177, 570 (1990);[16] E.L. Vold, F. Najmabadi and R.W. Conn, Nucl. Fusion 32, 1433 (1992).
[17] A.K. Prinja and C. Koesoemodiprodjo, J. Nucl. Materials 196-198, 340 (1992).
[18] W.A. Houlberg, K.C. Shaing, S.P. Hirshman and
M.C. Zarnstorﬀ, Phys. Plasmas 4, 3230 (1997).

5. Summary and discussion
We have developed the one-dimensional ﬂuidtype neutral transport model well compatible withTASK/TX. The main characteristic of the model isto separate all the neutrals into three groups from the
aspect of their energy and their origin: the slow, thermal and halo neutrals. The diﬀusivity for the thermal neutrals with long mean free path, mainly moving from the edge to the core, is carefully modelledconsidering the neutral motion in a two-dimensional
space, as shown in ﬁg. 2. The neutral density proﬁles predicted by TASK/TX with the model developed were compared with that by a Monte Carlo code
TOPICS/NT in the cases with the low or high density
in the transient or steady state phase, and the very
good agreement has been obtained. It means that themodel is useful and suitable for cooperation with atime-dependent transport code, which must cover thewide range of plasma parameters.
There is room for improvement of the model. If it
were applied to a plasma with a signiﬁcantly high density and a low temperature near the edge region, themodel prediction would deviate from a Monte Carloprediction. In such a plasma the mean free path ofﬁrst-generation thermal neutrals is shorter than thatin the above-considered cases. Due to the rather shortmean free path, these neutrals may undergo second
charge exchange in the outer core region before ionization and then the resultant second-generation thermal
neutrals will have higher temperature compared withthe ﬁrst-generation one, leading to longer mean freepath. If we introduce the second-step or multi-stepcharge exchange process into the model, its predictions will agree with Monte Carlo calculations withhigh accuracy even in such an extreme case.

Acknowledgements
One of the authors (MH) would like to thank Dr.
T. Ozeki for his continuous encouragement and support. This work was supported by the Grant-in-Aid
for Young Scientists (Start-up) (No. 20860092) and
Scientiﬁc Research (S) (No. 20226017) from Japan
Society for the Promotion of Science (JSPS).

534