Study of Self-Consistent Particle Flows in a Plasma Blob with Particle-in-Cell Simulations

The self-consistent particle (cid:13)ows in a (cid:12)lamentary coherent structure along the magnetic (cid:12)eld line in scrape-oﬀ layer (SOL) plasma (plasma blob) have been investigated by means of a three-dimensional electrostatic particle-in-cell (PIC) simulation code. The presence of the spiral current system composed of the diamagnetic and parallel currents in a blob is con(cid:12)rmed by the particle simulation without any assumed sheath boundary models. Furthermore, the observation of the electron and ion parallel velocity distributions in a blob shows that those distributions are far from Maxwellian due to modi(cid:12)cation with the sheath formation and that the electron temperature on the higher potential side in a blob is higher than that on the lower potential side. Also, it is found that the ions on the higher potential side are accelerated more intensively along the magnetic (cid:12)eld line than those on the lower potential side near the edge. This study indicates that particle simulations are able to provide an exact current closure to analysis of blob dynamics and will bring more accurate prediction of plasma transport in the SOL without any empirical assumptions.


I. INTRODUCTION
Early experimental studies of scrape-off layer (SOL) plasmas in magnetic confinement fusion devices implied the existence of fast convective plasma transport to the main chamber wall. 1 Recent experimental investigations of boundary layer plasmas have observed intermittent spikes of the plasma density fluctuation and intermittent filamentary coherent structures along the magnetic field lines. [2][3][4][5] Such observations are thought to be evidence of the convective transport across the magnetic field lines. The propagation mechanism of such coherent structures, which are called "blobs," was suggested 6 and many theoretical, numerical, and experimental investigations of the blob dynamics have been performed. [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21] However, most theoretical and numerical studies are based on two-dimensional reduced fluid models and closures for the parallel currents and kinetic effects (such as sheath formation between a SOL plasma and a divertor plate, velocity difference between ions and electrons, and other issues) are treated under some assumed and parameterized models in such fluid limits. Thus, we have developed a three-dimensional electrostatic particle-in-cell (PIC) code with particle absorbing boundaries in order to study blob dynamics. 22 Preliminary PIC simulations have shown that blobs propagate to the first wall across the magnetic field lines and that a blob evolves into a mushroom-shaped object. 23 In this study, we have investigated self-consistent particle flows in a blob by means of a three-dimensional (three space coordinates and three velocity components) electrostatic (electric field is solved and magnetic field is constant in time) PIC simulation code. In the code, full plasma particle (electron and ion) dynamics including the Larmor gyration motion for all particles in a blob and background plasma and the self-consistent electric field formed by the charge density which is obtained from all particles are solved with the equation of motion and Poisson's equation, respectively. 24 In Sec. II, we briefly describe the configuration and parameters of simulation. Also we compare results of preliminary simulations with the fluid theory. In Sec. III, we present results of the particle simulation. The results show that the self-consistent and self-organized structures of potential, flows, and current are formed in a blob. We then investigate the properties of particle velocity distributions in a blob. Finally, a summary of our work is given in Sec. IV.

A. Configuration and parameters
We study self-consistent particle flows in a blob by a three-dimensional electrostatic particle simulation code with full electron and ion dynamics. A slab geometry with x the counter radial direction, y the poloidal direction, and z the toroidal direction (parallel to the magnetic field B) is applied to the code, as shown in Fig. 1. The external magnetic field strength B increases in the positive x direction as where L x , L y , and L z are the system size in the x, y, and z directions and B Lx is the magnetic field strength at x = L x . At z = 0, L z (corresponding to the divertor plates), and x = 0 (corresponding to the first wall) (i.e., shaded plates in Fig. 1), the particle absorbing boundary where the electric potential ϕ is set as ϕ = 0 is placed. On the other hand, periodic boundary condition is applied in the y direction. At x = L x , the reflecting boundary condition, in which ∂ϕ/∂x = 0, is used. Although the system has no poloidal flows at the core side boundary under this boundary condition, we apply such a simple condition in order to study fundamental mechanisms on blob dynamics.
The initial plasma density distribution which includes the blob and background plasma is given by where n 0 is the density of the background plasma, n b0 is the initial density amplitude of the blob, and δ bx and δ by are the blob sizes in the x and y directions. Equation (2) means that the blob is initially located as a column along the ambient magnetic field at around (x, y) = (x b0 , y b0 ). The initial temperature in the blob is equal to that of the background plasma. Though the initial velocity distribution is given by Maxwellian, the velocity distribution after the start of calculation will be self-consistently modified according to sheath formation. This system has no particle and heat sources. Although the initial density distribution (also density at any time) is not in equilibrium with the magnetic field, this assumption is appropriate in the low beta limit.
The common simulation parameters are as follows. The number of spatial cells is v T i /v T e = 0.05. Thus, the initial ion-to-electron temperature ratio is T i /T e = 0.25. The external magnetic field strength is given as Ω i /ω pi = 1 where ω pi is the ion plasma frequency in the background plasma. The initial density ratio of the blob to the background plasmas is n b0 /n 0 = 2.7. The initial blob size in the x direction is δ bx = 4 ∆ g . The initial position of the blob is (x b0 , y b0 ) = (3L x /4, L y /2). Using the above configuration and parameters, we find that this blob is in the sheath-interchange (C s ) regime shown in Ref. 15 because collisionless plasma is considered.

