Ion finite Larmor radius effects on the interchange instability in an open system

A particle simulation of an interchange instability was performed by taking into account the ion finite Larmor radius (FLR) effects. It is found that the interchange instability with large FLR grows in two phases, that is, linearly growing phase and the nonlinear phase subsequent to the linear phase, where the instability grows exponentially in both phases. The linear growth rates observed in the simulation agree well with the theoretical calculation. The effects of FLR are usually taken in the fluid simulation through the gyroviscosity, the effects of which are verified in the particle simulation with large FLR regime. The gyroviscous cancellation phenomenon observed in the particle simulation causes the drifts in the direction of ion diamagnetic drifts.

Ion finite Larmor radius effects on the interchange instability in an open system I. Katanuma, a) S. Sato, Y. Okuyama, S. Kato, and R. Kubota Plasma Research Center, University of Tsukuba, Tsukuba, Ibaraki 305-8577, Japan (Received 10 September 2013; accepted 24 October 2013; published online 14 November 2013) A particle simulation of an interchange instability was performed by taking into account the ion finite Larmor radius (FLR) effects.It is found that the interchange instability with large FLR grows in two phases, that is, linearly growing phase and the nonlinear phase subsequent to the linear phase, where the instability grows exponentially in both phases.The linear growth rates observed in the simulation agree well with the theoretical calculation.The effects of FLR are usually taken in the fluid simulation through the gyroviscosity, the effects of which are verified in the particle simulation with large FLR regime.The gyroviscous cancellation phenomenon observed in the particle simulation causes the drifts in the direction of ion diamagnetic drifts.V C 2013 Author(s).All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License.[http://dx.doi.org/10.1063/1.4829682] The interchange instability is perhaps the most fundamental and basic magnetohydrodynamic (MHD) instability for magnetically confined plasma.There are some methods to stabilize the interchange instability.One is to make use of the line tying effects. 1,2The magnetic surfaces created in the torus, 3 for example, stabilize the interchange instability with the rotation of the magnetic field lines lying on the magnetic surface, which is the results of line tying effects.
The ion finite Larmor radius (FLR) effects are also expected to stabilize the interchange instability.The first work of the FLR effects on an interchange instability was seen in Ref. 4, where Rosenbluth, Krall, and Rostoker have found the stability effects of FLR with the help of Vlasov equation.Roberts and Taylor derived the same stability condition as Ref. 4 by using the extended MHD equations. 5ere, they used the generalized Ohm's law and they included the viscosity in the equation of motion.The extended MHD equations have been widely used to study interchange instability, because they are able to be applied to the complicated magnetic confined closed systems. 6,7[10][11][12] The generalized Ohm's law and gyroviscosity in the extended MHD equation were derived theoretically. 5,13ecently, it has been reported that the complete FLR stabilization of the interchange mode may not be attainable by the gyroviscosity or generalized Ohm's law alone in the frame work of extended MHD. 7The stability boundary of an interchange instability with FLR was determined by the kinetic analysis. 4,14Thus, it is worth verifying with the help of a particle simulation that the interchange instability can be really stabilized by FLR effects completely.
The traditional electrostatic particle-in-cell code (explicit PIC code) uses the equation of ion and electron motions with the Poisson equation. 15,16The mesh interval D and the time step Dt have to be taken smaller than the Debye length k De and the inverse of the plasma oscillation x À1 pe and electron cyclotron frequency X À1 ce , respectively, because all of the electrostatic oscillations are included in the explicit code.
The interchange instability by using the explicit PIC code was carried out by Goede, Humanic, and Dawson,17 where two-dimensional spatial 64 Â 64 grid was used.In order to follow such a slow time scale of interchange instability, they used the very small ion/electron mass ratio M i /M e ¼ 1.The linear growth rate of an interchange instability observed in the simulation did not compare well with the growth rate derived from the Vlasov equation in a local approximation, although it showed the stabilization due to FLR effects.Recently, the interchange instabilities mainly in the magnetotail are simulated by the electromagnetic explicit PIC code. 18,19However, the comparison of the interchange instabilities itself observed in the simulation with those derived theoretically has not been made.
The increases of the time step Dt and the grid interval D enable to simulate a low frequency phenomena in the plasma on a large scale.The implicit time integration scheme has a potential to increase Dt and D keeping the numerical stability of the simulation.1][22] Watanabe et al. have described the implicit algorithm and applied to the two-dimensional plasma with external magnetic field, where ion and electron cross-field motions are assumed only E Â B drifts. 23Barnes et al. have developed the implicit algorithm which can be able to applied to the two-dimensional plasma with the external magnetic field directly, and they demonstrated the simulation on an interchange instability (without FLR). 24his paper uses the uniform gravitational field g ¼ gê x shown in Fig. 1, where êx is the unit vector along x-axis.Here, the centrifugal force due to the non-zero magnetic field line curvature is replaced by the gravitational force.The uniform external magnetic field B ¼ Bê z is applied along z-axis.The FLR effects on the interchange instability are investigated in the geometry of Fig. 1.The PIC code used in this paper has adopted the implicit algorithm by Barnes et al. 24

