Proceedings of ITC/ISHW2007

A Method for Ion Distribution Function Evaluation Using Escaping
Neutral Atom Kinetic Energy Samples
P.R. Goncharov, T. Ozaki, E.A. Veshchev, S. Sudo
National Institute for Fusion Science, 322-6 Oroshi-cho, Toki, 509-5292, Japan
A reliable method to evaluate the probability density function for escaping atom kinetic energies is required forthe analysis of neutral particle diagnostic data used to study the fast ion distribution function in fusion plasmas.
Digital processing of solid state detector signals is proposed in this paper as an improvement of the simple
histogram approach. Probability density function for kinetic energies of neutral particles escaping from the plasma
has been derived in a general form taking into account the plasma ion energy distribution, electron capture and loss
rates, superposition along the diagnostic sight line and the magnetic surface geometry. A pseudorandom number
generator has been realized that enables a sample of escaping neutral particle energies to be simulated for given
plasma parameters and experimental conditions. Empirical probability density estimation code has been developedand tested to reconstruct the probability density function from simulated samples assuming Maxwellian andclassical slowing down plasma ion energy distribution shapes for different temperatures and different slowing downtimes. The application of the developed probability density estimation code to the analysis of experimental dataobtained by the novel Angular-Resolved Multi-Sightline Neutral Particle Analyzer has been studied to obtain thesuprathermal particle distributions. The optimum bandwidth parameter selection algorithm has also been realized.
Keywords: neutral particle analysis, ion distribution, statistical data processing, empirical probability density, kernel
bandwidth selection
1. Introduction
Measurements of kinetic energy distributions ofneutral atoms escaping from magnetically confinedplasma are used in controlled fusion experiments as amethod to investigate the ion component distributionfunction and its evolution due to the application ofvarious plasma heating schemes. The ion distributionfunction reflects the kinetic effects, the single particle
confinement properties depending on the particularmagnetic configuration, the finite β effects such as MHD
induced fast ion losses, radial electric field effects, etc.
The nuclear fusion reaction rate is determined by the ion
distribution and thus its studies at suprathermal energies
near the rate coefficient curve maximum are of primary
importance. Advanced neutral particle diagnostics basedon solid state detectors with high energy resolution, e.g.
[1], are used to study the suprathermal ion distributionfunction. Statistical data processing is required to obtain asmooth normalized probability density function forparticle energies using the measured random samples [2].

2. Escaping Neutral Particle Energy Distribution
The probability density function (PDF) f(E) for
kinetic energies of neutral H0 particles escaping from theplasma of a magnetic confinement fusion device in ageneral form is given byauthor’s e-mail: pavel@nifs.ac.jp∫ Q ( ρɶ ) λmfp ( E , ρɶ ) d ρɶ 1
ρ−f ( E ) = Ae−1

∫ g ( E, ρ ) ×min

ρ min−1− ∫ Q + ( ρɶ ) λmfp
( E , ρɶ ) d ρɶQ + ( ρ ) e ρmin−ρ−1− ∫ Q − ( ρɶ ) λmfp
( E , ρɶ ) d ρɶ−Q ( ρ ) e ρmin dρρ(1)where A is the normalization constant. The source
function for H0 atoms of energy E within the plasmag ( E , ρ ) = ni ( ρ ) fi ( E , ρ )∑ l n(l ) ( ρ ) σ v

(l )

(2)

is expressed via the local plasma proton distributionni ( ρ ) f i ( E , ρ ) and the sum of rates over all targets fortheelectron

captureprocess. The derivatives
Q ( ρ ) = d Λ d ρ > 0 and Q − ( ρ ) = d Λ d ρ < 0 ofthe sight line distance Λ along the two intervals between
ρ = 1 and ρ = ρ min are obtained from the known+

structure of magnetic surfaces ρ = const . The neutral
flux attenuation enters in the form of Poisson exponents,
where lmfp(E, ρ) is the H0 mean free path with respect toall electron loss reactions.

Proceedings of ITC/ISHW2007

