Charge Separation and Electric Field at a Cylindrical Plasma Edge

1. Introduction

There are number of physical situations where plasmas neutrality breaks down through boundary layers called plasma sheaths, which are either free or in contact with a wall. The plasma sheaths transition problems are at the heart of an industrial revolution whose theme is the design of matter on the molecular scale.The study of the charge separation at a plasma edge requires generally the solution of the kinetic equations of plasmas which, for a collisionless plasma, usually reduce to the well-known Vlasov equation. Some examples for the solution of the Vlasov equation for sheaths transition problems have been presented in Shoucri, 2008a, 2009a. A problem of interest is the problem involving the generation of radial electric fields and poloidal flows to achieve radial force balance at a steep density gradient in the presence of an external magnetic field. This problem is of great importance in the steep density gradients pedestal of the high confinement mode (H-mode) in tokamaks, since it largely affects the edge physics of the H-mode. In the present work, we shall study the problem of the generation of a charge separation and the associated electric field at the edge of a cylindrical plasma column, in the presence of an external magnetic field directed along the cylinder axis. In previous publications on this problem (Shoucri et al., 2003, 2004, 2008b, 2009b), we have considered the case where the electrons were frozen by the magnetic field lines, with a constant density profile which changes rapidly along the gradient over an ion orbit size. Along the gradient the electrons bound by the magnetic field cannot move across this field to exactly compensate the ion charge which results from the finite ions’ gyroradius. This effect is especially important for large values of the ratio$$, where $$ is the ions’ gyroradius and $$is the Debye length. Accurate calculation of the charge separation is important for the accurate calculation of the self-consistent electric field. This requires also an accurate calculation of the exact ion orbits using a kinetic equation. In the present work, we use an Eulerian Vlasov code to study the evolution of the ion distribution function for the problem of the charge separation at the edge of a cylindrical plasma, in the presence of an external magnetic field directed along the cylinder axis. Eulerian codes have the advantage of a very low noise level and make it possible to measure accurately a very small charge separation (Shoucri et al., 2003, Shoucri et al., 2004 & Shoucri, 2009a,b), and allows accurate results in the low density regions of the phase-space.

It was pointed out in the analysis of H-mode power threshold in a tokamak by Groebner et al., 2001, that the changes in the electron density n_e and n_e in the transition to H-mode are small, and changes in T_e are barely perceptible in the data. Electrons and ions have a density profile which varies rapidly along the gradient over an ion orbit size. In previous publications (Shoucri et al., 2003, 2004, Shoucri, 2002, 2008a,b, 2009a,b) the electrons density and temperature profiles were kept constant. In the present work we will allow the electrons to move, and it will be sufficient for the purpose of our study to describe the motion of the electrons, having a small gyroradius, by a guiding center equation (Shoucri et al., 1997). This allows a more accurate description of the contribution of the electrons, with respect to the approximation previously used in Shoucri et al., 2003, 2004, 2009a, where the profile of the electrons was assumed constant in time. The electrons motion across the magnetic field in the gradient region is limited by the guiding center equation, and the magnetized electrons cannot move sufficiently across the magnetic field to compensate the charge separation which results from the ions motion due to the finite ion orbits. To determine this charge separation at the plasma edge along the gradient, it is important to calculate the ion orbits accurately by using an Eulerian Vlasov code. The larger the ions’ gyroradius, the bigger the charge separation and the self-consistent electric field at the edge. (Hence the important role played by even small fractions of impurity ions, because of their large gyroradii, in enhancing the electric field at the plasma edge).