II. LINEAR GROWTH RATES OF AN INTERCHANGE INSTABILITY WITH FLR
Rosenbluth and Simon derived the dispersion relation of an interchange instability with the FLR effects 14 where x is the wave frequency, k y is the wave number, E y is the perturbed electric field of the interchange instability, and c is the light speed.The steady state electric field E 0 , gravitational acceleration g and equilibrium density gradient dq=dx are assumed to be a function of x and their vectors are in the x direction.The uniform external magnetic field B is assumed to be applied in the z direction.The coefficients T and S in Eq. ( 1) are defined as where x pi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 4pn i e 2 =M i p is the ion plasma frequency, X ci eB=M i c is the ion cyclotron frequency, e is the unit charge, M i is the ion mass, and n i is the ion number density.The mass density q and pressure P of ions are given as Here, f i is the ion distribution function in the equilibrium state as The distribution in the equilibrium state should be a function of constant of motion, so that X is a canonical momentum in the y direction which is the same as the guiding center position in the case of uniform external magnetic field, i.e., X ¼ x þ v y =X ci , where x is the real position of ion in the x direction.Henceforth, the ion guiding center density profile n g (X) is assumed to be where h(x) is the Heaviside step function [h(x) ¼ 0 for x < 0 and h(x) ¼ 1 for x > 0].The guiding center density n g (X) of Eq. ( 5) lead to the real density n i (x) of where q i ffiffiffiffiffiffiffiffiffiffiffi ffi Figure 2 plots the linear growth rates of an interchange instability with FLR in the case of k y g=X 2 ci ¼ 10 À4 , which FIG. 1.Initial ion guiding center positions, gravitational acceleration vector (g), and external magnetic field.Vertical (horizontal) axis is the x axis (y axis) with the scale in mesh numbers.

FIG. 2. Dispersion relation of the interchange instabilities with FLR. (a)
is the linear growth rates c as a function of X ci =x pi , where each number in (a) denotes the magnitude of k y q i .(b) is the linear growth rates as a function k y q i , where each number in (b) denotes the magnitude of X ci =x pi .
was obtained by solving Eq. ( 1) in the range of 0 k y x 2p with the boundary condition that w ¼ 0 at k y x ¼ 0 and w ¼ 0 at k y x ¼ 2p and k y L H ¼ p.The vertical axis c in Fig. 2 is the linear growth rate of the interchange instability, i.e., x ¼ x r þ ic, where i ffiffiffiffiffiffi ffi À1 p .The linear dispersion relation in the geometry of Fig. 1 for k y q i ¼ 0 is given by 25 which is plotted in Fig. 2(a) in case of k y q i ¼ 0. Figure 2 indicates that the interchange instability is stabilized more strongly by FLR for the smaller X ci =x pi and larger k y q i .

