Fusion Res.
Res. SERIES,SERIES, Vol.
Vol. 6 (2004) 000-000
618–621
J. Plasma Fusion

Computer Simulation of Fishbone Oscillation
SHIOZAKI Yutaka and TODO Yasushi1
Graduate University for Advanced Studies, Toki, 509-5292, Japan
1
National Institute for Fusion Science, Toki, 509-5292, Japan(Received: 9 December 2003 / Accepted: 9 August 2004)

Abstract
Particle-magnetohydrodynamic (MHD) hybrid simulations were carried out to investigate nonlinear evolution ofprecessional fishbone instability using parameters similar to the PDX experiment. Spatial profile of the fishbone mode
is different from that of the kink mode. The fishbone mode frequency shifts downward at saturation of the instability,
while it is close to the experimental value in the linear growth phase. The saturation level of the plasma displacement
is comparable to the radius of the q = 1 surface.
Keywords:
precessional fishbone instability, computer simulation, frequency chirping

1. Introduction

study. Electromagnetic field is given by the MHD
description. This approximation is reasonable under thecondition that the energetic ion density is much less than thebulk plasma density. The ideal MHD equations are,Fishbone instability was first observed during nearperpendicular neutral beam injection in the PDX tokamak[1]. The instability takes place in recurrent bursts,
accompanied by frequency decrease during each burst. Twotypes of fishbone modes, precessional fishbones [2] andthermal ion diamagnetic drift fishbones [3], weretheoretically predicted. Nonlinear simulations of the thermal
ion diamagnetic drift fishbones were carried out in Ref. [4]where a perturbative description was applied. On the otherhand, little is known about self-consistent nonlinear evolutionof the precessional fishbones. Perturbative approach is not
valid for the precessional fishbones, since the mode
frequency emerges from the shear Alfvén continuous
spectrum and depends intrinsically on the energetic ion
distribution. This makes the problem more difficult andinteresting than the thermal ion diamagnetic drift fishbones.
We carried out particle-MHD hybrid simulations withparameters similar to the PDX experiments. This is the firstsimulation of the precessional fishbone that demonstrates thesaturation of the instability and decrease in the mode
frequency.
In Sec. 2, simulation model and method are described.
Simulation results are reported in Sec. 3. Section 4 is devoted
to disccussion and summary.

2t= - d $ _ t vi,
2t
t 2 v + t v $ dv = - dp + n1 ^d # Bh # B,
2t
0
2B = - d # E,
2t
2p
= - d $ _p vi - _ c - 1i p d $ v,
2t
E =- v # B

(1)
(2)
(3)
(4)
(5)

where n 0 is the vacuum magnetic permeability, c is theadiabatic constant, t is the mass density, v is the velocity,and p is the pressure. The MHD equations are solved using afinite difference scheme of 4th order accuracy in space andtime.
The drift-kinetic description is employed for theenergetic ions which are generated by the neutral beam
injection. The guiding-center velocity u is given byu = v* + vE + vB,v
v* = B 8 B + t B d # b B,
1 E# b ,vE = B
6@vB = q1B 8- n d B # b B,
hmh v
t = qB ,hdv
mh v= v* $ 8 qh E - n dB B,dt

2. Simulation model
In the particle-MHD hybrid simulation model [5] plasma
is divided into two parts, bulk plasma and energetic ions. Thebulk plasma is described by the ideal MHD equations. Twofluid effects which might be important to the kink andfishbone instabilities are not taken into account in the present

(6)
(7)
(8)
(9)
(10)
(11)

©2004 by The Japan Society of Plasma
Science and Nuclear Fusion Research

Corresponding author’s e-mail: shiozaki@toki.theory.nifs.ac.jp1618

Shiozaki Y. and Todo Y., Computer Simulation of Fishbone Oscillation

