Development of Momentum Conserving Monte Carlo Simulation Code for ECCD Study in Helical Plasmas

Parallel momentum conserving collision model is developed for GNET code, in which a linearized drift kinetic equation is solved in the five dimensional phase-space to study the electron cyclotron current drive (ECCD) in helical plasmas. In order to conserve the parallel momentum, we introduce a field particle collision term in addition to the test particle collision term. Two types of the field particle collision term are considered. One is the high speed limit model, where the momentum conserving term does not depend on the velocity of the background plasma and can be expressed in a simple form. The other is the velocity dependent model, which is derived from the Fokker–Planck collision term directly. In the velocity dependent model the field particle operator can be expressed using Legendre polynominals and, introducing the Rosenbluth potential, we derive the field particle term for each Legendre polynominals. In the GNET code, we introduce an iterative process to implement the momentum conserving collision operator. The high speed limit model is applied to the ECCD simulation of the heliotron-J plasma. The simulation results show a good conservation of the momentum with the iterative scheme.


Introduction
Efficient current dive by the electron cyclotron heating (ECH) is an important issue in tokamaks and, also, helical devices.In helical plasmas, because of three dimensional magnetic configuration, behaviors of supra-thermal electrons are complicated and the radial transport of those electrons would alter the ECH physics.Thus it is necessary to study ECH physics including the complex electron orbits in helical plasmas.We have developed GNET code [1], which can solve a linearized drift kinetic equation in the five dimensional phase-space (two velocity and three real spaces).We have studied the radial transport of supra-thermal electrons by the ECH [1,2] and the electron cyclotron current drive (ECCD) [3] in helical plasmas.
However, in the present GNET code, a linear Monte Carlo collision operator is applied.This operator consider the collisional effect between test particle and background particle only as the pitch angle scatter and energy scattering.The change of the background particle distribution and the momentum transferred from the test particle to the background are ignored in this operator.The conservation of the momentum during particle collisions is an important issue in studying the ECCD, the neoclassical transport and etc.
It is pointed out that conserving momentum in electron-electron collisions may have the effect for the current [4], it is necessary to modify the operator to conserve momentum.The Ray tracing code shows a large impact of parallel momentum conservation for ECCD simulation [5,6].However, in the previous study the finite orbit and radial drift effects are not considered because of a radially local assumption.
In this paper, in order to study ECCD quantitatively, we develop the model conserving the parallel momentum for GNET.In order to conserve the parallel momentum, we introduce the field particle collision term in addition to the test particle one.Two models of the field particle collision term are considered.One is the high speed limit model, where the momentum conserving term does not depend on the velocity of the background plasma, and is easy to implement.The other is the velocity dependent model, which is derived from the Fokker-Planck collision term directly.
In the GNET code, we introduce an iterative process to implement the momentum conserving collision operator.We apply the developed model to the ECCD in the Heliotron J plasma with the momentum conserving operators.

Momentum conserving model
The parallel momentum conserving model is developed in order to implement to the GNET code, which is based on the Monte Carlo method.In the GNET code, we consider the change of the electron distribution from the Maxwellian distributuion, δ f , by the ECH as and we solve the linearized drift kinetic equation for δ f in the five dimensional phase-space as where C lin is the linear collision term, and S ql and L orbit are the quasi-linear ECH heating term and the orbit loss term, respectively.In order to consider the parallel momentum conserving model we assume the particle collision term as where is the linear collision term and C( f max , δ f ) is the field particle collision term which represents the collision effect for the background particles.
In this study we consider two momentum conserving models.One is the high speed limit model.In the high speed limit, the relative velocity between the test and background particles is given only by the test particle one.Then the collision frequency does not depends on the background particle velocity and the momentum conserving term does not depend on the velocity of the background plasma.The other is the velocity dependent model, which is derived from the Fokker-Planck collision term directly.

High speed limit model
In the high speed limit model, we assume that ECH heated supra-thermal electron velocity is much faster than thermal one and the field particle collision term C( f max , δ f ) can be expressed as where p(x, v) is a function of the real space coordinate, x, and the velocity, v. p(x, v) is determined by the conservation laws of energy and momentum.After some calculations, we obtain where n 0 means the density of background electrons.If we take a gyro-phase average of p(x) the perpendicular velocity component cancels out and we obtain Once p(x, v) is obtained from Eq. ( 5), we can calculate C( f max , δ f ) which compensates the lost momentum and energy from test particle.Then we can consider C( f max , δ f ) as a new source-sink term.
In the GNET code, if we iteratively calculate until δ f converges, we obtain a final profile of C( f max , δ f ).We label the steady state solution obtained by using S ql as δ f 0 and use C( f max , δ f 0 ) which becomes a new source term.The steady state solution of this source term is δ f 1 .Obtaining δ f 1 , we can consider the conservation of momentum when we calculate δ f 0 .However the test particle lost the momentum due to the collision with background plasma in the process of calculating δ f 1 .Therefore we iteratively calculate δ f n (n is the natural number) as until the lost momentum approaches 0. At the same time we evaluate the momentum which the test particle lost and calculate the momentum loss rate from them.We stop the iteration when the momentum loss rate becomes small sufficiently.Finally, we obtain the distribution function, δ f , with the momentum conserving collision operator by calculating