In full toroidal geometry, there are “neoclassical” effects which can play a role in this problem, such as the neoclassical enhancement of the classical ion polarization drift, or the neoclassical damping of poloidal flows (Stix, 1973, Hirshman 1978, Waltz et. al., 1999). We focus for simplicity in the present work on a cylindrical geometry for the problem of the generation of an electric field and poloidal flow at a plasma edge due to the finite ions’ gyroradius, when the external magnetic field is applied along the axis of the cylinder. The solution we present is a two-dimensional (2D) solution in a cylindrical geometry, with the external magnetic field assumed uniform along its axis. We compare the radial electric field calculated along the gradient with the macroscopic values calculated from the kinetic code for the gradient of the ion pressure in the radial direction, and we find that this quantity balances the radial electric field fairly well, a result similar to what has been presented in one-dimensional geometry (Shoucri, 2002, Shoucri et al. 2003, 2004). The contribution of the Lorentz force term along the gradient is negligible. For the parameters we use in the present work, the solution allows for a small value of $$ to exist, especially on the high field side of the cylinder as it will be discussed below, which causes a small oscillation in the radial direction.

As mentioned above, the ions are described by a kinetic Vlasov equation and the electrons are described by a guiding center equation. These equations are solved numerically using a method of characteristics. There have been important advances in the last few decades in the domain of the numerical solution of hyperbolic type partial differential equations using the method of characteristics. The application of the method of characteristics for the numerical solution of the kinetic equations of plasmas and for the guiding center equations has been recently discussed in several publication (see for instance Shoucri, 2008a,c, 2009a). These methods are Eulerian methods which use a computational mesh to discretize the equations on a fixed grid, and have been successfully applied to different important problems in plasma physics involving kinetic equations, such as laser-plasma interaction (Ghizzo et al., 1990, 1992, Strozzi et al., 2006, Shoucri, 2008d), the calculation of an electric field at a plasma edge (Shoucri, 2002, 2008b, 2009b, Shoucri et al., 2000, 2003, 2004), the applications of gyro-kinetic codes to study edge physics problems in plasmas (Manfredi et al., 1996, Shoucri, 2001, Shoucri et al., 2005, Pohn & Shoucri, 2008) and to collisional plasmas (Batishchev et al., 1999). These methods present the great advantage of having a low noise level, and allow accurate results in the low density regions of the phase-space (Shoucri, 2008c). In the applications presented (Shoucri, 2008a,c, 2009a), the computation was usually done on a fixed grid, so no dynamical grid adjustment was necessary, and interpolation was restricted to the use of a cubic spline, which compared favourably with other methods (see, for instance, Pohn et al., 2005), so altogether the method was accurate and remained relatively simple. Interpolation in several dimensions using a tensor product of cubic B-splines has been also successfully applied (Sonnendrücker et al., 1999, Shoucri, 2008a, 2009a, 2011). The method of characteristics has been also successfully applied in fluids to problems having shock wave solution (Shoucri & Shoucri, 2007). Large Courant-Frederichs-Levy (CFL) computation parameter is possible, and therefore the time-step numerical limitation by large velocities can be removed, if the physics makes it possible. A more complete study on splines can be found in the book of Ahlberg et al., 1967, and an important theoretical study on the method of characteristics can be found, for instance, in the book of Abbott, 1966. More applications to plasma physics problems can be found in the several references we have cited.

2. The relevant equations and the numerical method