B. Comparison with fluid theory
We have performed some preliminary simulations in order to compare particle simulation results with the fluid theory of blob dynamics. The simulation parameters are shown in Table I. Here, ∆t found in this table is the time step. Figure 2 shows the time variations of  x nec on the x-y plane at z = L z /2. Here, x nec is the x component of position of the electron center of mass, which is obtained in the area where the electron density, n e , is higher than n 0 + n b0 /10. That is, x nec represents the position of a blob. Also, we show the relation between the poloidal blob size, δ by , and the blob propagation speed in the −x direction, v b , in Fig. 3. Here, v b is calculated from x nec from Ω i t = 80 to 140. The solid line in Fig. 3 represents the theoretical radial blob propagation speed given by (3), we assumed that the potential term and the poloidal derivative are approximated as Figure 3 indicates that the observed radial propagation speeds are in good agreement with the fluid model.

III. SIMULATION RESULTS
In this section, we now show simulation results where the number of spatial cells in the system is set as N x ×N y ×N z = 64×64×512. Although we have given a small value to N z of the preliminary simulations shown in Sec. II B to save the computation time and computer resources, that of the simulation described in this section has a larger value in order to enable a long calculation because particles escape to the end plate along the magnetic field with a velocity of ∼ c s , i.e., the time of escape is estimated at ∼ L z /2c s . Other parameters (∆ g , ∆t, δ bx , and δ by ) are the same as those in Run 2 in Sec. II B. Thus, the system size is  Figure 4 presents the electron density distributions on the x-y plane at z = L z /2 at various times. From Fig. 4, we find that the blob propagates to the first wall and that the blob evolve to steepened shape and lower amplitude. The reason why the amplitude of blob is decreased is that plasma particles escape to the end plates.

A. Blob propagation
In Fig. 5, we show the time variation of position of the electron center of mass. Figure 5 indicates that the blob propagates in the radial direction with v b = 0.0551 c s which is calculated from x nec from Ω i t = 80 to 250 and that the poloidal velocity of the blob is quite small. its absolute value is much smaller than that of the diamagnetic current. In the fluid theory, the diamagnetic current and the current by the grad-B drift are given as and and In the case of this simulation, j gymax /j dymax ≈ 0.09. However, the net poloidal current contributed by j gy along the radial direction is larger than that of j dy because of the integrations ∫ x 1 x 0 j dy (y = y b0 ) dx = e n 0 c s ρ s and where we assumed that n is given by Eq. (2) and that x 0 ≪ x b0 ≪ x 1 . On the other hand, the polarization current is approximated as Therefore, the polarization current is also much smaller than the diamagnetic current.

C. 3D flows and current structures
In order to understand three-dimensional structures of flows, we show three-dimensional streamlines of the electric current, j, and the electron and ion fluxes, F e and F i , in Figs. 7-9.
Here, the flux and the current are defined by and respectively, where the subscript s refers to electrons (e) or ions (i), v is the velocity, f is the distribution function in the six-dimensional phase space, and q is the charge. In the simulation, F s is calculated by where x α , y β , and z γ are the components of position of grid whose number is (α, β, γ), x s,i , and i, respectively, and S is the form-factor of the finite-size particle. 24 Thus, F s includes almost all particle motions and drifts.
The electric current in the blob draws trajectories like two-dimensional vortices, as shown in Fig. 7. Such a configuration is caused by the composition of the parallel current and the poloidal current dominantly due to the diamagnetic current. The angle between the plane on which the vortex exists and the x-y plane becomes larger as z closes to the edge. component to other components of the ion flux is larger than that of the electron flux near the edge. This is because the electron diamagnetic drift speed is faster than that of the ion (i.e., the electron temperature is larger than the ion temperature). Furthermore, end points of streamlines of the electron and ion fluxes tend to concentrate in the lower and higher potential sides, respectively. This difference causes the inclination of the current vortex mentioned above.
To confirm the presence of the net poloidal current in the blob, which includes the diamagnetic, grad-B drift, and polarization currents, we calculate the net poloidal currents I y0 and I y1 from the simulation result, where and that is, (πδ bx δ by /2) I y1 (z) means the total of the poloidal current j y (x, y nemax , z ′ ) which pass through the area x l ≦ x ≦ L x , L z /2 ≦ z ′ ≦ z (see the schematic diagram in Fig. 10 (a)). Also, (πδ bx δ by /2) I y0 (z) is considered as the total of the poloidal current on the background plasma. Here, x l = L x /4. In Fig. 10 (b), we show I y1 (z) and I y0 (z) at t = 87.14 Ω −1 i . From the curve of I y1 (z) which is represented by the thick line, it is found that the net poloidal current flows from the lower to the higher potential sides. This net current is thought to be brought mainly by the grad-B drift current because of the discussion in Sec. III B.
Then, we compute the net toroidal current in the lower potential side, i.e., (πδ bx δ by /2) I z (z) is the total of the toroidal current j z (x, y, z) that penetrates the area x l ≦ x ≦ L x , y nemax − L y /2 ≦ y ≦ y nemax , as shown in the schematic diagram in Fig. 10 (a).
Thus, I z (z) is considered as the averaged toroidal current density on the lower potential side in the blob. The curve of I z (z) is also presented as the thin line in Fig. 10 (c). The thick line in Fig. 10 (c) represents I y1 (z) − I y0 (z). The panel (c) in Fig. 10 and the current continuity equation indicate that the toroidal current from the divertor plates is almost converted to the poloidal current. Also, the magnitude of the averaged toroidal current density is ∼ en 0 c s near the edge.