Velocity dependent model
Although the high speed limit model conserves the momentum, it does not include the exact information in the velocity space.We derive the more exact model, the velocity dependent model, from the Fokker-Planck collisional term directly.This model includes exact information more than the high speed limit one.The field particle operator can be expressed using Legendre polynominals P n (cos θ) as where v is the total velocity of an electron and θ represents the pitch angle.Introducing the Trubnikov-Rosenbluth potential and define u = cos θ to simplfy [7,8], we can de- Λ e/e represents the amplitude of field particle term and in this paper it is assumed as Λ c e 4 /m 2 e 2 0 , where Λ c is coulomb logarithm, e is charge, m e is mass of an electron and 0 is permittivity in vacuume.
In order to obtain the field particle term C n ( f max , δ f (n) (v)) which is determined by the obtained pertubed distribution function δ f (n) (v).We can iteratively calculate δ f (n) (v) in the same way with the simple model case.
After the iterative method, we calculate the complete collision operator according to Eq. ( 11).

Simulation result
In this study we apply the momentum conserving mode assuming the same magnetic configuration, heating and plasma parameters as the previous paper [3].We consider the Heliotron J configuration with the lower bumpy magnetic configuration b = 0.01 at the magnetic flux surface ρ = 0.67.The radial heating point is set to (r 0 /a, φ 0 , θ 0 ) = (0.1, 45 • , 0 • ).ECRF heating power is deposited at the top of the ripple in this configuration.We assume the quasi-linear ECH heating term as where  We set the parameters describing the EC resonance condition as follows: EC wave frequency is 70 GHz, 2ω ce /ω = 0.98, n = 0.44 and Δ = 1.0 × 10 −3 .
We run the GNET iteratively and obtain the steady state solution, δ f .Fig. 1 (a) shows the firstly obtained distribution function.The distribution becomes asymmetric in v at the high energy region.This is because many ECRH accelerated electrons hardly become trapped and the collisional relaxtion of the electron deficit in low energy region is faster than that of the accelerated electrons.As a result, the excess of electrons with positive v occured and it is found that the negative toroidal current is driven by the Fisch-Boozer effect.Fig. 1 (b) shows the sourcesink term to conserve the momentum using the steady state solution δ f 0 .Then the steady state solution δ f 1 is evaluated using the this source-sink term (Fig. 1 (c)) and again the next source-sink term is evaluated (Fig. 1 (d)).In the two source-sink terms we can see the larger distribution in the positive v region and this means the lost momentum have large effect in this region.
Fig. 2 shows the momentum loss rate at the each iterative calculation in the simple model.We evaluate the momentum loss at each calculation, and define the momentum loss rate as p loss = (p 0 − p n )/(p 0 ), where p 0 is the momentum lost by test particle at first calculation and p n represents one at the n th iterative caluculation.Fig. 2 shows the momentum loss decreases as the iterative calculation advanced and dropped less than 5% of lost momen-tum than initial simulation.The calculated ECCD current of the simple model is -20.1kA and the non-conserving one is -18.4kA.We can see the calculated ECCD current is larger by 9.2% than that of non-conserving one.

Conclusion
In order to study the ECCD physics on helical plasmas, we have developed the parallel momentum conserving collision model for GNET code, in which a linearized drift kinetic equation is solved in the five dimensional phasespace using Monte Carlo method.we have introduced the field particle collision term to conserve the parallel momentum in addition to the test particle collision term.Two types of the field particle collision term have been considered.One is the high speed limit model, where the momentum conserving term does not depend on the velocity of the background plasma and can be expressed in a simple form.The other is the velocity dependent model, which is derived from the Fokker-Planck collision term directly.Using the field particle collision term we apply an iterative method to implement the momentum conserving collision model.The high speed limit model has been applied to the ECCD simulation of the heliotron-J plasma.The simulation results has shown a good conservation of the momentum with the iterative scheme.

Figure 2 .
Figure 2. The momentum loss rate at each iterative calculation.