We consider the cylindrical coordinate system ($$. The plasma is uniform in the z direction. The radial direction in the cylindrical plasma is normal to a vessel surface which is located at r = r_max=R. The constant magnetic field is in the z direction which represents the toroidal direction, and is the poloidal direction. The constant magnetic field in the z direction is given by:

$$E1

where B₀ is the magnetic field along the axis of the cylinder at = /2, and $$ where $$ is the major radius of the tokamak. We assume that at r = r_max=R we have $$=0.2. In this case we can also write$$. We consider a deuterium plasma,$$. The ions are described by the normalized 2D Vlasov equation for the ions distribution function $$

$$E2

where from Eq.(1) the ion cyclotron frequency$$. The electrons are described by the normalized 2D guiding-center equation:

$$E3

where the drift velocity$$. Equation (3) can be developed in our 2D system to give the following equation:

$$E4

where $$ ;

$$E5

This system is coupled to Poisson’s equation:

$$$$

$$E6

Time is normalized to the inverse plasma frequency$$, velocity is normalized to the acoustic velocity$$, and length is normalized to the Debye length$$. T_e is the electron temperature and m_i is the ion mass. The potential is normalized to$$, and the electric field is normalized to$$, and the density is normalized to the peak initial central density. The ion cyclotron frequency $$as previously defined is normalized to$$. The system is solved over a length $$= 175 Debye lengths in front of the vessel surface, with an initial ion distribution function for the deuterons over the domain $$ given in our normalized units, by:

$$E7

$$and $$ are the initial ion and electron density profiles respectively, with $$ in the initially neutral system. R is the radius of the cylinder as we previously mentioned. We also use the following parameters:

$$E8

In the present calculations the parameter$$. With the initial distribution function for the ions in Eq. (6), we have $$ spatially constant, the factor 2T_i in Eq. (7) in the calculation of the gyroradius takes into account the fact that the perpendicular temperature is$$.

Equations (2) and (3) are solved by a method of fractional step (Yanenko, 1971), first applied for a Vlasov equation by Cheng & Knorr, 1976, and Gagné & Shoucri, 1977, coupled to a method of characteristics. For the general case where several dimensions are involved, the fractional step technique allows the reduction of the multi-dimensional equation to an equivalent set of reduced equations. To advance Eq. (2) and Eq.(3) for one time-step t, the splitting of the equations is applied as follows.

1. We solve for t/2 the equations:

$$E9

$$E10

Equations (8) and (9) are solved by interpolation along their characteristics, to be described below. We then solve Poisson’s equation in Eq.(5) to calculate the new electric field.

1. We solve for t the equations:

$$E11

$$E12

Again Eq.(10) is solved by a method of characteristics, to be described below.

1. We repeat Step1 for t/2, and then solve Eq.(5) to calculate the new electric field$$,$$.

This completes a one time-step cycle t.

2.1. The solution for Step 1

For the solution of Eq. (8), we first solve the characteristic equations:

$$E13

(12)

at a given$$ in velocity space. The solution of Eqs (12) originating at ($$) at time t, and reaching the grid point ($$) at $$ can be written as follows, for a half time-step t/2:

$$

$$E14

For v_r t/2 « 1, the second equation reduces to$$. Therefore the solution of Eq. (8) can be written, for half a time-step t/2, as follows:

$$E15

The right hand side of Eq. (14) is calculated by interpolation using a tensor product of cubic B-splines, in which is periodic (Shoucri et al., 2004, Shoucri, 2008a, 2009a). For each $$ and$$, we write for the interpolation function$$:

$$E16

taking into account that is periodic. $$and$$are cubic splines.The cubic B-spline is defined as:

$$E17

and $$ otherwise. The grid size h (either$$or $$in our notation), is assumed uniform. For the calculation of the coefficients $$ of the cubic B-spline interpolation function $$ for periodic $$ in Eq.(15) see details in Shoucri et al., 2004 and Shoucri, 2008a, 2009a.

We then solve Eq.(9) for a half time-step t/2 along the characteristics:

$$E18

Equations (17) are solved using an iterative process. We assume that at the time-step$$, r is at the grid point $$ and $$is at the grid point$$. The following leapfrog scheme can be written for the solution of Eq.(17):

$$E19

$$E20

where $$ is the point where the characteristic is originating at $$(not necessarily a grid point). Put:

$$

$$E21

we can rewrite Eqs.(18,19) as follows:

$$E22

$$E23

which are implicit equations for$$, and which are solved by iteration. Usually, two iterations are sufficient for convergence. We start with$$. We get at the first iteration:

$$

$$E24

$$

$$E25

Second iteration:

$$E26

$$

$$E27

The quantities $$ and $$are calculated by a cubic B-spline interpolation, similar to what is described in Eq.(15). Finally to advance Eq.(9) by$$, we calculate:

$$E28

where the right hand side of Eq.(27) is calculated again using a cubic B-spline interpolation as indicated in Eq.(15).

2.2. The solution for Step 2

We go to Step2 and solve for $$ Eqs.(10,11). We first solve the equation:

$$E29

The electric field $$remained small, since$$, and r ~ R =5000 towards the edge in our calculation.

$$E30

Following the same iterative steps as described for Eq.(27), to an order $$the solution of Eq. (29) yields the following solution to Eq.(28):

$$E31

Again the 2D interpolation in Eq. (30) is done using a tensor product of cubic B-spline, as explained in Eq.(15), and a and b are calculated with similar iterative steps as in Eqs.(23-26) and are given, to an order$$, by the expressions:

$$E32

$$E33

where$$;

$$E34

The density$$ is calculated by integrating Eq.(11) as follows:

$$E35

where $$ is calculated in Eq.(27). To calculate $$and$$, we then repeat Step 1 for the solution of Eqs. (8,9) for $$ using $$ and $$ calculated from Eq.(30) and Eq.(33). The electric field is then updated to calculate $$ and $$ by solving Poisson equation.

2.3. The solution of Poisson’s equation

Equation (5) is solved by first Fourier transforming the equation in the periodic direction$$:

$$E36

and the resulting equation is then discretized in the radial direction over a uniform grid following the method of Knorr et al., 1980 (also rediscussed more recently in Crouseilles et al. 2009):

$$E37

where

$$E38

and $$ ;$$;

$$E39

Equation (35) is solved using a tridiagonal algorithm. To get the boundary conditions, we assume in the present calculation that the deuterons and electrons currents hitting the cylindrical wall surface at$$are collected by a floating potential cylindrical vessel. Therefore we have the relations for the charge collected:

$$E40

where

$$E41

and $$ is defined in Eq.(17). $$from Eq.(37) is used to obtain for the potential the following boundary condition at$$:

$$E42

By integrating the relation$$, we get a relation between the charge collected on the cylindrical wall and the charge appearing in the initially neutral plasma:

$$E43

where $$ is the width of the plasma slab, which extends from $$to $$at the edge of the cylindrical plasma of radius$$. $$is the charge appearing in the system. From Eq.(40), we get the following relation:

$$E44

We assume that the plasma particles are allowed to enter or leave at the boundary at$$. So the difference between the electric fields at the boundary $$ and at r = R must be such that Eq. (40) must be satisfied at every time-step in every direction$$. In the present simulation, $$is calculated from Eq. (37), which defines the derivative of the potential at the right boundary to be used in the solution of Eq. (35). We fix the potential to be zero at the boundary at$$ and solve Poisson equation in Eq.(35) for the potential. Then the resulting electric field at$$, calculated at $$ from Eq. (5), must satisfy Eq. (41).

The initial density profiles at the neutral plasma edge are given by:

$$E45

The profiles in Eq. (42) situate the steep gradient to be centered at a distance of $$from the wall of the vessel, which put the plasma relatively close to the floating wall of the vessel. The wall of the vessel will collect the charge coming from the plasma, especially due to the large ions gyroradius. The system is solved for an edge thickness$$, and R=5000 for the radius of the cylinder (we expect this radius to be large in a tokamak, so in the domain$$, the variation of the quantity $$remains very close to 1.

We use N=250 grid points in space in the radial direction, and 128 grid points in the azimuthal direction. 80 grid points are used in each velocity direction. The velocity extrema for the ions velocities are $$ in our normalized units, with T_e/T_i =1 in the present simulations.

3. Results

The code was executed for a sufficiently long time to reach a steady state. We noted around t=500 that the code has indeed reach a steady state. Since $$ (normalized to$$), then the gyroperiod$$. This means that the code has executed more than 5 gyroperiods. Figure(1) shows the electric field at = 0 (full curve), = /2 (broken curve) and = (dashed-dotted curve), at time t = 495 (left figure) and t=500 (right figure). At the edge and along the gradient the electric field $$is directed towards the center to the interior of the plasma. The electric field at the floating vessel wall at r = R was calculated using Eq. (37) at 128 points over the 2 circle. We see from Fig. (1) that the electric field at = 0 is higher (in absolute value) than the one at =, and we note that the curves for 0<R-r<75 have reached a steady state stable equilibrium, with the charge collected at r=R remaining constant, while towards the center for 75<R-r<175, the curve at = shows a small steady state oscillation around zero. The charge collected on the floating wall of the cylindrical vessel at = 0, /2 and are respectively -0.1455, -0.1389 and -0.063 at t=500. For the solution of Poisson equation, we set the value of the electric field$$ at the cylindrical vessel wall exactly equal to the charge collected on the wall according to Eq.(37). If the charge collected on the wall is equal to the charge appearing in the system, then according to Eq.(41) :

$$E46

which is what we see for the curves at =0 and /2 in Fig.(1). Indeed the code calculate for the charge appearing in the system $$ the values of -0.1475 and -0.1399, very close to the values of -0.1455 and -0.1389 calculated for $$ at =0 and /2 as we previously mentioned. And the quantity$$was negligible by an order of magnitude ($$and $$at =0 and /2 respectively).

In our code, we allow at $$for the possibility of plasma to flow across the boundary. Note that for our present set of the parameters the electric field $$calculated remained negligible from =0 to /2 (see Fig.(6) below), and small around = .

Around = ,the charge appearing in the system is not exactly equal to the charge collected at the cylindrical vessel wall, but shows a small oscillation due the fact that there is a small plasma circulation at$$. In this case the electric field $$ is calculated from Eq.(41). At t=500, the value calculated for the right hand side of Eq.(41) is$$, while the value calculated by the code is $$(which is the value we see in the right figure in Fig.(1) for the = curve ). The agreement is very good.

Figure 1.

Electric field E_r at $$ (full curve), $$(broken curve), and $$ (dashed-dotted curve), at t=495 (left figure) and t=500 (right figure).

Figure(2) shows the potential at = 0 (full curve), = /2 (broken line) and = (dashed-dotted line), at time t = 495 (left figue) and t=500 (right figure). The curves at = 0 and = /2 show negligible oscillation, while the curve at = shows a small oscillation. Since as indicated in Fig.(1) the E_r profiles are essentially constant, especially for 0<R-r<75 along the gradient, then the oscillation of the curves at = (dashed-dotted line) is taking place in such a way that the slope of the potential $$ remains constant for 0<R-r<75. Only for 75<R-r<175 does the small oscillation of the curve at = in Fig.(2) (dashed-dotted line) translates into a small oscillation of the electric field E_r..

Figure(3) gives at time t = 500 and for = 0 the electric field E_r (full curve). Also at = 0 the dashed-dotted curve in Fig.(2) gives the Lorentz force, which in our normalized units is given by$$, and the broken curve gives the pressure force$$, $$, with the following definition:

$$E47

$$E48

Figure 2.

Potential at $$ (full curve), $$(broken curve) and $$ (dashed-dotted curve) at t=495 (left gigure) and t=500 (right figure).

Figure 3.

Plot at $$ of the electric field E_r (full curve), the Lorentz force $$ (dashed-dotted. curve), and the pressure force $$ (broken curve). The curve -n_i/2 is plotted for reference ( dashed- 3 dotted curve). At time t =500.

The $$ term (broken curve) shows a very good agreement along the gradient with the solid curve for$$, and the Lorentz force appears negligible along the gradient. In a region of about two gyroradii from the wall (around 40 Debye lengths from the wall), we have small irregular oscillations in space (and time), the accuracy of $$ being degraded by the division with a very small value of the density $$ appearing close to the vessel surface. To avoid this problem, we plot in Fig. (4) the quantities$$, $$at $$ ( note that$$). We see that there is a very nice agreement for the relation $$ along the gradient (the density $$ is also plotted in Fig.(4) to locate the different profiles with respect to the gradient). The Lorentz force term remains negligible in the gradient region, and remains very small in the bulk at 75<R-r<175. The $$ term is zero in the bulk since n_i and T_i are essentially flat in that region. We have seen in this case that the total charge $$ appearing in the system and calculated by the code by integrating the charge as in Eq. (40) and by calculating $$from Poisson equation, is essentially equal to the electric field$$. We note that curves similar to Fig.(4) can be calculated at different angles$$, showing $$ essentially balanced by the$$.

Figure 4.

Plot at $$ of $$ (solid curve), $$(dash-dot curve), and $$ (broken-curve). The curve –n_i/10 is plotted for reference (dashed-3 dotted curve). At time t=500.

The initial density profiles are given in Eq.(42), so the initial charge is zero. Fig.(5) shows the charge density $$ at the time t = 500 and $$ (full curve),$$ (broken curve) and $$(dashed-dotted curve). It is interesting to note that the code was able to maintain the steep density gradients (see Fig.(4) for a plot of$$). The stable in time steep density profile changes in space rapidly along the gradient over an ion orbit size ($$), and the relaxation of the steep gradients during the simulation determines the charge density$$. The charge density is important along the gradient at the plasma edge. The electrons, described by the guiding center equation given in Eq.(3) cannot compensate along the gradient the charge separation caused by the finite ion gyroradius, which results in the charge density we see in Fig.(5).

Figure 5.

Charge density at $$ (full curve), $$(broken curve) and $$ (dashed-dotted curve), at t=500.

We have noted in Fig.(1) the constant value of the $$profiles, especially for 0<R-r<75 along the gradient, translates the constant slope of the potential. However, the potential profile at = , although keeping essentially the same slope for 0<R-r<75, shows an oscillation in time. This results in a variation with an electric field $$ which, as we mentioned before, is small for the parameters we are using. Figures (6) presents this electric field $$ at t=485 and t=495. It shows no oscillation $$ in at = 0 (full curve, essentially equal to zero), a negligible oscillation at = /2 (broken curve), and a small oscillation for the = (dashed-dotted curve), on the high field side of the cylinder. The curve at = is interesting, it shows the oscillating curve for$$ having a constant flat value in the gradient region, and a linear variation in the inner region.

Figure 6.

Electric field $$ at $$ (full curve), $$(broken curve), and $$ (dash-dot curve), at t=485 (left figure) and t=495 (right figure).

Figure (7) gives at t = 500 the temperature T_ir as defined in Eq. (43), for = 0 (full curve), = /2 (broken curve) and = (dashed-dotted curve). The dashed-3 dotted curve is for n_iT_ir, which is essentially the same for all angles . Figure (8) presents similar results for T_i as defined in Eq.(43). Note in Fig.(7) and Fig.(8) the division by the very small density at the edge gives irregular oscillations. However the dashed-3 dotted curve in Fig.(7) and Fig.(8) is very smooth, and is for n_iT_ir and n_iT_i respectively, which removes the problem of the division by the very small density at the edge. We present in Fig.(9) the total pressure $$ for = 0 (full curve), = /2 (broken curve) and = (dashed-dotted curve). We see that these curves are essentially identical, for the parameters we are using in the present simulation there is no variation in for the total pressure term.

Figure 7.

Temperature T_ir at $$ (full curve), $$(broken curve) and $$ (dashed-dotted curve). The dashed -3 dotted curve is for n_iT_ir, which is essentially the same for all$$. Time t=500.

Figure 8.

Temperature $$ at $$ (full curve), $$(broken curve) and $$ (dash-dot curve). The dashed-3 dotted curve is for$$, which is essentially the same for all$$. Time t=500.

Figure 9.

Pressure $$at $$ (full curve), $$(broken curve) and $$ (dashed-dotted curve), at t= 500.

Figure(10) presents the plot of the component of the$$ drift, which in our normalized units is written$$, at = 0 (full curve), = /2 (broken curve) and = (dashed-dotted curve). The curves for 0<R-r<75 have reached a stable equilibrium as discussed in Fig.(1) for E_r, while towards the center for 75<R-r<175 the curves at = shows a small steady state oscillation around zero, as previously mentioned for E_r in the right figure in Fig.(1). Note in Fig. (10) that this azimuthal (i.e. poloidal)$$ drift is flat close to the edge, and is below the acoustic speed at = at the edge, and above the acoustic speed at = 0 and = /2 at the edge (velocities are normalized to the acoustic speed C_s). We plot in Fig. (11) the total poloidal current (in the direction)$$, where the diamagnetic drift is$$. We see that the total current is essentially zero, i.e. the profile is adjusting itself so that the $$ drift and the diamagnetic drift are essentially equal and opposite (in our units, the total poloidal current is$$, essentially equal to zero, showing a small oscillation around zero in Fig.(11), at the left boundary). So the $$ drift is balanced fairly well by the diamagnetic drift. This is the result we get if we calculate the poloidal current:

$$E49

We get indeed a negligible value for$$. (note that$$).

We note that due to the small$$field (see Fig.6), the $$drift has a small oscillating component $$ in the radial direction, with a curve similar to what is presented in Fig.(6). This radial oscillation is negligible around $$ and$$, and is small on the high field side around$$, as previously discussed for $$ in Fig.(6).

Figure 10.

Plot of $$ for $$ (full curve), $$(broken curve) and $$ (dashed-dotted curve). At time t=500.

Figure 11.

Plot of the total poloidal current $$ for $$ (full curve), $$(broken curve) and $$ (dashed-dotted curve). At time t=500.

4. Conclusion

We have presented in this work the self-consistent kinetic solution for the problem of the generation of a charge separation and an electric field at a plasma edge, under the combined effect of a large ratio of the ions’ gyroradius to the Debye length $$ (equal to 20 in the simulation we have presented) and a steep density gradient, when the electrons which are bound to the magnetic field cannot compensate along the gradient the charge separation due to the finite ions’ gyroradius. In the cylindrical geometry considered, a fully kinetic equation has been used to describe the ions, and a guiding center equation has been used to describe the electrons bound to the magnetic field. These equations have been solved using the method of characteristics (Shoucri, 2008a,b,c,d, 2009a). The numerical method used for the solution in cylindrical geometry, based on an integration of the equations along the characteristics coupled to a two-dimensional interpolation (Shoucri et al., 2004) applied successively in configuration space and in velocity space, is producing accurate results.

The problem of the formation of a charge separation is of great importance in the study of the H-mode physics in tokamaks. We have considered the case where the gradient in the density profiles is located in front of a floating cylindrical vessel. So the charge appearing in the system is essentially equal to the charge collected on the walls of the floating vessel. The solution shows in the radial direction that the electric field along the gradient is balanced by the radial gradient of the pressure, and the total poloidal current is essentially zero. The present results where electrons are described by a guiding center equation are close to what has been previously reported when assuming the electrons stationary and frozen by the magnetic field, and the code shows that a solution with a steep gradient, maintaining a charge separation, where the electron and ion densities vary rapidly over an ion gyroradius, is possible. Also the calculation with the present set of parameters which allow a small value of the poloidal field $$to exist, shows the presence of a small oscillation of $$on the high field side around $$(see Fig.(6)), to which is associated a small radial oscillation due to the radial component of the $$drift, equal to $$ in our normalized units.

The present code is a step closer to a code which will include “neoclassical” effects due to toroidal geometry, which can play a role in this problem, such as the neoclassical enhancement of the classical ion polarization drift, or the neoclassical damping of poloidal flows (Stix, 1973, Hirshman 1978, Waltz et. al., 1999).