D. Kinetic properties in a blob
Although we have reduced kinetic particle data to the macroscopic (fluid) quantities and mentioned their properties regarding flows on the blob in the above discussion, we should see some kinetic properties in order to investigate kinetic effects on blob dynamics.
Then, we calculate the parallel velocity distributions on various parts in the blob. The thick line in Fig. 11 (a) represents the electron parallel velocity distribution on the higher potential side around z = L z /2 (averaged in the partial volume: 9/16 ≦ x/L x ≦ 13/16, 1/2 ≦ y/L y ≦ 11/16, 63/128 ≦ z/L z ≦ 65/128). On the other hand, the thin line in Fig. 11 (a) shows that on the lower potential side around z = L z /2 (averaged in the partial volume: 9/16 ≦ x/L x ≦ 13/16, 5/16 ≦ y/L y ≦ 1/2, 63/128 ≦ z/L z ≦ 65/128). Although both lines in Fig. 11 (a) indicate the loss of the fast particles, it is found that the electron temperature on the higher potential side is higher than the lower potential side. We obtain the effective parallel electron temperature on the higher and lower potential sides as T ez+ = 0.70 T e and T ez− = 0.58 T e , respectively, from these distributions. Since the electrons on the higher potential side are more strongly prevented from running away to the divertor plate by the electric potential, this temperature difference arises. Figure 11 (b) shows the ion parallel velocity distributions on the higher potential side (the thick line) and the lower potential side (the thin line) around z = L z [averaged in the partial volumes: (9/16 ≦ x/L x ≦ 13/16, 1/2 ≦ y/L y ≦ 11/16, 63/64 ≦ z/L z ≦ 1) and (9/16 ≦ x/L x ≦ 13/16, 5/16 ≦ y/L y ≦ 1/2, 63/64 ≦ z/L z ≦ 1), respectively]. Figure 11 (b) indicates that the ions on the higher potential side are accelerated more intensively than those on the lower potential side. From these distributions, we get the averaged ion parallel velocities on the higher and lower potential sides as ⟨v iz+ ⟩ = 1.06 c s and ⟨v iz− ⟩ = 0.97 c s , respectively. Because the sheath potential difference on the higher potential side is larger than that on the lower potential side, this property appears. Furthermore, we obtain the effective parallel ion temperatures on each side as T iz+ = 0.28 T i and T iz− = 0.25 T i .

IV. SUMMARY AND DISCUSSION
In this paper, we have studied self-consistent particle flows in a plasma blob with a first principle 3D calculation including particle parallel dynamics and Larmor motions. The simulation has self-consistently demonstrated the formation of the potential structure, the particle flows, the electric current, and the temperature structure in a blob without any artificial sheath boundary models. In most of the previous works about blob dynamics, the current closure in a blob is provided from simple and ideal models. 20 However, this study shows that particle simulations have the possibility of giving an exact current closure to the investigation of blob dynamics. Furthermore, the temperature structure and the velocity difference between the higher and lower potential sides (i.e., velocity shear) have the probabilities of inducing blob spin 13 and some kind of instability, 25,26 respectively. These problems are topics for future research. Though we show the obvious self-consistent structures in a blob in this paper, we will investigate temporal dynamics of coherent structures (e.g., some instabilities) in the following works. Although the code does not include any additional collision models (e.g., additional Coulomb collision, 27,28 collision with neutrals, 29 etc.) and the smaller mass ratio m i /m e , the smaller connection length L z , the smaller temperature ratio T i /T e , the larger grad-B than those in typical SOL, and the identical temperatures between the blob and background plasma are applied in the present work in order to investigate fundamental kinetic dynamics caused by particle motion, we plan to study the dependence on those parameters and perform more specific simulations where the parameters of real devices, additional collision models, and detachment dynamics will be employed. The direct comparison between the simulation results and experiments (including consideration of methods to compare) is also an important topic in future work.