where v is the velocity parallel to the magnetic field, n is
the magnetic moment, m h and q h are mass and charge of anenergetic ion, and b is a unit vector paralell to B. The effectsof the energetic ions on the bulk plasma are taken into
account in the MHD momentum equation:t 2 v+ t v $ dv = _ Q - Qh i E + c n1 d # B- jh m # B-dp , (12)2t0where Q and Q h are the total charge density and energetic-ion
charge density, and j h is the energetic-ion current density.
The total charge density Q is negligible in the MHD
approximation where the quasineutrality is satisfied. Paying
attention to the fact that - Q h E cancel out with the Lorentz
force due to E # B current of the energetic ions, Eq. (12) isrewritten tot 2 v + t v $ dv = c n1 d # B - j'h m # B - dp,
2t0Fig. 1 Initial beta profiles and q-profile.

(13)
radius is a = 0.43 m. The magnetic field at the axis is 1.5 T.
The plasma density is taken to be uniform 5 # 1019 m - 3 . The
bulk plasma is hydrogen and the energetic ion is deuterium.
The initial distribution of energetic ions in the velocity space
is a slowing down distribution with the maximum energy
50 keV. The slowing down distribution is taken to be
isotropic in the velocity space for simplicity. The spatial
profile of energetic ion pressure is a Gaussian with pressure
gradient scale length 0.15 m at the q = 1 surface, where q is
the safety factor. Two runs where the central energetic ion
beta values b h (0) are, 0 % and 3 %, respectively, werecarried out.
Figure 1 shows the initial beta profiles and q-profile forb h (0) = 3 % . Figure 2 shows spatial profiles of the radial
velocity v r in the linear growth phase. Solid and dashedcurves represent the sine and cosine part of the m/n = 1/1
harmonic, respectively. Here m is the poloidal mode number
and n is the toroidal mode number. Dashed-dotted curves
show the q -profile. In Fig. 2 the phase is chosen so as tomaximize the sine part amplitude. The well-known top-hat
structure of the kink mode is seen in Fig. 2(a). This mode is apurely MHD mode. For the central energetic ion beta valueof 3 % a fishbone mode is destabilized. We have investigatedenergy evolution. Only the energetic ion energy is decreasing
while the magnetic energy, MHD fluid kinetic energy, andMHD thermal energy are increasing. This clearly indicatesthat the instability is driven by the energetic ions. We canconclude that the instability is the fishbone instability. It canbe seen in Fig. 2(b) that the spatial profile of the fishbone
mode is different from that of the kink mode shown inFig. 2(a).
For the fishbone instability, amplitude evolution ofm/n = 1/1 harmonics of the radial velocity v r is shown inFig. 3(a). The frequency evolution of m/n = 1/1 harmonics
of the radial magnetic field fluctuation is shown in Fig. 3(b).
Frequency in the linear growth phase is 18 kHz whichcorresponds to ~ + 3.5 # 10 - 2 vA/R 0 and is close to the
fishbone frequency in the PDX experiments, where ~ is the
fishbone mode frequency, and vA is the Alfven velocity.

j'h / # d v* + vB n fd 3 v + d # M

1 P d # b - P d lnB # b - d # c Ph = bm, (14)=jh + B_ hih=B
M =- # n b f d 3 v ,(15)

where d # M is the energetic ion magnetized current,
Ph and Ph 9 are the energetic ion parallel pressure andperpendicular pressure, respectively. For the initialdistribution, the energetic ions are assumed to be isotropic in
the velocity space. Equation (13) has an exact simple form
for isotropic energetic ion distributions,t 2 v + t v $ dv = n1 ^d # Bh # B - dPh - dp.
2t0(16)The df particle simulation method [6] is applied to theenergetic ions. Time evolution of the j -th particle’s weight,
w j , is described byRVd ~ =- (1 - ~ ) S v + v d b $ d + e d v o 2 Wln f , (17)_ij S E
dt j
dt 2v WW 0S
1TXndv
qh= b + t d # bB $ < mh E F + d b $ =- mh dB G ,e(18)dt o 8
1d b = b - b0 ,(19)where f 0 is the initial distribution which is a function of themagnetic surface and energy. Using this weight, the energetic
ion current j lh in Eq. (14) is evaluated throughPh = Ph 0 + ! ~ j mh v 2 j S _ x - X j i,(20)Ph = = Ph = 0 + B ! ~ j n j S _ x - X j i,