III. ELECTROSTATIC PIC CODE
The electrostatic PIC code is used in order to research the stabilization effects of FLR on an interchange instability.The PIC code used in this paper adopts the implicit algorithm, 1,2,24 so as to remove the unnecessary high frequency electrostatic oscillations such as the electron cyclotron waves.
The equation of motion makes use of the modified leapfrog differential scheme where the superscript n means the time step nDt and the subscript a indicates the particle species.The electrostatic potential / n is solved by the following equations: where the second equation in Eq. ( 10) is the Poisson equation.
The velocity v nþc 0 j in the right hand side of Eq. ( 9) is defined by If c 0 ¼ 0 and / n instead of / n are used, Eq. ( 9) becomes the normal centered leap-frog scheme with the second order accuracy of time step Dt and is numerically unstable if jX ca jDt > 1.The detailed algorithm of the implicit PIC code is given in Ref. 24.
In this paper, c 0 ¼ 0.1 has been adopted for electrons and the time step is chosen x pe Dt ¼ 2 for the simulation of an interchange instability, where x pe ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 4pn e e 2 =M e p is the electron plasma frequency.However, c 0 ¼ 0 is chosen for ions so that the ion Larmor motion is followed with the second order accuracy of time step Dt.

IV. SIMULATION RESULTS
The code uses the two-dimensional 128 Â 128 spatial meshes in xy and three velocity components v x , v y , and v z .The (4 Â 128) 2 ¼ 262, 144 ions, and electrons each are included in the simulation.The geometry used in the simulation and analysis is plotted in Fig. 1.The ion guiding centers are uniformly distributed in the region x < 64D, where D is the spatial mesh interval, at t ¼ 0 and electron real positions are distributed at each ion real position.The uniform magnetic field B is applied in the z direction and gravitational acceleration g is in the x direction (g ¼ gê x ), as shown in Fig. 1.Henceforth, x pi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4phn i ie 2 =M i p and x pe ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4phn i ie 2 =M e p are used as ion and electron plasma frequencies, where hn i i ¼ n 0 =2 is the ion (electron) density averaged over the entire simulation space.
The interchange instability is observed in the simulation because the geometry in Fig. 1 is unstable to the interchange modes.Figure 3 plots the linear growth rates of interchange instability with (1, 1) mode measured in the linearly growing phase of the simulation.Here the (1, 1) mode means the one with the wave numbers k y ¼ 2p=128D and k x ¼ 2p=128D.The solid circles in Fig. 3 are the results without FLR, where c 0 ¼ 0.1 has been chosen for ions in Eqs. ( 9) and (11).The solid line arrowed by a symbol s A is the theoretically calculated linear growth rate Eq. ( 8).
The solid squares in Fig. 3 are the simulation results with FLR, where c 0 ¼ 0 is set for ions in Eqs. ( 9) and ( 11).The solid line arrowed by a symbol s B is the theoretically calculated linear growth rate in case of k y q i ¼ 0:16 in Fig. 2(a), where k y q i =ðk y g=X 2 ci Þ 1=4 ¼ 1:6.The simulation parameters are that k y hq i i ¼ 2:845 Â 10 À1 ; X ci =x pe ¼ 2:2155 Â 10 À3 ; k y g=X 2 ci ¼ 1:0 Â 10 À3 , where k y ¼ 2p=128D and Note that all quantities in Eqs. ( 1) and (2) are normalized as where x E k y cE 0 =B, although E 0 ¼ 0 is assumed throughout this paper.The solid squares in Fig. 3  changing the ion mass M i and the initial ion temperature T i in the simulation.The agreements between theory and simulation are good in both cases with and without FLR in Fig. 3.It is found that the flute instability is stabilized at X ci =x pi ' 0:5 in Fig. 3. Figure 4 plots the simulation results with parameters that M i =M e ¼ 1800; k y g=X 2 ci ¼ 9:94 Â 10 À4 ; X ci =x pi ¼ 0:943; X ce =x pe ¼ 40.The ion thermal Larmor radius hq i i in Fig. 4 is normalized by the spatial mesh interval D. The solid circles which are the linear growth rate observed in the simulation are obtained by changing initial ion temperature T i .The solid line in Fig. 4 is obtained by solving Eq. ( 1) with X ci =x pi ¼ 0:943.The agreements between the linear growth rates obtained analytically and by using simulation are good.There is the apparent systematic (but small) discrepancy observed in Fig. 4 at small values of the ion Larmor radius hq i i < 2. The linear growth rates (solid circles) in Fig. 4 have been measured from the wave form of field energy as shown in Fig. 5. Thus, the discrepancy is beyond the measuring error, but we do not know the reason.This systematic discrepancy can be seen in Fig. 3 at hq i i ¼ 0, which is consistent with those in Fig. 4.
The time evolutions of field energy j/ k j 2 of (1, 1) mode are plotted in Fig. 5.The dotted straight lines in the figure were used to obtain the linear growth rates (solid circles) in Fig. 4. Here, each origin of field energy has been shifted upward (or downward) at the amount of arbitral magnitudes in the figure, so that the relative magnitudes of j/ k j 2 with different hq i i are meaningless in Fig. 5.
Figure 6 plots the ion and electron real positions at x pet ¼ 11 000, which is the time when the interchange instability enters the nonlinear growing phase from the linear phase as is seen in Fig. 5.The ion real positions are plotted in the upper figures and the electrons are in the lower figures, respectively.It is found that the quasi-neutral condition is satisfied very much in all cases of hq i i.Here, Fig. 6(a) is the simulation result with hq i i ¼ 0:5, Fig. 6(b) is hq i i ¼ 1:5, Fig. 6(c) is hq i i ¼ 3:0, and Fig. 6(d) is hq i i ¼ 6:0, respectively.When the thermal ion Larmor radius hq i i, where normalized by a spatial mesh interval D, is small, there can be seen many interchange instabilities with high mode numbers, while in the case of large hq i i an interchange instability with only the lowest k y is seen.The interchange instability is almost stabilized by FLR for hq i i ¼ 6:0.The ions begin to behave like a viscous fluid in the presence of the magnetic viscosity in Figs.6(c) and 6(d).Equation (8), which is the linear growth rate in case q i ¼ 0, indicates that the linear growth rate is c / ffiffiffiffiffiffiffi jk y j p .In the simulation with hq i i 6 ¼ 0, ions and electrons have the initial Maxwellian velocities, so that the thermal fluctuations excite the interchange modes with various k y .In the case of hq i i Շ 2 the interchange mode with high k y are more unstable than that with low k y , so that many saturated high k y interchange instabilities with mushroom-shaped front 1,6 are observed in Figs.6(a) and 6(b).On the other hand, the high k y interchange modes are stabilized by the FLR effects in case of hq i i տ 2 and so the only interchange instability with the lowest k y is observed in Figs.6(c) and 6(d).
Figure 7 plots the time evolution of j/ k j 2 of (1, 1) mode in the case of k y g=X 2 ci ¼ 10 À3 ; X ci =x pi ¼ 1:0; hq i i ' 5:796.It is found that an interchange instability grows with different growth rates in the linear phase (solid straight line) and in the nonlinear phase (dotted broken straight line).Here the growth rate (solid straight line) in the linear growing phase has been calculated by the linear theory Eq. ( 1).The growth rates in the nonlinear phase are slower than those in the linear phase in Figs. 5 and 7.
As is seen in Fig. 8, in the linear phase (x pe t ¼ 4000 $ 10 000 in Fig. 7) the boundary between ions and vacuum has a sinusoidal shape.However, in the nonlinear phase, the shape of the boundary is distorted from the sinusoidal curve greatly.Here, Figs.8(a)-8(d) plot the equi-contour lines of ion density at x pe t ¼ 8000, 11 000, 14 000, and 17 000, respectively.Figures 8 and 7 are the results of the same simulation.As is seen in Fig. 7, the time 7000 Շ x pe t Շ 12 000 is the linearly growing phase of the interchange instability, so that the equi-contour lines (boundary between plasma and FIG. 4. Dispersion relation of the interchange instabilities in theory and simulation.The solid circles are the growth rates of the (1, 1) mode measured in the linearly growing phase of the interchange instability in the simulation with FLR.The solid line is the theoretical linear growth rates with FLR.FIG. 5. Time evolution of field energy j/ k j 2 with (1, 1) mode.The dashed straight lines in the figure are drawn in order to determine the linear growth rates in the simulation.The numbers in the figure denote the magnitudes of thermal ion Larmor radius hq i i, which are normalized by the mesh interval D. Here, each origin of field energy j/ k j 2 has been shifted upward (or downward) in order to display the growing phase clearly, so that the relative magnitudes of j/ k j 2 with different hq i i are meaningless in Fig. 5. vacuum) in Figs.8(a) and 8(b) can be approximate sinusoidal function well.On the other hand the time 12 000 Շx pe t Շ18 000 is the nonlinearly growing phase of the instability, so that the equi-contour lines (boundary between plasma and vacuum) in Figs.8(c) and 8(d) deviate from the sinusoidal function greatly.
The remarkable feature in Fig. 8(d) is that the front of the flute instability drifts in the -y direction, which is the ion diamagnetic direction.It is known that this drifts result from the gyroviscosity.The equation of motion in the extended MHD is written to be 5,7 The last term r Á P i is the viscosity term containing the twofluid effects and ion gyroviscosity effects.Frequently, the gyroviscous stress is approximated by 6,26 r Á P i ' Àqu Ã Á ru: Here, u Ã is often assumed to be the ion diamagnetic drift velocity in the uniform magnetic field, i.e., Moving the gyroviscosity term to the left hand in Eq. ( 13) in the approximation of Eq. ( 14) yields which is known as the "gyroviscous cancellation."The simulation has not solve the time evolution of the magnetic field because the electrostatic PIC code are used.That is, the uniform magnetic field remains until the end of the simulation.The equilibrium ion drift velocity, therefore, comes from the gravitational drift which is given by u g ¼ cM i g Â B=qB 2 ¼ À4:5 Â 10 À4 x pe Dê y in the parameters used to obtain Fig. 8. Here, Fig. 8 reveals that the interchange instability drifts in the -y direction.The drift speed u f ront of wave front of the interchange instability during 14 000 Շ x pe t Շ 17 000 in Fig. 8 It is found ju Ã j ) ju f ront j ) ju g j, although the wave front of the interchange instability drifts in the direction of the ion diamagnetic drift (that is, the gyroviscous cancellation phenomenon has been observed in this particle simulation with uniform magnetic field).