3. Numerical Experiment
Ideally, the passive diagnostic data is an array(E1, …, EN) of energies of escaped neutral particles
measured along a certain observation direction, and N is
the total number of particles collected during a certaintime interval. This array is a sample of realizations of therandom variable E distributed according to the law (1).
Such form of data is achievable with solid state detectors
by using pulse heighttechniques, while the other analysis
analyzers, e.g. E || B ones, intrinsically form ahistogram of the incoming particle energies over a certain
number of subintervals called energy channels. Technicaldetails may be found in [1]. The formulation of theproblem considered here is to obtain an estimate f(*)(E) ofthe unknown exact probability density function f(E) ofneutral particle energies from the experimental data. Thesought function preferably should satisfy a specified
precision criterion. The obtained PDF estimate is then to
be used to reconstruct the ion distribution for further
analysis.
Assuming a predefined theoretical PDF f(E) one cancarry out a numerical experiment by generating a sampleof escaped atom energies for given plasma parametersand experimental conditions. We apply the inverse
cumulative distribution function (CDF) approach. First, asample of pseudorandom numbers (u1, …, uN) uniformlydistributed within the 0,1) interval is generated using
an algorithm from [3]. Then, the energy values arecalculated as solutions of the equation[whereF (E j ) = u j ,(3)EF ( E ) = ∫ f ( Eɶ )dEɶ(4)0

is the CDF. These simulation results can be supplied as
input data for the PDF estimation procedure to test its
performance, since the original exact f(E) used in thesimulation is known.
Two typical ion energy distribution laws have beenused in the numerical simulation, namely, (a) Maxwellian
distribution with ion temperature Tifi ( E ) =

2 1

Eπ Ti

Ti

exp ( − E T i ) ,
γ 3, E ,π  2 Ti 

2Fi ( E ) =(5)(6)

xwhereγ (α , x ) = ∫ tα −1e−t dt0

is the lower

incomplete gamma-function; and (b) the classical slowing
down distribution for a delta-like fast ion source function−S ( v − v0 ) =( v − v0 )2

 v* ( v, t ) − v0 S0 τ s  v − v0   (8)− erf erf 33 e
8π v + vc  e  where the slowing down time τ s =(7)3m pTe3 24 2π ne e 4 Λme1 2,

the critical velocity vc3 =3 2π Te3 2, Λ is the Coulomb
2m p me1 2logarithm, and v* ( v, t ) =(( v + v ) eThe ion velocity v =33c3t / τ s− vc31/ 3)[4].

2 E / m p , v0 is the injection

velocity, S0 and e determine the source rate and width,
and t is the time. Fig. 1 (a) shows the Maxwellian PDF
for two different Ti values and Fig 1 (b) shows theclassical slowing down PDF at t = 0.8 s for injectionenergy E0 = 150 keV and two different pairs of the target
plasma ne and Te values. Histograms of the corresponding
pseudorandom number samples governed by these PDFs
are shown in Fig. 1 (c) and (d).

4. Data Processing Method
As an improvement of the neutral particle diagnostic
data analysis, we have applied the probability density
estimation using kernel smoothing techniques, e.g., [5, 6].
The kernel PDF estimate

f (K ) ( E ) =1 N  E − Ej ∑K ,Nh j =1  h h>0(9)is determined by the kernel function K(z) and the kernel
bandwidth h. The performance criterion of this method isthe value of the mean integrated squared error+∞

2MISE ( f ( K ) , f ) = ∫  f ( K ) ( E ) − f ( E )  dE , (10)0where averaging is over different samples of N
realizations (E1, …, EN), and its “asymptote” for N >> 1(large sample approximation)

c1 c2c32 h4,(11)+Nh42w h e r e c1 = ∫ K 2 ( z ) dz , c2 = ∫ ( f ′′( z ) ) dz a n d
AMISE ( f ( k ) , f ) =c3 = ∫ z 2 K ( z )dz . The optimum kernel derived in [7] isK ( z) =31 − z 2 ) I ( −1, 1) ( z ) , w h e r e(4I ( −1, 1) ( z )equals unity within (-1, 1) and equals nought outside.
However, it is emphasized [5, 6] that the choice of thekernel function shape has a small influence on the method
performance, while the bandwidth parameter choice ismore important. Therefore, Gaussian kernel

K ( z) =e2S0 e
4π v 2 e πfi ( v ) =1 −z2 2e2π

has been used, since it has continuous derivative.

(12)

Proceedings of ITC/ISHW2007(a) 10-1

n e = 0.7x1013cm -3 ; T e = 3.5 keV(b)
Ti = 5 keV
Ti = 10 keVn e = 4.5x1013cm -3 ; T e = 0.5 keV-1f(E), keV
1010τs = 1 s ; t = 0.8 s-1-2f(E), keV-1101010-210-3τs = 0.01 s ; t = 0.8 s-3E0 = 150 keV-40501001502000250501100CountsCounts(d) 10210050100150200

250150200250

E, keVE, keV(c) 101001502002101100250050100

E, keVE, keV-1(e) 10(f)Ti = 5 keV Kernel PDF-21010-210-3

τs = 0.01 s Kernel PDF(E), keV(K)(K)-310ffτs = 1 s Kernel PDF-1-110(E), keV-1

Ti = 10 keV Kernel PDF-410050100

150200250050E, keV

100150200250

E, keV

Fig. 1. (a) Maxwellian PDF for Ti = 5 keV (green) and Ti = 10 keV (blue); (b) classical slowing down PDF for τs = 1 s(green) and τs = 0.01 s (blue); (c) histograms of pseudorandom number samples distributed according to thelaws shown in (a); (d) histograms of pseudorandom number samples distributed according to the laws shown in (b);(e) kernel PDFs calculated from Maxwellian law pseudorandom number samples shown in (c); (f) kernel PDFs
calculated from slowing down distribution law pseudorandom number samples shown in (d).

A reliable practical method for optimum h selection
was proposed in [8] and revisited recently in [9]. The
bandwidth is sought by solving the equation15
1h = 
 ,2πNφ(η(h))
4


(13)

φr (h) =N N
 Ei − E j 
1K (r ) 
 . (15)( r +1) ∑∑N ( N − 1)hi =1 j =1
 h 

The function
17

where the function φr is expressed via kernel derivativeK ( r ) ( z ) = ( −1) r Her ( z ) K ( z )as follows(14) −6 2φ4 (a )  5 7
η (h) =  h φ6 (b) depends on the values(16)Proceedings of ITC/ISHW2007
10310ARMSNPA Detector Segm. #2ICH, NBI #1 185 keV NBI #3 151 keV
ne = 1.5x1013cm-3 ; Te = 3 keV
2101100ARMSNPA Detector Segm. #2
NBI #3 150 keV
13-3
ne = 0.5x10 cm ; Te = 1.75 keV

H/Ar target plasma

102101100H/He target plasmaCountsCounts

10350100150200250

50100150200

250

E, keVE, keV
Gaussian Kernel PDF(AMISE optimal bandwidth)

Gaussian Kernel PDF(AMISE optimal bandwidth)-110-1f (E), a.u.

f (E), a.u.

10-2

10-210-3-31010

50

100

150200250

50E, keV

100

150200250

E, keV

0Fig. 2. Experimental H energy spectra (upper) and PDF estimates (lower) for two different ion heating schemes.
17 16 2 a =  σˆ
 5N 

19

 480 2 ˆand b =  105 N  σ , (17)
NReferences[1]N11Ei − E ) , E = ∑ Ei . (18)(∑N i =1N − 1 i =1

[2]

Direct implementation of formula (15) is slow. An
approximate fast calculation technique is given in [9].

[3]
[4]

5. Application to Experimental Data

[5]

where σˆ =2

These methods have been tested by reconstructing theprobability density function from the generated pseudorandom
number samples assuming Maxwellian and classical slowing
down plasma ion energy distribution shapes for different
temperatures and different slowing down times. The test results
are shown in Fig. 1 (e) and (f). The analysis of passive
chord-integrated experimental data obtained with theAngular-Resolved Multi-Sightline Neutral Particle Analyzer [1]on Large Helical Device is illustrated in Fig. 2. Kernelsmoothing methods require a certain choice of the bandwidthparameter. An automatic choice described in Section 4 ispreferable for routine data processing.

[6]
[7]
[8]
[9]

E.A. Veshchev, T. Ozaki, et al., Rev. Sci. Instrum., 77,
10F129 (2006).
P.R. Goncharov, T. Ozaki, et al., Proc. 34th EPS Conf.
Plasma Phys., P5.081 (Warsaw, 2007).
S.K. Park, K.W. Miller, Comm. ACM, 31, 1192 (1988).
J.G. Cordey, M.J. Houghton, Nucl. Fusion, 13, 215(1973).
A.W. Bowman, A. Azzalini, Applied smoothing
techniques for data analysis, Clarendon (1997).
M.P. Wand, M.C. Jones, Kernel Smoothing, Chapman andHall, London (1995).
V.A. Epanechnikov, Theory Probab. Appl., 14, 153(1969).
S.J. Sheather, M.C. Jones, J. R. Statist. Soc. B, 53, 683(1991).
V.C. Raykar, R. Duraiswami, Very fast optimal
bandwidth selection for univariate kernel density
estimation,UniversityofMaryland,
CS-TR-4774/UMIACS-TR-2005-73 (2006).