(21)

jj

where S (x - X j ) is the shape factor of each super particle.

3. Simulation results
Results of particle-MHD hybrid simulations of the kinkinstability and the fishbone instability are reported in thissection. In the simulation, parameters are similar to the PDX
experiments. The major radius is R 0 = 1.43 m and the minor2619

Shiozaki Y. and Todo Y., Computer Simulation of Fishbone Oscillation

Fig. 2 Comparison in radial velocity spatial profiles of (a)kink mode for b h (0) = 0 % and (b) fishbone mode forb h (0) = 3 % .
Solid and dashed curves represent the sine andcosine part of the m/n = 1/1 harmonic, respectively.
Dashed-dotted curves show the q -profile.

Fig. 3 (a) Amplitude evolution of m/n = 1/1 harmonic of theradial velocity vr and (b) frequency evolution ofm/n = 1/1 harmonic of the radial magnetic fluctuationdBr .

The frequency shifts downward at saturation. The saturationlevel of the total velocity is v/v A + 2 # 10 - 3 . This
gives the plasma displacement at the saturation-2
p + v/~ + 5.7 # 10 R 0 + 0.18a + 0.6rq = 1 , where rq = 1
denotes the minor radius of the magnetic surface where thesafety factor is q = 1. Thus, the saturation level of the plasma
displacement is comparable to the radius of the q = 1 surface.
This saturation level is a factor of 7 larger than that foundin the simulation of the thermal ion diamagnetic drift
fishbone [4].

clear as Ref. [7] predicts. We should notice that Ref. [7]assumed that no energetic particles are present around theq = 1 surface. In the simulation presented in this paper,
energetic ions were distributed in a Gaussian form withregard to the minor radius and were present around the q = 1surface. The energetic particles around the q = 1 surface inthe present simulation might modify the two step structure.
Another simulation where the energetic particles are absent
from the q = 1 surface will be useful to test the conjecture.
In this paper we have investigated the preccessional
fishbone instability using parameters similar to the PDXexperiment. The spatial profile of the fishbone mode isdifferent from the top-hat structure of the kink mode. The
saturation of the fishbone instability and the frequency
downshift have been demonstrated. The saturation level ofthe plasma displacement was comparable to the radius of theq = 1 surface. More work is needed to clarify the mechanismof saturation and frequency downshift.

4. Discussion and summary
It was theoretically predicted in Ref. [7] that thefishbone mode has a two step structure and the width of theindividual steps is a factor c/~ smaller than the distance
between them, where c and ~ are the growth rate and thereal frequency of the mode, respectively. In the fishbone
simulation results they are c + 8 # 10 - 3 ~A and-2~ + 4 # 10 ~A , where ~A = vA/R 0 . The ratio is c/~ + 0.2 .
We can see a two step structure in Fig. 2(b), but it is not so3620

Shiozaki Y. and Todo Y., Computer Simulation of Fishbone Oscillation

References[4] J. Candy, H.L. Berk, B.N. Breizman and F. Porcelli,Phys. Plasmas 6, 1822 (1999).
[5] Y. Todo and T. Sato, Phys. Plasmas 5, 1321 (1998).
[6] S.E. Parker and W.W. Lee, Phys. Fluids B 5, 77 (1993).
[7] A. Ödblom, B.N. Breizman, S.E. Sharapov, T.C. Henderand V.P. Pastukhov, Phys. Plasmas 9, 155 (2002).

[1] K. Mcguire, R. Goldston, M. Bell et al., Phys. Rev. Lett.
50, 891 (1983).
[2] L. Chen, R.B. White and M.N. Rosenbluth, Phys. Rev.
Lett. 52, 1122 (1984).
[3] B. Coppi and F. Porcelli, Phys. Rev. Lett. 57, 2272
(1986).

4621