V. SUMMARY
The interchange instabilities in the geometry of Fig. 1 were investigated by using the electrostatic implicit PIC code.The growth rates of the interchange instability in the linearly growing phase in the particle simulation agree well with the theoretical linear growth rates with FLR which were obtained by solving Eq. (1).It is found that the interchange instability with the large FLR grows in two phases, that is, the linearly growing phase and the subsequent nonlinearly growing phase.The growth rate in the nonlinear phase is slower than that in the linear phase, and the interchange instability grows exponentially in both phases.The wave front of an interchange instability deviates from a sinusoidal shape in the nonlinear phase.The effects of gyroviscosity on the interchange instability seems to play an important role in its growth.The gyroviscous cancellation phenomenon has been observed in the particle simulation.The drift speed of wave front of the interchange instability, however, is much slower than the ion diamagnetic drift velocity.
The simulations were performed in the parameter of k y hq i i < 1, which has been assumed in all theoretical works with FLR. 4,5,8,9,11,12,14The particle simulation has introduced no approximations in the basic laws of mechanics and electricity.So the FLR is expected to stabilize the interchange instability completely in the real magnetically confined plasma.In the future, the particle simulations with k y hq i iտ1 will be performed to research the FLR effects on an interchange instability in the range X ci =x pi ) 1.

FIG.
FIG.Contour plots of ion number density.Here, the figure is the simulation result in the case k y g=X 2 ci ¼ 10 À3 ; X ci =x pi ¼ 1:0; hq i i ' 5:796.(a) is the contour plot of the ion density at x pe t ¼ 8000, (b) is at x pe t ¼ 11 000, (c) is x pe t ¼ 14 000, and (d) is x pe t ¼ 17 000, respectively.
The time evolution of field energy j/ k j 2 with (1, 1) mode.The solid straight line denotes the linear growth rate which was calculated analytically by using Eq.(1).The dashed straight line is drawn to fit the growth rate to the wave form in the nonlinear growing phase of j/ k j 2 .Here, the figure is the simulation result in the case k y g=X 2 ci ¼ 10 À3 ; X ci =x pi ¼ 1:0; hq i i ' 5:796.