Contents

• Cross-field charge transport by the diocotron instability in pulsar magnetospheres with gaps
• 1 Introduction
• 2 Model and underlying assumptions
• 3 Method of solution for non-linear simulations
• 4 Non linear simulations
• 5 Non linear development of the diocotron instability in a charge-fed disk
• 6 The quasilinear model
• 7 Results of the quasilinear analysis
• 8 General conclusions
• References
• 9 Online Material
• Appendix A: Angular momentum conservation and charge confinement
• Appendix B: Numerical methods for 2D simulations
• Appendix C: Numerical solution of the quasilinear equations

A&A 411, 203-213 (2003)
DOI: 10.1051/0004-6361:20031239

Cross-field charge transport by the diocotron instability in pulsar magnetospheres with gaps

J. Pétri ¹ - J. Heyvaerts ¹ - S. Bonazzola ²

1 - Observatoire de Strasbourg, 11 rue de l'Université, 67000 Strasbourg, France
2 - DARC, Observatoire de Meudon, place Jules Jansen, 92195 Meudon, France

Received 28 June 2002 / Accepted 30 July 2003

Abstract
In a previous work, we have shown by linear analysis that a thin charged disk in differential rotation in the magnetosphere of a neutron star with vacuum gaps is unstable to a collisionless instability induced in non-neutral plasmas by differential rotation, the diocotron instability. In this paper we study the long-time-scale evolution of this instability in the non-linear regime by means of both direct numerical simulations and a quasilinear model. We show that, when the disk is externally fed with charged particles produced by a moderate pair creation activity in the magnetosphere, the diocotron instability causes diffusion of the charged particles across the magnetic field lines outwards. An equatorial cross-field electric current is observed to form, carrying a net charge flux radially outwards. This constitutes a hitherto ignored charge transport mechanism in the pulsar magnetosphere. We briefly discuss how this turbulent charge transport mechanism could bear on the problem of electric current closure in pulsar's magnetospheres.

Key words: stars: pulsars: general - plasmas - magnetohydrodynamics (MHD) - methods: numerical

   
1 Introduction

   
1.1 Gaps in the closed pulsar magnetosphere

The electromagnetic activity of a pulsar is due to particle acceleration and pair creation cascade in its magnetosphere. The vicinity of the neutron star crust just above the polar caps, at which open field lines are rooted, is one particular region where such processes occur, but it is not the only one. Vacuum gaps are also likely to open in the vicinity of the null surface in the magnetosphere, where the Goldreich Julian charge density (Goldreich & Julian 1969) vanishes (Holloway 1973; Michel 1979). Smith et al. (1998, 2001) stress the fact that the concept of an entirely plasma filled closed magnetosphere is, for this reason, untenable since gaps would unstably widen. Closed magnetospheres with possible large vacuum gaps must then be considered. These gaps could be the source of the very high energy pulsed gamma ray emission presently observed up to tens of GeV in seven high energy pulsars (Thompson 2001). These energies are not yet high enough, however, to rule out the possibility that this radiation be produced at polar caps. The geometry of neutral-surface gaps and the electric tension across them self-consistently adjust so that, at equilibrium, there is no tendency for charge extraction out of the surfaces of the plasma-filled volumes. If the gaps are locally subject to pair creation activity, their thickness should adjust to marginal stability to pair avalanche creation. An entirely plasma-filled magnetosphere, or one which contains only very small vacuum gaps, is (almost) corotating and charged with the Goldreich-Julian corotation charge density. It is therefore electrically equivalent to the interior of the star. In our work, we regard such regions as being parts of a "renormalized'' star. The "renormalized star radius'' is the radius of the plasma-filled region about the star. However, vacuum sparking by pair avalanches cannot be active at any distance from the star. At some distance, a region should be met, where pair creation activity becomes weak, or inexistent, and large vacuum gaps are indeed present (Michel 1991; Pétri et al. 2002a). In-between these two regions, an intermediate region exists where there still is pair creation activity, but at a rate unable to maintain the flux tubes entirely filled with lepton plasma. If pairs can be removed from these flux tubes by some cross-field diffusion process, an electric circuit could be set up inside the closed magnetosphere. This current circulation could connect to the current system associated with the pulsar wind. This idea extends the notion of outer gaps proposed by Cheng Ho & Ruderman (1986) to the closed magnetospheric region.

   
1.2 Current circulation in the pulsar magnetosphere and open problems

When the magnetic field of the neutron star is dipolar and the magnetic and rotation axis are almost aligned, the emission of a charged plasma wind from the polar caps constitutes a net loss of charge. It is not granted that the system remains electrically neutral (Jackson 1979) at the scale of the star or at the scale of the light cylinder. A return current, which compensates for this loss of charge must self-consistently develop somewhere in the magnetosphere to ensure a stationary state, maintaining the total charge of the neutron star and its electrospheric environment constant. A number of ideas have been proposed to understand how this return current could physically arise. Mestel et al. (1979) proposed that charged particles cross magnetic field lines near the light cylinder as a result of the radiation reaction force. Beskin et al. (1983) have shown that in a boundary layer near the light surface of a magnetosphere filled with an electron-positron plasma, the electric drift approximation breaks down and particles are strongly accelerated cross-fields. This produces a current which closes the field-aligned current system flowing in the magnetosphere. More recently Contopoulos et al. (1999) have determined the MHD structure of the pulsar magnetosphere and of its relativistic wind by demanding that the wind solution be regular at the light cylinder. Their result features a charge-carrying wind emanating from the polar caps. This current returns in a current sheet which forms in the equatorial plane and bifurcates along the interface between open and closed field lines inside the light cylinder. No specific physical mechanism is proposed by these authors to explain the physical origin of this return current. If it is to consist of outflowing particles, there should be some region in the return current path that manages to provide charges of the appropriate sign, possibly at pair creation regions. Such regions would be akin to the outer gaps that Cheng et al. (1986) have suggested. These authors argue that the loss of charges along some field line near the light cylinder which are of a sign opposite to that of the Goldreich Julian charge density at the foot point of this line must induce the formation of a sparking gap on this open field line located in some vicinity of the light cylinder. The existence of such outer gaps would be necessary to feed the pulsar wind with particles of the appropriate sign in these regions. Little has been done up to now to describe self-consistently the electric and flow structure of such regions, although the analysis of the radiation that could be emitted by them has been developed in some detail (see for example Harding's review 2001).

   
1.3 On the possible role of turbulent cross-field charge transport

If charge transport would, for some reason, occur across closed magnetic surfaces, the return current mentioned above could be provided by particles originating from the closed magnetosphere, in particular from pair-creating gaps located inside this region. Such a transport, which could result from instabilities, would cause particles to diffuse to the edge of the closed region and then escape. This relatively new idea has also been suggested by Spitkovsky & Arons (2002). Specifically, a diffusive cross-field charge transport such as this could result from the presence of electrostatic turbulence induced by the development of the so-called diocotron instability in the differentially rotating part of the electrosphere.

As by now, velocity shear instabilities in non-neutral electrospheres of neutron stars have still been given little attention, although they are potentially very important, just as they are in laboratory devices (Davidson 1990). The purpose of this paper is to show that cross-field transport of charges in a closed magnetosphere with large gaps is a natural consequence of the differential rotation of the electrically disconnected part of its electrosphere. Electrically disconnected regions are defined here as being regions separated from corotating regions by vacuum gaps along the magnetic field lines which join them to the star. The differentially rotating structure is subject to the diocotron instability (Spitkovsky & Arons 2002; Pétri et al. 2002b). This instability is specific to non-neutral plasmas. It taps on the charge density inhomogeneity. In a non-neutral plasma there can only be equilibrium when some plasma motion provides the Lorentz force necessary to balance the space charge force. Some angular velocity profile is then necessarily associated with any equilibrium density profile. In the presence of a non axisymmetric perturbation, the perturbed charge density distribution is associated with some displacement of the isocontours of the density profile. Near the radius where this perturbation corotates with the equilibrium drift velocity of the non-neutral plasma, the charge density evolves, in the corotating frame, essentially secularly, which may cause instability (Davidson 1990). Once charge overdensities are formed, the charged particles are trapped in their vicinity because the electric drift motion carries them along equipotentials. The diocotron instability is of a collisionless and electrostatic nature. It is not physically similar to the Kelvin Helmholtz instability because plasma inertia plays no role in it. However, it is described in its linear stage by an equation which is identical to the one which describes the linear development of the Kelvin-Helmholtz instability.

Our purpose in this paper is to show that, in its non-linear stage, the diocotron instability can indeed give rise to an average cross-field turbulent transport of charge. We find that this happens when some regions of the electrospheric equatorial belt are fed with charge from some external source, which may occur due to, for example, moderate pair creation activity in the gaps overlying it. The construction of a global current circulation model incorporating this effect is an ambitious undertaking that is left for future work.

Spitkovsky & Arons (2002) have studied the diocotron instability in connection with the pulsar wind formation by way of numerical simulation, using a particle-in-cell code that they have adapted for this purpose. Their simulation has analyzed the development of the diocotron instability in a domain which begins quite far from the neutron star, at about 0.1 light radius, and extends to the light cylinder. The inner spherical boundary surface in their calculation acts as a source and sink of charged particles, feeding the electrospheric system in the simulation domain. Their initial condition consist of an electrosphere with large gaps, which then develops the diocotron instability, causing outwards extension of the charge distribution across field lines. This increases the amount of charge present at larger radii and reduces the size of the gaps as compared to the initial situation. Partial filling of these gaps results from the diffusion of charged particles. A complete filling of the magnetosphere by this process is however unlikely and is in fact not observed in their simulations, where some difference persists between the corotation density and the actual density in the equatorial belt. Complete filling would bring the differential rotation, which causes this transport, to an end. Spitkovsky and Arons suggest that the instability allows charges to diffuse up to the light cylinder and to fill a sizeable fraction of the electrosphere. Their simulation is carried out in the electrostatic approximation. The inclusion of magnetic field perturbations would be desirable near the light cylinder.

In this paper we want to examine the development of the diocotron instability in its non-linear stage under conditions which differ somewhat from those considered by Spitkovsky and Arons.

We studied the outcome of the instability both when the equatorial belt is electrically isolated and when it is charge-fed. Simulations carried in the case of an equatorial belt left in isolation have shown that no significant radial current caused by the diocotron fluctuations develops. Therefore we have choosen not to report on these simulations in this paper. We have observed that the electric charge gathers in this case in large localized entities of strong density, which eventually saturate and orbit the neutron star. This behaviour of an isolated equatorial belt reflects the limitations to charge transport implied by upper bounds to the amount of charge that could possibly leave the system under such conditions. These bounds have been derived analytically, in agreement with similar results previously obtained in the field of non-neutral plasma physics in the context of electronic tube physics (O'Neil 1980; O'Neil & Smith 1992; Davidson 1990). In Appendix A we extend these classical results to the present geometry, which is slightly different from that considered by O'Neil (1980). The rather stringent limits derived apply only to an isolated system. They do not restrict charge transport when the differentially rotating belt structure is somehow fed with an external source of charges. We therefore concentrate in this paper on this latter situation, which is astrophysically more relevant and also more likely to be met. We have found that in this case secular cross-field charge transport occurs and that the disk is maintained in a diocotron unstable state.

We did not follow Spitkovsky and Arons in assuming that the corotating part of the magnetosphere freely provides plasma to the non-corotating part. We view charge feeding as the result of moderate pair creation activity in the vicinity of the innermost part of the equatorial belt, for which we have derived a simple representation. We adopt a simplified model that is amenable to a direct numerical solution of the equation of motions, avoiding the need for using particle simulation methods. The charge concentrations and vortices produced by the instability are then more clearly apparent than they are in a particle simulation and they can be studied more easily. The price to pay is a reduction of the dimensionality of the system: we reduce the near-equatorial plasma system to a thin disk. We comment on the implications of this in Sect. 2.

The calculations of Krause-Polstorff & Michel (1985) as well as our own work (2002a) have shown that an electrosphere in axisymmetric electrostatic equilibrium consists of corotating domes above polar caps, charged with a sign which depends on that of the scalar product  $\vec{\Omega} \cdot \vec{B}$ at the pole, and a belt charged with the opposite sign about the equator. This belt consists of an inner corotating part and an extended tail separated, along field lines, from the corotating electrosphere by a large vacuum gap. The tail is therefore in differential rotation. Its radial extent depends on the total charge of the system and varies between a few times the renormalized radius of the star and infinity, when the total charge exceeds three times the reference charge $Q_{\rm c} = {{4 \pi}\over 3}~\epsilon_0~\Omega_*~B_*~R_*^3$. The differentially rotating belt is thin at moderate and large distances from the star. Initially at least, it is safe, in these regions, to regard it as a flat structure. Although our study considers a range of distances from the star, our electrostatic approximation does not allow us to closely approach the light radius.

We have conducted a 2D simulation to show how turbulent charge diffusion proceeds and have resorted to quasilinear theory to conduct a longer time scale analysis. This quasilinear approach averages out the transport over the azimuthal direction and describes the evolution of the mode amplitudes by linear growth corresponding to the instantaneous azimuthally averaged charge density profile.

   
1.4 Organization of the paper

In Sect. 2 we present our simplified picture of charge distribution in the neutron star's electrosphere with gaps, discuss its implications and present our mathematical description of the non-linear development of the diocotron instability. In Sect. 3 we describe the main aspects of our numerical methods, the technicalities being presented in Appendix B. In Sect. 4 we present the system that we have simulated and the diagnostics that have been performed. Our simulations of a charge-fed system are presented in Sect. 5. The analysis of its long term behaviour required calculations on a longer time scale than could be achieved by direct 2D simulations. Therefore we have developed a quasilinear model which is described in Sect. 6. Its solution has been obtained by the numerical methods described in Appendix C. The results are presented in Sect. 7. Our general conclusions concerning the conditions under which the diocotron instability can be effective in transporting charges across magnetic field lines in the pulsar's electrospheres and the possible consequences of this are presented in Sect. 8.

   
2 Model and underlying assumptions

Here, we study the long-time evolution of the diocotron instability (Pétri et al. 2002b) due to non-linear effects by means of numerical simulations, simultaneously solving the Poisson and the continuity equations. We make the same assumptions as in Pétri et al. (2002b). Namely, we assume that:

• the disk is thin and finite in extent;
• the particles are well within the light cylinder (non relativistic approximation);
• the magnetic field $\vec{B}_0$ is dipolar and constant;
• consequently the perturbations are only of an electrostatic nature;
• the motion of charges, with velocity $\vec{v}$, is treated at the electric drift approximation;
• these motions are restricted to the disk's plane.

The assumption of a thin disk has been adopted for the reasons explained in the introduction. It is supported by results concerning the geometry of the equilibrium electrosphere (Krause Polstorff & Michel 1985; Pétri et al. 2002a). These equilibrium solutions have a charged equatorial belt separated by large gaps from polar corotating domes. The belt has an inner corotating part and an extended tail. This latter part is separated along field lines from the corotating electrosphere by a large vacuum gap. This is the reason for its differential rotation. The radial extent of this tail region varies with the total charge of the system. It ranges between a few times the renormalized star radius, i.e. the radius of the region which is entirely plasma-filled, and infinity (or, more physically, the light radius) when the total charge exceeds three times the reference charge $Q_{\rm c} = {{4 \pi}\over 3}~\epsilon_0~\Omega_*~B_*~R_*^3$.

The electric potential in the gaps of the system consisting of an equatorial belt of a finite thickness surrounding a neutron star differs from the potential in the simplified system where the neutron star is supposedly surrounded by a flat disk, meant to represent the finite thickness belt. The difference is the largest in regions where the belt has an extension perpendicular to the equatorial plane comparable to the radial distance to the star, for example near the junction of the differentially rotating belt with the corotating magnetosphere. We adopt for the flat disk an initial surfacic charge density distribution identical to that which exists in the 3D equilibrium model. The surfacic charge distribution in the 3D structure is defined by integrating the volume charge distribution along field lines. The initial electric potential is then made consistent with this flattened charge distribution. As a result the rotation profile differs from that in the three dimensional equilibrium model in regions where the belt is thick. This difference may be large and counterrotation may even be locally seen in the innermost region of the disk. A perfect matching of the disk rotation to the star's rotation cannot be achieved in the simplified geometry. We have therefore left a little vacuum space between star and disk near its inner boundary. In the thin regions of the belt, which are distant from the star, the alteration of the 2D velocity profile as compared to the 3D equilibrium one remains limited. The description of these regions, which are most relevant to charge transport, by a flat model is then satisfactory.

As far as the study of perturbations is concerned, a thin disk approximation is appropriate for modes of an horizontal wavelength much larger than the disk thickness. This is the case for the low and moderate values of the azimuthal number m, which correspond to the linearly most unstable perturbations and develop in the non linear regime. The analysis of plasma motion in the flat disk approximation amounts to performing a latitudinal average over the belt structure, or, quite similarly, to discussing the effects of only those modes which have an azimuthal wavelength significantly larger than the belt's thickness. These modes must be the ones which also have the largest scale length in the radial direction. They are therefore expected to be those which contribute most to charge transport. We believe that keeping detailed track of the three-dimensional structure of the charge distribution in the belt is not essential in assessing the effects of the diocotron instability for charge transport in the thin regions of this structure. It would of course be more important if one were to provide a full geometrical description of the structure of the electrosphere. Hereafter, the equatorial belt will be refered to as being "the equatorial disk''. The neglect of diamagnetism is justified in regions of the disk where plasma motion is far from being relativistic. This assumption breaks down in the vicinity of the light cylinder.

The surfacic charge density distribution in the disk, $\sigma$, evolves in time according to the following system:

    

$$\frac{\partial \sigma}{\partial t} + {\rm div}\ (\sigma ~ \vec{v}) = 0$$ (1)

$${\varepsilon_0}~ \Delta \Phi + \sigma ~ \delta(z) = 0$$ (2)

$$\vec{E} = -~ \nabla \Phi$$ (3)

$$\vec{v} = \frac{\vec{E} \wedge \vec{B}_0}{B_0^2}\cdot$$ (4)

The electric potential and the electric field are denoted by $\Phi$and $\vec{E}$ resp. They are split into two parts, the first one, associated with the star, being labeled by a star, *, and the second one, associated with the disk, being labeled by D. Thus  $\Phi=\Phi_*+\Phi_D$ and  ${\vec E}={\vec E}_*+{\vec E}_D$.

   
3 Method of solution for non-linear simulations

   
3.1 The problem

The system (1)-(4), which contains eight unknown scalar quantities, can be reduced to just the Poisson Eq. (2) for the potential $\Phi$and the continuity Eq. (1) for the surfacic charge density $\sigma$ by substituting Eqs. (3) and (4) in (1). In cylindrical coordinates, the latter can be written as:

$$% \frac{\partial \sigma}{\partial t} + \frac{1}{r} \frac{\p... ...{r} \frac{\partial }{\partial \varphi}(\sigma~ v_\varphi) = 0.$$ (5)

The total electric field is expressed in terms of the potential as  $\vec{E} = -\nabla \Phi_D - \nabla \Phi_*$ where the stellar potential of the star, $\Phi_*$, is given, in the presence of the equatorial disk with total charge Q_D, by

$$% \Phi_*(R,\theta) = \frac{Q - Q_D - Q_{{\rm im}D}} {4~\pi~\... ...{3}~\Omega_\ast~B_\ast~R_\ast^5~\frac{P_2(\cos~\theta)} {R^3}$$ (6)

$Q_{{\rm im}D}$ is the charge of the image of the disk in the spherical star surface. This expression is valid even for non-axisymmetric charge distribution in the disk because the stellar potential remains axisymmetric, whatever this charge distribution may be. The components of the drift velocity, given by Eq. (4), can be written as:

  

$$\displaystyle{-\frac{1}{r~B_0}~ \frac{\partial\Phi}{\partial\varp... ... \displaystyle{2~\frac{r^2}{B_*~R_*^3}~ \frac{\partial\Phi_D}{\partial\varphi}}$$                                     v_r = (7)

$$v_\varphi \displaystyle{\frac{1}{B_0}~ \frac{\partial \Phi}{\partial r} }$$ =

$$\displaystyle{ -2~\frac{r^3}{B_*~R_*^3}~ \frac{\partial\Phi_D}{\p... ...Q_{{\rm im}D}}{2~\pi\varepsilon_0~B_*~R_*^3}~r + \Omega_*~\frac{R_*^2}{r}}\cdot$$ = (8)

The stellar magnetic field  $\vec{B}_0$ is dipolar and constant in time. In the equatorial plane, it is given by

$$% {\vec B}_0({\vec r}) = - \frac{B_*~R_*^3}{2~r^3}~{\vec e}_z.$$ (9)

The solution of the Poisson equation can be expressed in terms of the Green's function G, the specific expression of which is presented in Sect. B.2 of Appendix B, as

$$% \Phi_D(\vec{r}) = \int\!\!\! \int_D G(\vec{r}\vert\vec{r}~')\ \sigma(\vec{r}~') \ {\rm d}S_D.$$ (10)

The integration in Eq. (10) is over the entire disk. For an axisymmetric system, the charge distribution in the disk is time-independent and the radial drift velocity vanishes. Any evolution must result from non-axisymmetric perturbations.

The problem boils down to solving Eqs. (10), (5) with Eqs. (7)-(8). Our numerical methods to achieve this are described in Appendix B. The Sect. B.1 of this appendix presents the efficient and highly accurate method that we have devised for solving Poisson's equation by performing the integration in (10). To exactly conserve the total electric charge when solving Eq. (5), we have used a finite volume method rather than a finite difference scheme. Details are given in Sect. B.4 of Appendix B.

   
3.2 Boundary and initial conditions

Appropriate boundary conditions must be imposed in the azimuthal and radial directions. All physical quantities are periodic functions of the azimuthal angle $\varphi$. Boundary conditions in the radial direction are however more difficult to specify. To allow charge to spread out of the initial boundaries of the disk, we define a computational domain which is radially more extended, by a factor two or three, than the initial size of the charge distribution. This extra extension is chosen in view of the desired duration of the simulation. Taking a grid domain larger than the disk's extension causes some loss of accuracy and some waste of calculation time because some part of the CPU time is used in advecting vacuum. But it avoids the need to precisely follow the evolution of the edge of the charge distribution, which develops a complicated shape. At the outer boundary of the grid, we impose outgoing wave conditions to avoid the propagation of any signal from the exterior to the interior of the grid and to prevent parasitic reflexion at this boundary. This allows charges to escape out of the computational domain. At the inner boundary of the grid, at r=R_*, the potential of the disk is strictly zero, which means that the radial drift velocity vanishes and that particles cannot cross this surface. So, we take a solid wall condition at this boundary. This is at variance with the work of Spitkovsky & Arons (2002) who feed their computational domain with charge at the inner boundary. Boundary conditions are implemented with help of ghost cells. More details on this in relation to the CLAWPACK code are to be found in LeVeque's book (1997).

We finally specify initial conditions. For this we need to know the function  $\sigma_0(r)$, the surfacic charge density in the disk, at equilibrium which we take from our 3D equilibrium calculations (Pétri et al. 2002a). We initiate unstable motions by adding some non axisymmetric initial perturbation. We weakly excite initially a few modes m and let the system evolve for a time long enough to observe a significant redistribution of the charge density on the disk.

3.3 Testing of the numerical scheme

We have checked our algorithm by comparing the results which it provides in the linear regime with our earlier results (Pétri et al. 2002b) which have been obtained by an entirely different method. The eigenvalues found by running the present code agree to within a few percent with those found by our earlier approach. Some more details are given in Appendix B, Sect. B.5.

   
4 Non linear simulations

   
4.1 The simulated system

We have simulated a disk which suffers charge feeding from an external source with particles of the sign of its own surface charge, which in our simulations is positive. By this we mean to represent the effect of electron positron pair creation in the gaps existing between the domes and the equatorial belt (Pétri et al. 2002a). Feeding the disk with charges by this process could maintain its structure diocotron unstable at a low level, such that the lepton plasma production would give rise to a cross-field charge flux, even in the inner magnetosphere. Introducing a source of charges in this model is a simple first step towards the study of the interplay between pair creation and the diocotron instability.

The charge deposition on the disk is accompanied by the creation of an equal amount of charge of opposite sign. If it is assumed that the system suffers no net charge loss, this complementary charge would appear as surface charge on the stellar surface due to the electrostatic influence exerted by the disk. If, by contrast, it is assumed that this complementary charge, created simultaneously with the charge feeding the disk, is lost in some pulsar wind emitted from the polar caps, there would be a net charging of the star-electrosphere system, until the disk has grown so large that it reaches the light cylinder and the deposited charge is also lost at this outer boundary to the wind. In our simulations this eventual equilibrium regime has not been reached, although some substantial charge escape through the disk occured on the long term (see Sect. 7).

The assumption made concerning the fate of the negative charge, leading to a constant or increasing charge of the system, has no effect on its non-linear dynamics, as can be seen from Eqs. (7) and (8). The total charge of the system only affects the azimuthal velocity by adding a constant value to the rotation rate. We have studied the influence of the initial disk extension on the development of the diocotron instability. Several initial density profiles have been considered, along with different initial radial dimensions of the disk.

   
4.2 Physical quantities used for diagnostic

We monitored the time and space variations of a few physical quantities which have been expressed in terms of Fourier components of the surfacic charge density and potential. We follow the variations of the rotational speed and charge density, but also of other quantities, such as the mean radial current, the electrostatic potential energy, the charge flux or the mean radial extension of the disk. The mean radius  $R_{{\rm av}}$ of the charge distribution is defined here as the average value of r weighted by the normalized charge density distribution. We similarly define an average value of r², $R_{{\rm 2,av}}$, and, from it and  $R_{{\rm av}}$, a standard deviation, E_R, which is a measure of the fluctuations in the distance to the origin of the isocontours of the charge density distribution. The mean radial current across a circle  $\mathcal{C}$ of radius r at time tis:

$$% <i>(r,t) = \int_0^{2\pi}\!\! \sigma(r,\varphi,t) ~v_r(r,\varphi,t) ~r~{\rm d}\varphi .$$ (11)

The total charge having crossed the circle  $\mathcal{C}$ during the total time T of the simulation is the integral from zero to T of this mean radial current. We also considered the electrostatic potential energy:

$$% E_p = \frac{1}{2}~\int \! \! \! \int\sigma~\phi~{\rm d}S.$$ (12)

   
5 Non linear development of the diocotron instability in a charge-fed disk

In the study of Spitkovsky & Arons (2002), charges are injected at the inner boundary of the computational domain. This, however, cannot occur self-consistently when the electric potential is maintained at the prescribed, and axisymmetric, corotation value, as it would be, for example, inside the solid star (see Sect. 3.1). It is our view that charge-feeding of the disk is more likely to be the result of marginal pair creation activity at the intermediate region separating the fully filled part of the magnetosphere from the part which supports large vacuum gaps. This intermediate region is expected to be situated close to the renormalized star's surface. Very close to the fully filled region, the gap sustains but a small potential drop. At larger distances, it becomes broad, but inactive to pair-creation. We have determined the extent of this charged particle deposition region from an analysis of the opacity of the gaps of our equilibrium models to the gamma ray curvature photons emitted by particles accelerated in these gaps (Pétri et al. 2002a). The energy reached by these particles may be limited by radiation reaction. To simulate this pair creation effect near the inner boundary of the differentially rotating belt, we have introduced a source of charge which is assumed to feed the disk with particles of the same sign as its own charge density and with a given radial profile. The radial distribution of the adopted source function s(r)(Eq. (13)) is shown in Fig. 1.

$$% s(r) = 0.2 \ \left(\frac{1}{3} - \frac{0.5}{r}\right) \exp(-0.2~r).$$ (13)

The parameters in Eq. (13) have been chosen to represent charge deposition in a region the extent of which has been determined from our analysis of pair creation opacity. This source term is added on the right of the advection Eq. (5), which becomes:

$$% \frac{\partial \sigma}{\partial t} + \frac{1}{r} \frac{\pa... ...\frac{\partial }{\partial \varphi}(\sigma~ v_\varphi) = s(r).$$ (14)

Figure 1: Radial profile of the source term s(r) in the continuity equation.

The disk initially extends from  r₁=1.5 R_* to  r₂=50 R_*, where R_* is the renormalized star radius. The grid domain extends from R₁=1.01 R_* to  R₂=100 R_*. We could run the simulation from t=0 to t=8, with a time step  $\Delta t=1$. The results are shown in Fig. 2. On this time scale no small scale structure is observed to form and the system very quickly evolves to a non-axisymmetric configuration in which all the charge accumulates in a small region of strong density. The overall charge distribution shows no obvious spreading during the duration of the simulation, as can be seen from the analysis in Fig. 3. Most particles remain in the region  $[R_*, \ 20~R_*]$. The rotation profile flattens at distances larger than 20 R_* while, at distances less than 20 R_*, the gradient steepens due to the growth of the charge density at this location. The mean current is very intense when the charge redistribution begins to take place. It becomes much smaller later on. The total charge of the disk grows linearly with time due to charge feeding, indicating that negligible charge is lost through the boundaries of the computational domain during this early phase. As a consequence, the potential energy grows too. The flux of charge, which remains large in the region of charge deposition, is zero at both boundaries, indicating the absence of any escaping charge.

Figure 2: Non linear evolution of the diocotron instability for a total charge of the system equal to $Q_{\rm c}$ and with an external source of charges s. Initially, the disk extends from  r1=1.5 R* to  r2=50 R* while the grid extends from  R1=1.01 R* to  R2=100 R*. The pictures are taken from t=0 to t=8 with a time step  $\Delta t=1$, from left to right, and from top to bottom. The inner circle corresponds to r=R* while the outer one to r=33 R*.

Figure 3: Main characteristics of the non linear evolution of the diocotron instability corresponding to the case presented in 2.

This situation cannot be maintained on the long term. It is clear that the disk has not reached a quasi stationary state when our simulation terminates. All simulations mentioned above have been made with a grid of  $50\times256$ cells. Some extra runs have been made with an improved resolution of  $200\times256$ cells with identical results.

To observe the eventual stationary state, in which significant loss of charges through the grid boundaries is to occur, the duration of the simulation should have been be increased by a factor of order 10 or more. This implied too long a computing time and too large data files. We therefore had to turn to a more efficient approach to the long term evolution of the system. We compromised by developing a quasi linear model as an alternative to the 2D non linear simulation. The quasilinear model reduces the description of the system to just a few azimuthally averaged quantities. In the present geometry it is one-dimensional in space and therefore much less expensive in computer time than the non-linear 2D model (we gained a factor 10 in CPU time). In a number of cases this model very well reproduces the results of the full non-linear 2D calculation in the time interval of overlap.

   
6 The quasilinear model

When only a few unstable modes are excited, the instability evolves to the formation of coherent non linear structures in the charge density. When enough modes are excited and their phases are sufficiently incoherent, however, it is possible to describe the evolution in an approximate way in terms of a few azimuthally averaged macroscopic quantities only, such as the charge density and the electric potential. In this section, following Davidson (1990), we formulate a quasilinear model of the diocotron instability.

   
6.1 Derivation of the quasilinear equations

The quasilinear model only describes approximately the evolution of a few azimuthally averaged quantities. Such a reduction is justified when the correlation length of the electric field fluctuations is small as compared to the radius, or when their correlation time is small as compared to the characteristic evolution time of the density profile. In the present situation, the applicability of the quasilinear model might be questioned in view of the fact that structures with a certain degree of spatial coherence are observed to form in our non linear simulations. However their characteristic "life time'' appears to be short as compared to the characteristic time observed for the evolution of the average density structure. We checked a posteriori that in the region of time overlap the results of the fully two-dimensional simulation and of the quasilinear model agree well enough. The agreement in fact turned out to be very good, due to this relatively small value of the correlation time.

In the quasilinear model, any physical quantity  $\psi(r,\varphi,t)$ is split into two parts, its azimuthally averaged value  $<\psi>(r,t)$, which we generally denote by bracketing, like <>, and its fluctuation with respect to this average, $\delta\psi(r,\varphi,t)$, which is just the difference between $\psi(r,\varphi,t)$and  $<\psi>(r,t)$. This splitting is applied to the potential and to the surface charge density. The angular averaging operator has the following properties:

$$\frac{\partial}{\partial \varphi}<\psi> ~ = ~ <\frac{\partial \psi}{\partial \varphi}> ~ = 0$$ (15)

$$<\chi ~\frac{\partial \psi}{\partial \varphi}> ~ = - <\psi ~\frac{\partial \chi}{\partial \varphi}>.$$ (16)

The potential and the surface charge density are then split into average value and fluctuations as indicated, and these split expressions are introduced in the Poisson and charge conservation equations, which are then azimuthally averaged. This gives:

$$% \frac{1}{r}~\frac{\partial}{\partial r}\left( r~\frac{\par... ...}{\partial z^2} = -\frac{<\sigma>}{\varepsilon_0}~\delta(z)$$ (17)

$$% \frac{\partial <\sigma>}{\partial t} = \frac{1}{r}~\frac{\... ...igma~ \frac{\partial \delta\phi}{\partial \varphi}> \right).$$ (18)

The effective charge flux, which is the average value of the radial current density at radius r, is seen to equal $-1/B_0~\times$ $<\delta\sigma~(\partial \delta\phi/\partial \varphi)>$. It is the most interesting and meaningful quantity in Eq. (18). It results from nonlinearities involving the fluctuations. The azimuthally averaged components of the electric drift velocity are given by:

$$% <v_r>~=-\frac{1}{r~B_0}<\frac{\partial <\phi>}{\partial \varphi}> ~= 0$$ (19)

$$% <v_\varphi>~=\frac{1}{B_0}~\frac{\partial<\phi>}{\partial r}\cdot$$ (20)

Equations (17)-(20) describe the non linear evolution of the mean potential, density and velocity profiles in response to the fluctuations  $\delta\sigma$and  $\delta\phi$. The equations describing the non linear evolution of these fluctuations are obtained by subtracting Eqs. (17) and (18) from Eqs. (2) and (5) resp. This gives:

$$% \Delta \delta\phi = -\frac{\delta\sigma}{\varepsilon_0}~\delta(z)$$ (21)

 

$$% \left( \frac{\partial}{\partial t} + \frac{<v_\varphi>}{r}~ \fr... ...\frac{\partial}{\partial r} \left( \frac{<\sigma>}{B_0} \right) = \qquad \qquad \frac{1}{r} \left[ \frac{\partial\delta\phi}{\partial \varphi} \f... ...} <\delta\sigma~ \frac{\partial \delta\phi}{\partial \varphi}> \right) \right].$$ (22)

These equations give a closed description of the non-linear evolution of the instability, which, at this stage, is still fully equivalent to the set of the exact continuity and Poisson equations. As the density profile evolves, the perturbations grow and $<\sigma>$ readjusts according to (18) in response to these unstable perturbations.

The quasilinear approximation consists in neglecting the bilinear terms on the right of Eq. (22). This is valid, for example, for small enough fluctuations. The linearized continuity equation for  $\delta\sigma$ is then solved in conjunction with Poisson's equation for  $\delta\phi$ and the resulting expressions substituted in Eqs. (17) and (18). We so obtain a closed equation for the average density profile $<\sigma>$. Fluctuations $\delta\sigma$and  $\delta\phi$ are found as in the linear analysis, by complex Fourier expansion with respect to the azimuthal variable $\varphi$ (as in the first equation of (B.1)). These Fourier components satisfy the following equations, which result from Eqs. (23)-(24) by dropping the terms bilinear in  $\delta\phi$ and  $\delta\sigma$. This simplification is the essence of the quasilinear approximation:

  

$$\frac{1}{r}~\frac{\partial}{\partial r}\left( r~\frac{\partial \d... ... \frac{m^2}{r^2} \delta\phi_m = -\frac{\delta\sigma_m}{\varepsilon_0}~\delta(z)$$ (23)

$$\left( \frac{\partial}{\partial t} + {\rm i}~m~\frac{<v_\varphi>}... ...ta\phi_m ~ \frac{\partial}{\partial r} \left( \frac{<\sigma>}{B_0} \right)\cdot$$ (24)

The slow evolution of the mean density $<\sigma>$ can then be expressed in terms of these Fourier variables as

$$% \frac{\partial<\sigma>}{\partial t} = - \frac{1}{r}~ \fra... ...^{+\infty} {\rm i}~ m ~ \delta\sigma_m^*~\delta\phi_m \right).$$ (25)

Equations (23) and (25) form a nonlinear coupled system jointly describing the evolution of $<\sigma>$,  $\delta\phi$ and  $\delta\sigma$, at the quasilinear level of approximation. To analyze the behavior of its solutions, we restrict the analysis to only slowly growing modes ( $\gamma_m>0$), with a WKB type of time dependence of the form:

 

$$% \delta\sigma_m(r,t) = \delta\sigma_m(r) ~ {\rm e}^{-{\rm i}~\int_0^t (\omega_m(t') + {\rm i}~\gamma_m(t')) {\rm d}t'}$$ (26)

$$\delta\phi_m(r,t) = \delta\phi_m(r) ~ {\rm e}^{-{\rm i}~\int_0^t (\omega_m(t') + {\rm i}~\gamma_m(t')) {\rm d}t'}.$$ (27)

This is unlike linear analysis because the eigenfrequency $\omega_m$and the growth rate $\gamma_m$ now vary slowly in time in response to the change of the mean profile $<\sigma>$. Substituting (26) in the continuity Eq. (24), we get:

$$% (\omega_m + {\rm i}~\gamma_m - m<\Omega>)~\delta\sigma_m = ... ...rac{\partial}{\partial r} \left( \frac{<\sigma>}{B_0} \right)$$ (28)

where  $<\Omega>(r,t)$ is the mean angular velocity in the disk at time t. Integrating Poisson's equation by means of Green's functions, we finally obtain for the Fourier components of the density perturbation a Fredholm integral equation of the third kind:

$$% \omega~\delta\sigma_m = m~<\Omega>~\delta\sigma_m + \frac... ...t_{R_1}^{R_2} G_m(r,r')~\delta\sigma_m (r') ~ r' ~ {\rm d}r'.$$ (29)

This eigenvalue problem has been met in our linear analysis (Pétri et al. 2002b). Now, however, the density and angular velocity profile are allowed to vary in time. The eigenvalues $\omega$and the eigenfunctions  $\delta\sigma_m$ also vary. Inserting (26) in the continuity Eq. (25) we finally obtain a one-dimensional diffusion equation for the azimuthally averaged density:

$$% \frac{\partial}{\partial t}<\sigma>\ = \frac{1}{r}~\frac{\... ...ial r} \left( \frac{<\sigma>}{B_0} \right) \right) + s(r,t).$$ (30)

The source term s(r,t) has been added on the right of this equation. The symmetry properties of the eigenvalues and eigenfunctions, i.e. $\omega_{-m} = - \omega_m$ and $\gamma_{-m} = \gamma_m$, allow us to express the diffusion coefficient $\mathcal{D}$ as

 

$$% \mathcal{D}(r,t) \!= \! \frac{2}{B_0(r)} \! \sum_{m=1}^{+\infty... ...delta\phi_m(r,t)\vert^2} {(\omega_m(t) - m \! <\Omega>(r,t))^2 + \gamma_m(t)^2}$$ (31)

while  $\vert\delta\phi_m(r,t)\vert^2$ evolves according to the wave kinetic equation:

$$% \frac{\partial}{\partial t}\vert\delta\phi_m(r,t)\vert^2 = 2~ \gamma_m(t) \vert\delta\phi_m(r,t)\vert^2.$$ (32)

The mean radial flux of electric charge, which is the surfacic current density, is then given by:

$$% <f_{\rm s}>(r) = - \mathcal{D}(r,t) ~ \frac{\partial}{\partial r} \left( \frac{<\sigma>}{B_0} \right)$$ (33)

and the mean radial current through the circle of radius r is:

$$% <i_{\rm s}>(r) = - 2~\pi~r~ \mathcal{D}(r,t) ~ \frac{\partial}{\partial r} \left( \frac{<\sigma>}{B_0} \right)\cdot$$ (34)

The quasilinear evolution of the disk can be followed by numerically solving the linear integral Eq. (29) and the nonlinear diffusion Eq. (30). The numerical procedures employed are described in Appendix C.

Figure 4: Characteristics of the quasilinear evolution of the diocotron instability with an additional external source. Initially, the system has a total charge given by  $Q_{{\rm tot}}=3~Q_{\rm c}$. Time is in units of  $\Omega _*^{-1}$, the distance in units of the stellar radius R*, the growth rate in $\Omega _*$, the current is expressed in  $\varepsilon_0~\Omega_*^2~B_*~R_*^3$ and the diffusion pseudo-coefficient  $\mathcal{D}$ in  $\Omega_* B_*~R_*^2$.

   
7 Results of the quasilinear analysis

To compare the results provided by this quasilinear model with those of the non-linear simulation, the same parameters have been adopted in both cases. The initial conditions are the same as in the non-linear simulation. The only difference lies in the boundary conditions, which are outgoing wave conditions for the non-linear simulation, and fixed boundary conditions for the quasilinear model. This boundary condition is better suited for the diffusion equation. The instability is initiated by the excitation of all modes, with an initial amplitude of  $10^{-5}~\sigma_0(r)$.

We chose the same source as in Sect. 5. The velocity gradient is enhanced by the effect of the source of charges (see Fig. 4). Due to the instability, the particles are transported outwards, causing an expansion of the disk. The ratio of the diffusion coefficient to  $\Omega_* B_* R_*^2$, which is shown in Fig. 4, stabilizes at a fairly large value of order 10^-2. The associated current increases continuously with time. Although all modes were initially excited, most of them progressively disappear and only a few survive when the state of marginal stability is reached, maintaining diffusion across the magnetic field lines. An interesting phenomenon appears near t=10. The charge deposited on the disk begins to escape out in significant amount and irreversibly. Here, we find a situation in which the system fed with charges succeeds in generating, by the effect of the diocotron instability, a sustained flux of charges outwards. No stationary state has been reached, since, at the end of the simulation, charge is still accumulating over the disk. Since a continued growth of the accumulated charge is impossible, we expect that at extremely long times a truely stationary state would be reached. The non linear simulation shows charge structures which, at any given time, appear to be spatially rather coherent, with well defined charge condensations which retain their identity for a few rotational times. However, after a time of this order, these structure disappear and are replaced by similar ones, with unrelated phases. We believe that this poor coherence time, as opposed to the rather large spatial coherence, is the reason why a quasilinear description meets some success.

The increase in the charge content of the disk does not cause a breakdown of the thin disk assumption. This can be judged from the results obtained. Dividing the surfacic charge density by the differential Goldreich-Julian density calculated from the rotation profiles that we obtain from the flat model, we get an estimate of the corresponding 3D belt thickness h(r). This estimate shows that h/r remains less than 0.75 at distances larger than about 5 renormalized stellar radii R_* and decreases rapidly with increasing r. The region where h/r remains smaller than this is a bit larger when the charge of the system increases instead of remaining constant. Thus the flat disk approximation which we have adopted remains valid up to the end of our quasilinear simulation in the outer regions of the system.

   
8 General conclusions

In the absence of an external source of charges, calculations not reported in this paper have shown no significant current development in the equatorial plane, in agreement with the magnetic confinement property discussed in Appendix A, which places an upper limit to the fraction of the charge that can leave the system under these conditions.

The evolution of the system is entirely different in the presence of an external source of charges feeding the disk with charges of its own sign. This situation would be representative of neutron star magnetospheres with gaps which are partially pair-producing. These sources of charges would be located close to the region of contact of the differentially rotating belt with the corotating magnetosphere. The mean radial extension of the electrospheric disk first continuously increases due to charge diffusion induced by the diocotron instability. On the long term, a radial charge flux is organized and maintained, the charge diffusion coefficient keeping a sizeable value, due to continued charge deposition. The instability has the effect of redistributing the electric charge over the disk, which in return causes the growth rates of unstable modes to decrease from an initial value of the order of the star's rotation frequency. The system approaches marginal stability, but does not exactly reaches it, due to continued charge feeding. In this asymptotic situation, the growth rates become small, but remain just large enough to maintain diffusion of the charged particles across the magnetic field lines. For a sufficiently large duration of the simulation, a permanent outwardly oriented electric current has been observed to form in the equatorial plane, carrying away a large fraction of the charge delivered by the source to the disk. This type of current circulation in the closed magnetosphere, due to the presence of the diocotron instability, has been given but little attention before, although it could bear on the problem of current closure of the electric circuit of an active pulsar magnetosphere because it breaks the need for currents to be field-aligned in the closed magnetosphere. It has not been possible, during the finite time of our simulations, to observe our simulated system settling down to a stationary state in which the charge deposition by the source would be balanced by the radial current induced by the diocotron instability. The construction of an electric current model for an active pulsar involving this turbulent transport mechanism in magnetically closed zones with gaps and connecting self-consistently to magnetically open wind regions across the light-cylinder is an ambitious task which is left for future work. Our present contribution establishes, in the framework of the adopted model, the relevance of cross field conduction by non-neutral plasma turbulent processes. We believe that this idea of a diocotron-turbulent charge transport in the pulsar's electrosphere, when extended to a 3D electrospheric structure including the presence of charged polar domes, will open the way to a new understanding of the electrodynamics of the pulsar's magnetosphere.

Acknowledgements

We thank R. Walder and D. Folini for their help in finding and using the CLAWPACK package and H. Baty for help in the numerical solution of the quasilinear system.

References

 
• Abramowitz, M., & Stegun, I. A. 1965, Handbook of mathematical functions (N.Y.: Dover)  
• Beskin, V. S., Gurevich, A. V., & Istomin, I. N. 1983, ZhETF, 85, 401 NASA ADS  
• Cheng, K. S., Ho, C., & Ruderman, M. 1986, ApJ, 300, 500 NASA ADS  
• Contopoulos, I., Kazanas, D., & Fendt, C. 1999, ApJ, 511, 351 NASA ADS  
• Davidson, R. C. 1990, Physics of non-neutral plasmas (Addison-Wesley Publishing Company)  
• Goldreich, P., & Julian, W. H. 1969, ApJ, 160, 971  
• Harding, A. 2001, in High Energy Gamma Ray Astronomy, ed. F. A. Aharonian, & H. J. Völk (American Institute of Physics), 115  
• Holloway, N. J. 1973, Nat. Phys. Sci., 246, 6 NASA ADS  
• Jackson, E. A. 1979, ApJ, 227, 266 NASA ADS  
• Krause-Polstorff, J., & Michel, F. C. 1985, A&A, 144, 72 NASA ADS  
• Leveque, R. J. 1997, Computational Methods for Astrophysical Fluid Flow., Saas-Fee Lecture Notes  
• Mestel, L., Phillips, P., & Wang, Y. M. 1979, MNRAS, 188, 385 NASA ADS  
• Michel, F. C. 1979, ApJ, 227, 579 NASA ADS  
• Michel, F. C. 1991, Theory of Neutron Star Magnetospheres (Univ. Chicago Press)  
• O'Neil, T. 1980, Phys. Fluids, 23, 2216 NASA ADS  
• O'Neil, T., & Smith, R. A. 1992, Phys. Fluids B, 4, 2720 NASA ADS  
• Pétri, J., Heyvaerts, J., & Bonazzola, S. 2002a, A&A, 384, 414 NASA ADS  
• Pétri, J., Heyvaerts, J., & Bonazzola, S. 2002b, A&A, 387, 520 NASA ADS  
• Press, W., Teukolsky, S., Vetterling, W., & Flannery, B. 1997, Numerical Recipes in C (Cambridge Univeristy Press)  
• Smith, I. A., Michel, F. C., & Thacker, P. D. 2001, MNRAS, 322, 209 NASA ADS  
• Smith, I. A., Thacker, P. D., & Michel, F. C. 1998, Rev. Mex. Astronom. Astrofis, 7, 211  
• Spitkovsky, A., & Arons, J. [astro-ph/0201360]  
• Thompson, D. J. 2001, in High Energy Gamma Ray Astronomy, ed. F. A. Aharonian, & H. J. Völk (American Institute of Physics), 103

  
9 Online Material

   
Appendix A: Angular momentum conservation and charge confinement

   
A.1 The electromagnetic angular momentum of the electrospheric system

The effect of the diocotron instability on the structure of an isolated electrospheric disk is under the control of the conservation law of angular momentum applied to the system consisting of the particles and the electromagnetic field. O'Neil (1980), in a seminal paper, first derived an upper bound on the fraction of electrons that may reach the outer wall of a charge-separated cylindrical plasma column confined by an uniform magnetic field. The plasma is bounded by a perfectly conducting and particle-absorbing wall at some finite distance. The wall is assumed to pass no energy flux. The bound relies on the conservation of angular momentum and energy for the system consisting of the electrons and the electromagnetic field. The method consist in deriving upper bounds to the material and electromagnetic contribution to the total angular momentum. Here, we extend these results to the present geometry, which is different from that considered by O'Neil (1980) in that the magnetic field is not uniform, the distribution of electrons is in the equatorial plane, so that the electromagnetic contribution to the angular momentum expresses differently. Also, the system, although it is unbound externally, is limited by a perfectly conducting central sphere. We neglect however electron inertia and consider electrostatic perturbation fields only. When particle inertia is neglected, the total angular momentum of the system reduces to that of the electromagnetic field alone. In this appendix we show that this conservation law implies the constancy in time of the following integral:

$$% {\cal{I}} = \int_{R_1}^{R_2}\!\!\int_0^{2~\pi} \left( \frac... ...\frac{2}{R_*} \right)~\sigma^0(r)~r~{\rm d}r~{\rm d}\varphi .$$ (A.1)

From this we deduce in Sect. (A.2) below that the amount of charge $\Delta Q$ that can escape to infinity from an isolated electrospheric disk with initial charge Q is limited by

$$% \frac{\Delta Q}{Q} \leq 1 - \frac{1}{2~Q}~ \int_{R_1}^{R_2... ...ac{R_*^2}{r^3} \right)~\sigma^0(r)~r~{\rm d}r~{\rm d}\varphi .$$ (A.2)

If, initially, the charges are very close to the star, the double integral is close to 2 Q. The fraction of the charge allowed to escape to infinity is in this case very small. Conversely, if the disk initially extends far away from the star, this double integral becomes negligibly small and the maximal fraction approaches unity. The escape of a large fraction of the charges present in the electrosphere then becomes possible, which doesn't mean that it actually occurs.

Let us establish the constancy of the double integral (A.1) for isolated evolution. Let us recall that the momentum density of the electromagnetic field is:

$$% \vec{g} = \varepsilon_0~\vec{E}\wedge\vec{B}.$$ (A.3)

The angular momentum density of the electromagnetic field is then:

$$% \vec{l}=\varepsilon_0~\vec{r}\wedge(\vec{E}\wedge\vec{B}).$$ (A.4)

The system formed by the disk and the electromagnetic field is closed when there is no exchange of particles, nor any energy transfer with the exterior. The Poynting flux through the star surface is zero because the Poynting vector only has an azimuthal component there. Moreover, the electromagnetic field in the star, which results from the supposedly infinite conductivity of the crust of the neutron star, remains constant in time and thus plays no role in the energy balance. Neglecting the inertia of the charged particles, the conservation of the total angular momentum reduces to the conservation of the angular momentum of the electromagnetic field alone:

$$% \vec{L}_0 = \varepsilon_0\int\!\!\!\int\!\!\!\int_{\bar{*}} \vec{r}\wedge(\vec{E}\wedge\vec{B})~{\rm d}^3\vec{r}.$$ (A.5)

The integration is to be performed over all space outside the sphere of radius R_*. This region is denoted by $\bar{*}$. The magnetic field being axisymmetric, the only part of the electric field which contributes to Eq. (A.5) is the one which does not depend on $\varphi$. The axisymmetric part of the integrand in Eq. (A.5) is:

$$% \vec{r}\wedge(\vec{E}^0\wedge\vec{B}) - \vec{r}\wedge\left(\vec{\nabla}\Phi^0\wedge\left( \frac{1}{r}~\vec{\nabla}a\wedge\vec{e}_\varphi\right)\right)$$ =

$$(\vec{\nabla}\Phi^0\cdot\vec{\nabla}a)~ \vec{e}_z - \frac{z}{r} (\vec{\nabla}\Phi^0\cdot\vec{\nabla}a)~ \vec{e}_r.$$ = (A.6)

The potential $\Phi^0$ and the flux function a, defined by $\vec{B}=\frac{1}{r}\nabla a\wedge\vec{e}_\varphi$, are symmetric with respect to the equatorial plane, z=0. They are even functions of z, as is the gradient of the potential. Then, by parity, Eq. (A.5) reduces to:

$$% \vec{L}_0 = \varepsilon_0~\int\!\!\!\int\!\!\!\int_{\bar{*}... ...\nabla}\Phi^0\cdot\vec{\nabla}a)~{\rm d}^3\vec{r} ~ \vec{e}_z.$$ (A.7)

This vector is parallel to the polar axis z because the configuration is initially axisymmetric. The potential $\Phi^0$ can be split into two parts, the first one being associated with the star and the second one with the average radial distribution of the charge density, $\sigma^0(r)=\frac{1}{2~\pi}~\int_0^{2~\pi} \sigma(r,\varphi)~{\rm d}\varphi$. The expression (A.7) may be rearranged as:

 

$$% \vec{\nabla}\Phi^0\cdot\vec{\nabla}a {\rm {div}}~(\Phi^0~\vec{\nabla}a) - \Phi^0~\Delta~a$$ =

$${\rm {div}}~(\Phi^0~\vec{\nabla}a) - \Phi^0~\left(\Delta^*~a + \frac{2}{r}~\frac{\partial a}{\partial r}\right)$$ = (A.8)

where we have introduced the modified Laplacian, $\Delta^*$, defined by

$$% \Delta^* a = r~\frac{\partial}{\partial r} \left( \frac{1}... ...}{\partial r}\right) + \frac{\partial^2 a}{\partial z^2}\cdot$$ (A.9)

Since the magnetic field is assumed to remain potential (Sect. (2)) $\Delta^*~a=0$. Transforming the volume integral in the first term of (A.8) into a surface integral, we get

$$% \vec{L}_0 = \left( \varepsilon_0~ \int\!\!\!\int_{S_*}\Phi... ...ac{\partial a}{\partial r}~{\rm d}^3\vec{r}\right)~\vec{e}_z.$$ (A.10)

The electromagnetic field vanishes at infinity fast enough for the flux through the sphere at infinity to vanish. The surface integral in Eq. (A.10) has to be performed only on the sphere of radius R_*, denoted by S_*. The volume integral can be easily shown to be convergent because $\Phi^0\sim r^{-1}$, $\frac{1}{r}~ \frac{\partial a}{\partial r}\sim r^{-3}$. The surface integral is calculated by noting that the potential on S_*is the corotation potential and we get:

$$% \int\!\!\!\int_{S_*}\Phi^0~ \vec{\nabla}a\cdot {\rm d}\vec{S} =... ...}{4~\pi~\varepsilon_0} ~ \int\!\!\!\int_{S_*}\vec{\nabla}a\cdot {\rm d}\vec{S}.$$ (A.11)

The magnetic flux function a of the dipolar magnetic field is

$$% a(r,z) = \displaystyle{\frac{1}{2}~B_*~R_*^3~ \frac{r^2}{(r^2+z^2)^{3/2}}}$$ (A.12)

which implies that

$$\int\!\!\!\int_{S_*}a~\vec{\nabla}a\cdot {\rm d}\vec{S} = \frac{8~\pi}{15} ~B_*^2~R_*^5$$ (A.13)

$$\int\!\!\!\int_{S_*}\vec{\nabla}a\cdot {\rm d}\vec{S} = \frac{4~\pi}{3} ~B_*~R_*^3.$$ (A.14)

The volume integral I of (A.10) can be expressed from (A.12) and by using the sub-relativistic approximation to the Goldreich-Julian density as

$$% I = \frac{1}{\Omega_*}~\int\!\!\!\int\!\!\!\int_{\bar{*}} \Phi^0(\vec{r}) \rho_{{\rm GJ}}(\vec{r})~{\rm d}^3\vec{r}.$$ (A.15)

To calculate the contribution of the disk to I we express the disk potential by means of the Green function G. Exchanging the order of integration between $\vec{r}'$ and $\vec{r}$, we obtain:

 

$$\frac{1}{\Omega_*}~\int_{R_1}^{R_2}\!\!\int_0^{2~\pi} \left\lbrac... ...G(\vec{r}\vert\vec{r}') ~\rho_{{\rm GJ}}(\vec{r})~{\rm d}^3\vec{r}\right\rbrace$$                                      I_D =

$$\sigma^0_D(\vec{r}')~{\rm d}^2\vec{r}'.$$ (A.16)

The triple integral in curly brackets in Eq. (A.16) is the potential created at $\vec{r}'$ by a volumic charge distribution  $\rho_{{\rm GJ}}(\vec{r})$ outside the star S_*. This potential vanishes at infinity and on S_*. This triple integral is calculated by noting that when all space is filled with the Goldreich-Julian density, the magnetic field lines are equipotentials. We then find the corotation potential:

$$% \Phi_{{\rm GJ}}(\vec{r}) \int_0^{+\infty}\int_0^{2~\pi} \!\!\int_0^\pi L(\vec{r}\vert\vec{r}') ~\rho_{{\rm GJ}}(\vec{r})~{\rm d}^3\vec{r}$$ =

$$\Omega_*~a + {\rm cste} = \frac{1}{2}~\Omega_*~B_*~R_*^3~\frac{\sin^2\theta}{R} + {\rm cste}.$$ = (A.17)

L represents the Green function G without the image counterpart. For more details, see B.2. The potential created by the external part of the star is :

$$% \int_{R_*}^{+\infty} \int_0^{2~\pi}\!\!\int_0^\pi L(\vec{r}\ver... ...!\int_0^\pi L(\vec{r}\vert\vec{r}') ~\rho_{{\rm GJ}}(\vec{r})~{\rm d}^3\vec{r}.$$ (A.18)

The potential generated by the internal charge of the pulsar, which consist of the central point charge and the charge distributed with the density  $\rho_{{\rm GJ}}$ is

$$% \int_0^{R_*}\!\!\int_0^{2~\pi}\!\!\int_0^\pi L(\vec{r}\vert\vec... ...lon_0~R'} - \frac{1}{5}~\Omega_*~B_*~R_*^5~\frac{P_2(\cos\theta')}{{R'}^3}\cdot$$ (A.19)

Therefore, we deduce that:

$$% \int_{R_*}^{+\infty}\!\!\int_0^{2~\pi}\!\!\int_0^\pi L(\vec{r}\... ...~\varepsilon_0}~\frac{3}{5}~R_*^2\frac{P_2(\cos\theta')}{{R'} ^3} + {\rm cste}.$$ (A.20)

This potential does not satisfy the boundary conditions because its value on the surface S_* is

$$\begin{eqnarray*}-\frac{2}{5}~\frac{Q_{\rm c}}{4~\pi~\varepsilon_0}~\frac{P_2(\cos\theta)}{R_*} + {\rm cste}. \end{eqnarray*}$$

To obtain a solution which vanishes on S_* we add an harmonic potential  $\Phi_{{\rm CL}}$ vanishing at infinity and equal to  $(-\frac{2}{5}~\frac{Q_{\rm c}}{4~\pi~\varepsilon_0}~ \frac{P_2(\cos\theta)}{R_*}+{\rm cste})$ on S_*:

$$% \Phi_{{\rm CL}}(R,\theta) = -\frac{2}{5}~\frac{Q_{\rm c}}{... ...\varepsilon_0}~R_*^2~ \frac{P_2(\cos\theta)}{R^3}+{\rm cste}.$$ (A.21)

The potential created by the Goldreich-Julian distribution which satisfies the required boundary conditions is:

 

$$% \int_{R_*}^{+\infty}\!\!\int_0^{2~\pi}\!\!\int_0^\pi G(\vec{r}\vert\vec{r}') ~\rho_{{\rm GJ}}(\vec{r})~{\rm d}^3\vec{r} \frac{Q_{\rm c}}{4~\pi~\varepsilon_0} \left( \frac{R_*^2}{{R'}^3} - \frac{1}{R'}\right)$$ =

$$P_2(\cos\theta').$$ (A.22)

It vanishes at R_* and at infinity. Coming back to the complete volume integral (A.15):

$$I = \frac{Q-Q_D-Q_{{\rm im}D}}{4~\pi~\varepsilon_0~\Omega_*}~ \in... ...\int\!\!\!\int\!\!\!\int_{\bar{*}} a~\rho_{{\rm GJ}}(\vec{r})~{\rm d}^3\vec{r}.$$ (A.23)

The total charge enclosed in the spherical shell between Rand  $R+{\rm d}R$ is zero because of its angular distribution proportional to  $P_2(\cos\theta)$. Moreover

$$% \int\!\!\!\int\!\!\!\int_{\bar{*}} a~\rho_{{\rm GJ}}(\vec{r})~{\rm d}^3\vec{r} = \frac{2}{5}~B_*~R_*^2~Q_{\rm c}.$$ (A.24)

Restricting to the equatorial plane, we get:

$$% \int_{R_*}^{+\infty}\!\!\int_0^{2~\pi}\!\!\int_0^\pi G(\vec... ...} \left( \frac{1}{2~r'} - \frac{R_*^2}{2~{r'}^3} \right)\cdot$$ (A.25)

Thus, the component along the $\vec{e}_z$ axis of the total angular momentum of the electromagnetic field is

 

$$% L_0 = \frac{4}{5}~B_*~R_*^2~Q_{\rm c} + \frac{1}{3}~B_*~R_*^2~( Q - Q_{\rm c} - Q_D - Q_{{\rm im}D} )$$

$$+ \frac{1}{6}~B_*~R_*^3~ \int_{R_1}^{R_2}\!\!\int_0^{2~\pi} \left... ...rac{R_*^2}{r^3} - \frac{2}{R_*} \right)~\sigma^0(r)~r~{\rm d}r~{\rm d}\varphi .$$ (A.26)

The sum of the charge of the disk and of its image is

$$% Q_D + Q_{{\rm im}D} = \int_{R_1}^{R_2}\!\!\int_0^{2~\pi} \... ... \frac{R_*}{r} \right)~\sigma^0(r)~r~{\rm d}r~{\rm d}\varphi .$$ (A.27)

Isolating in Eq. (A.26) that part of L₀ which doesn't depend on the surfacic charge distribution on the disk, we find that the angular momentum conservation law for the electromagnetic field reduces to the constancy of the following integral:

$$% \int_{R_1}^{R_2}\!\!\int_0^{2~\pi} \left( \frac{3}{r} - \f... ...\frac{2}{R_*} \right)~\sigma^0(r)~r~{\rm d}r~{\rm d}\varphi .$$ (A.28)

   
A.2 A general limit to the loss of charge by an isolated electrospheric disk

The conservation law of angular momentum, applied to the system consisting of the particles and the electromagnetic field, imposes general restrictions on the plasma motion if the system is isolated. It has been shown above that this conservation law implies the constancy of the integral (A.1) or (A.28). From this we can deduce very general results concerning the non-linear evolution of the diocotron instability. We get more insight by regarding the system as a collection of individual particles, rather than a continuous charge distribution. The constant  $\frac{2}{R_*}$ in the integrand of Eq. (A.28) can be ignored because the total charge of the system is conserved. Let  ${\rm d}Q(r,\varphi)$ be the charge element at  $(r,\varphi)$. Eq. (A.28) then reduces to:

$$% \int_{R_1}^{R_2}\!\!\int_0^{2~\pi} \left( \frac{3}{r} - \frac{R_*^2}{r^3} \right)~{\rm d}Q(r,\varphi)\ =\ \rm {constant}.$$ (A.29)

The constancy of this quantity can be interpreted in terms of discrete particles carrying all the same charge q. Suppose that in the initial state the disk contains N charges, labeled i, localized at points  $(r_{i0},\varphi_{i0})$. The conservation of the angular momentum as expressed in the form Eq. (A.29) can be restated for this set of discrete particles as the following relation, valid for any time t>0

$$% \sum_{i=1}^N \left(\frac{3}{r_i(t)} - \frac{R_*^2}{r_i(t)^3... ...}^N \left(\frac{3}{r_{i0}} - \frac{R_*^2}{r_{i0}^3} \right) q.$$ (A.30)

Assume that $\Delta N$ particles escape to infinity. Their contribution to the total angular momentum will then be zero. To compensate for the decrease of their contribution, the remaining charges have to recess towards the star. In the worst case, they will all reach the star's surface, at r=R_*. In this final extreme situation:

$$% \left\lbrace \begin{array}{ccc} r_i(t) & = & R_* {\rm\ \ ... ...= & +\infty {\rm\ \ for\ \ i}>N-\Delta N. \end{array} \right.$$ (A.31)

The angular momentum conservation (A.30) applies and we can write:

$$% \sum_{i=1}^{N-\Delta N} \left(\frac{3}{R_*} \!-\! \frac{1}{... ...eft(\frac{3}{r_i(0)} \!-\! \frac{R_*^2}{r_i(0)^3} \right)\cdot$$ (A.32)

Thus, the maximum number $\Delta N$ of charges that the disk may lose is limited by:

$$% \Delta N = N - \frac{1}{2}~\sum_{i=1}^N \left(\frac{3~R_*}{r_i(0)} - \frac{R_*^3}{r_i(0)^3} \right)\cdot$$ (A.33)

The upper limit to the fraction of particles which could escape to infinity is given in terms of the initial conditions, by:

$$% \frac{\Delta N}{N} \leq 1 - \frac{1}{2~N}~\sum_{i=1}^N \left(\frac{3~R_*}{r_i(0)} - \frac{R_*^3}{r_i(0)^3} \right)\cdot$$ (A.34)

Returning to the continuous description in terms of surface charge density  $\sigma_0(r)$, this is expressed as Eq. (A.2). This upper limit overestimates the amount of charge actually lost, because we have assumed that remaining charges have all gathered on the stellar surface at r=R_*. Eq. (A.2) is a constant, time independent, limit, whatever the non linear evolution of the disk, complicated though it may be.

   
Appendix B: Numerical methods for 2D simulations

   
B.1 Solving Poisson's equation

In this section, we describe our numerical method to solve Poisson's equation, that is, to evaluate the double integral (10). We need a routine which calculates quickly but accurately this double integral. To do so, we Fourier transform the azimuthal dependence. The potential and the density are expanded in a real Fourier series (in "$\cos$'' and in "$\sin$''), while the Green's function is expanded in a complex Fourier series:

 

$$G(r,\varphi\vert r',\varphi') = \displaystyle{ \sum_{m=-\infty}^{+\infty} G_m(r,r') \ {\rm e}^{{\rm i}~m(\varphi-\varphi')}}$$ (B.1)

$$\sigma(r,\varphi) = \displaystyle{ \sum_{m=0}^{+\infty} \sigma_m^c(r) \cos~m\varphi + \sigma_m^s(r) \sin~m\varphi }$$ (B.2)

$$\Phi(r,\varphi) = \displaystyle{ \sum_{m=0}^{+\infty} \Phi_m^c(r) \cos~m\varphi + \Phi_m^s(r) \sin~m\varphi }.$$ (B.3)

Substituting these expansions in (10) and integrating over $\varphi$, we obtain the Fourier coefficients of the potential in terms of those of the charge density and the Green's function:

$$% \Phi_m^{c/s}(r) = 2~\pi~\int_{r_1}^{r_2} G_m(r,r') \sigma_m^{c/s}(r') ~ r' ~ {\rm d}r'.$$ (B.4)

The evaluation of expressions like (B.4) is complicated by the fact that the Green's functions G_m have a logarithmic singularity on the diagonal r=r'. We present our method for handling it. Consider the mode m=0. The integral (B.4) is regularized by the following reformulation:

 

$$% \Phi_0(r) 2~\pi~\int_{r_1}^{r_2} G_0(r,r') ( \sigma_0(r') - \sigma_0(r) ) ~ r' ~ {\rm d}r'$$ =

$$+ 2~\pi~\sigma_0(r) ~ \int_{r_1}^{r_2} G_0(r,r') ~ r' ~ {\rm d}r'.$$ (B.5)

The first integral on the right of Eq. (B.5) is regular because  $G_0(r,r') (\sigma_0(r') - \sigma_0(r) )$approaches zero when  $r\rightarrow0$ and can be calculated by usual methods. The second integral on the right of Eq. (B.5) can be performed analytically. Let us define the residual of the Green's function by

$$% RG_m(r) = 2~\pi~\int_{r_1}^{r_2} G_m(r,r') ~ r' ~ {\rm d}r'.$$ (B.6)

The residual of the Green function for m = 0 can be expressed in terms of elliptic integrals of the first king K and of the second kind E as:

 

$$\frac{1}{\pi~\varepsilon_0} \left[ r_2 ~ E\left(\frac{r^2}{r_2^2}\right) - r ~ E\left(\frac{r_1^2}{r^2}\right) \right.$$                           RG₀(r) =

$$+ \left. r ~ \left( 1 - \frac{r_1^2}{r^2} \right) ~ K\left(\frac{... ...\frac{R_*}{r} \left( r_1~ E\left(\frac{R_*^4}{r^2~r_1^2}\right) \right. \right.$$

$$- \left. \left. r_2~ E\left(\frac{R_*^4}{r^2~r_2^2}\right) \right)\right]\cdot$$ (B.7)

Thus the constant coefficient in the Fourier series of the potential is found. For all other modes, m>0, the integral is regularized by:

 

$$% \Phi_m^{c/s}(r) 2~\pi\int_{r_1}^{r_2} [G_m(r,r')-G_0(r,r')]~ \sigma_m^{c/s}(r') ~r'~{\rm d}r'$$ =

$$+ 2~\pi\int_{r_1}^{r_2} G_0(r,r')~ \sigma_m^{c/s}(r') ~r'~{\rm d}r'.$$ (B.8)

The integrand of the first term on the right of Eq. (B.8) is continuous in the interval [r₁,r₂]. We show in B.2 that the function  (G_m(r,r')-G₀(r,r')) is regular in the square  [r₁,r₂]², even on the diagonal r=r'. The second integral reduces to the same form as for the case m=0 by simply replacing $\sigma_0$ by $\sigma_m$. The Green functions for non-vanishing m's can be expressed in terms of hypergeometric functions (Abramowitz & Stegun 1965):

 

$$\frac{1}{4~\pi~\varepsilon_0} ~ \frac{\Gamma(m+\frac{1}{2})}{\Gamma(m+1)\Gamma(\frac{1}{2})}$$                           G_m(r,r') =

$$\left[ \frac{r_<^m}{r_>^{m+1}} \ _2F_1\left( m+\frac{1}{2}, \frac{1}{2}; m+1; \frac{r_<^2}{r_>^2}\right) \right.$$

$$- \left. \frac{R_\ast^{2m+1}}{(r~r')^{m+1}} ~ _2F_1\left( m+\frac{1}{2}, \frac{1}{2}; m+1; \frac{R_\ast^2}{(r~r')^2}\right) \right]\cdot$$ (B.9)

We have noted  $r_<={\rm min}(r,r')$and  $r_>={\rm max}(r,r')$. These functions have been numerically calculated by the integration of a second order ordinary differential equation (see Press et al. 1998).

   
B.2 Green functions

In this section we develop a numerical method to calculate quickly and efficiently, by using a FFT, the integrals  RG_m,jⁱ needed in the eigenvalue problem (29). These integrals are defined by:

$$% RG_{m,j}^i= 2~\pi~\int_{r_i}^{r_{i+1}} G_m(r_{j+1/2},r')~r'~{\rm d}r'.$$ (B.10)

In the equatorial plane the Green function is:

$$% G(r,\varphi\vert r',\varphi')= L(r,\varphi\vert r',\varphi'... ...}{r'}~ L\left(r,\varphi\vert\frac{R_*^2}{r'},\varphi'\right).$$ (B.11)

Here we have introduced the function L which gives the potential at the point  $(r',\varphi',z')$ generated by an unit point charge at  $(r,\varphi,z)$:

$$% L(r,\varphi,z\vert r',\varphi',z') \frac{1}{4~\pi~\varepsilon_0} \frac{1}{\sqrt{r^2+{r'}^2+(z-z')^2}}$$ =

$$\times \frac{1}{\sqrt{1-k^2\cos~(\varphi-\varphi')}}$$ (B.12)

$$\sum_{m=-\infty}^{+\infty} L_m(r,r') ~ {\rm e}^{{\rm i}m(\varphi-\varphi')}.$$ = (B.13)

The parameter $k^2=\frac{2 r r'}{r^2+{r'}^2+(z-z')^2}$ is less than unity. The coefficients of the Fourier series G_m of our Green function G are related to those of L by:

$$% G_m(r,r')= L_m(r,r') - \frac{R_*}{r'}~L_m\left(r,\frac{R_*^2}{r'}\right)\cdot$$ (B.14)

We need to calculate the integrals:

$$% RG_m(r) = 2~\pi~\int_a^b G_m(r,r')~r'~{\rm d}r'.$$ (B.15)

From Eq. (B.14), we get:

$$RG_m(r) = 2~\pi~\int_a^b L_m(r,r')~r'~{\rm d}r' + 2~\pi~R_* \int_a^b L_m(r,\frac{R_*^2}{r'})~{\rm d}r'.$$ (B.16)

These integration can be done by means of the Fourier transform of $L(r,\varphi\vert r',\varphi')$. Indeed,

$$% \int_a^b L(r,\varphi\vert r',\varphi')~r'~{\rm d}r' = \sum_{m=0}^{+\infty} \int_a^b L_m(r,r')~r'~{\rm d}r'~\cos m~\psi .$$ (B.17)

with  $\psi=\varphi-\varphi'$. The integrals  $\int_a^b L_m(r,r')~r'~{\rm d}r'$ are then the coefficients of the Fourier series of  $\int_a^b L(r,\varphi\vert r',\varphi')~r'~{\rm d}r'$. If this latter expression is known analytically, the integrals  $\int_a^b L_m(r,r')~r'~{\rm d}r'$ can be obtained by a FFT. A similar remark applies to  $\int_a^b L_m(r,\frac{R_*^2}{r'})~{\rm d}r'$. Performing the integration analytically, we get respectively:

$$% I_1(r,\psi) \int_a^b L(r,\varphi\vert r',\varphi')~r'~{\rm d}r'$$ =

$$\frac{1}{4~\pi~\varepsilon_0} \left[ \sqrt{r^2+b^2- 2~r~b \cos\psi} \right.$$ =

$$-\left. \sqrt{r^2+a^2- 2~r~a\cos\psi} \right. + r~\cos\psi$$

$$\left. \times \ln\left( \frac{\sqrt{r^2+b^2- 2~r~b \cos\psi} + b - r~\cos\psi} {\sqrt{r^2+a^2- 2~r~a \cos\psi} + a - r~\cos\psi}\right) \right]$$ (B.18)

$$% I_2(r,\psi) \!=\! R_* ~ \int_a^b L\left(r,\varphi\vert\frac{R_*^2}{r'},\varphi'\right)~{\rm d}r'$$

$$\!=\! \frac{R_*}{4~\pi~\varepsilon_0~r^2} \left[ \sqrt{r^2~b^2 \!+\! R_*^4 \!-\! 2~r~b~R_*^2 \cos\psi} \right.$$

$$\!-\! \left. \sqrt{r^2~a^2 \!+\! R_*^4 \!-\! 2~r~a~R_*^2 \cos\psi} \right . + R_*^2~\cos\psi$$

$$\left. \!\times\!\ln\left(\!\frac{\sqrt{r^2~b^2 \!+\! R_*^4 \!-\!... ...! 2~r~a~R_*^2 \cos\psi} \!+\! a~r \!-\! R_*^2~\cos\psi}\!\right) \right]\!\cdot$$ (B.19)

The integral I₁ diverges logarithmically when $r\in[a,b]$and $\psi=0$ because it represents the potential of a charged wire extending from a to b and directed radially along  ${\bf e}_{\vec r}$, the potential being evaluated on the wire itself. We must carefully study the behavior of I₁ in the neighbourhood of $\psi=0$. For r<aand r>b, I₁ is continuous and regular for all $\psi$, there is no singularity in the potential. However, for r>b there is a difficulty to evaluate this function numerically because it is given by an indeterminate expression of the form  $\frac{0}{0}$. Treating this indetermination analytically, we get for r>b:

$$% I_1(r,0) \!=\! \int_a^b L(r,r')~r'~{\rm d}r' \!=\! \frac{1}... ...\varepsilon_0} \left[ a-b + r~\ln\frac{r-a}{r-b} \right]\cdot$$ (B.20)

To deal with the case a<r<b, where the denominator of the logarithmic function approaches zero, we isolate the singular part

$$% r~\cos\psi~\ln\left( \sqrt{r^2+a^2- 2~r~a \cos\psi} + a - r~\cos\psi \right)$$ (B.21)

which behaves, for  $\psi\rightarrow0$, as

$$% r~\cos\psi~\ln\left(\frac{r^2}{r-a}~(1-\cos\psi) \right).$$ (B.22)

In this case we introduce the regularized function

$$% \int_a^b L(r,r')~r'~{\rm d}r'+ \frac{1}{4~\pi~\varepsilon_0} ~r~\cos\psi~\ln\left(\frac{r^2}{r-a}~(1-\cos\psi)\right)$$ (B.23)

whose FFT can be calculated. Its singular part is given by

$$% \frac{1}{4~\pi~\varepsilon_0} ~r~\cos\psi~\ln\left(\frac{r^2}{r-a}~(1-\cos\psi)~ \frac{1}{\sqrt{r^2+a^2-2~a~r~\cos\psi}+a-r~\cos\psi}\right)\cdot$$ (B.24)

It approaches zero when $\psi$ vanishes. In order to find the Fourier series expansion of I₁ we have to find those of the function:

$$% \frac{1}{4~\pi~\varepsilon_0} ~r~\cos\psi~ \ln\left(\frac{r^2}{... ...}{r\!-\!a}\right) + \frac{1}{4~\pi~\varepsilon_0} ~r~\cos\psi~ \ln(1-\cos\psi).$$ (B.25)

The task reduces then to evaluating the integrals:

$$% J_m=\int_0^{2~\pi} \cos\varphi ~ \ln(1- \cos\varphi) ~ \cos m\varphi ~ {\rm d}\varphi .$$ (B.26)

Straightforward calculation shows that:

$$J_0 = -2~\pi$$ (B.27)

$$J_1 = -\pi~\left(\ln2+\frac{1}{2}\right)$$ (B.28)

$$J_m = -\frac{2~m}{m^2-1} ~\pi \ \ \rm {for} \ \ m\ge2.$$ (B.29)

The second kind of integrals, I₂, corresponding to the potential created by the image of the charged wire, is regular for all r>R_*and for all $\psi$. Thus we can immediately perform the FFT on I₂.

   
B.3 The diagonal singularity

In this section, we show that the function  (G_m(r,r')-G₀(r,r')) is regular in the square  [r₁,r₂]², even on the diagonal r=r'. The singular term in (B.14) is contained in the function L_m(r,r') because for r and r' larger than R_*, $L_m(r,\frac{R_*^2}{r'})$ is regular. We first prove that  (L_m(r,r')-L_m-1(r,r')) is regular in  [r₁,r₂]².

$$L_m(r,r')-L_{m-1}(r,r') = \frac{1}{4~\pi~\varepsilon_0}~ \frac{1}... ...i}~m~\psi}-{\rm e}^{-{\rm i}~(m-1)~\psi}}{\sqrt{1-k^2~\cos\psi}} ~{\rm d}\psi .$$ (B.30)

Changing the variable by  $\psi=2~\varphi$, we get

$$L_m(r,r')-L_{m-1}(r,r') = \frac{1}{4~\pi~\varepsilon_0}~ \frac{1}... ...i}~(2~m-1)~\varphi}~\sin~\varphi}{\sqrt{1-k^2~\cos2~\varphi}} ~{\rm d}\varphi .$$ (B.31)

When r=r', the parameter k=1 and then:

$$lim_{~r'\to r}\; (L_m(r,r')-L_{m-1}(r,r') ) \frac{1}{4~\pi~\varepsilon_0~r}~ \frac{(-i)}{\pi} \int_0^{\pi} {\rm e}^{-{\rm i}~(2~m-1)~\varphi} ~{\rm d}\varphi .$$ = (B.32)

Finally, we find that this difference is real and finite and equal to:

$$% lim_{~r'\to r}\; ( L_m(r,r')-L_{m-1}(r,r')) = -~\frac{1}{2~m-1} ~ \frac{1}{2~\pi^2~\varepsilon_0~r}\cdot$$ (B.33)

By recurrence, it is straightforward to prove that the difference  L_m(r,r')-L₀(r,r') is regular on the diagonal r=r'. As a consequence, the difference  G_m(r,r')-G₀(r,r') is also regular on the diagonal r=r'.

Figure B.1: Some examples of density perturbations for the system with total charge equal to  $Q=Q_{\rm c}$. The Fourier coefficients in  $\cos m\varphi$ and in  $\sin m\varphi$ are presented for the modes m=2,5 and 10. The oscillation frequency and the shape of the envelop give the eigenvalues for the corresponding mode.

   
B.4 Solving the continuity equation

The continuity equation expresses charge conservation in the disk. To exactly conserve the total charge, we use a finite volume method instead of a finite difference scheme. Specific algorithms exist to solve this kind of problem. We used CLAWPACK, written by LeVeque, which is freely available on the web. It solves hyperbolic equations in one, two or three dimensions, using Godunov and MUSCL schemes with fractional step. CLAWPACK can handle, by means of Riemann solvers, general 2D advection equations of the form:

$$% \kappa(r,\varphi,t)~\frac{\partial q}{\partial t} + \frac{... ...artial r}(u~q) + \frac{\partial }{\partial \varphi}(v~q) = 0.$$ (B.34)

Equation (5) is of this form, with $q=\sigma$, u=r v_r, $v=v_\varphi$ and  $\kappa(r,\varphi,t)=r$. In the expression of the azimuthal component of the drift speed (8), the charge of the disk and of its image are given in terms of the Fourier coefficients of the density  $\sigma(r,\varphi)$. To evaluate these quantities, we first perform the integration in the azimuthal direction. Only the axisymmetric coefficient  $\sigma_{m=0}(r)$ survives this operation and we find

  

$$Q_D = \displaystyle{2~\pi~\int_{R_1}^{R_2} \sigma_{m=0}(r)~r~{\rm d}r}$$ (B.35)

$$Q_{{\rm im}D} = \displaystyle{-2~\pi~R_*~\int_{R_1}^{R_2} \sigma_{m=0}(r)~{\rm d}r}.$$ (B.36)

   
B.5 Testing the numerical scheme

We have checked our algorithm in the linear approximation. The growth rates and the eigenfrequencies obtained in the simulation have been compared to those found in our earlier linear analysis (Pétri et al. 2002b). We first tested the Godunov method, which is first order in space and time. The results have been found to be reliable for sufficiently small time steps. In the test, each mode mhas been made to evolve independently by removing the non linear terms from the equations. For too small time steps, the calculated growth rates have been observed to become too small, which is an effect of the numerical diffusion introduced by the space-time grid. The numerical diffusion coefficient is  $D_{{\rm num}} =({\Delta x^2}/{\Delta t})$. To prevent from this effect, we used a method of second order in space, for which numerical diffusion is negligibly small. This method, unlike the first order method, finds an approximation to the solution which is piecewise linear and not piecewise constant. Piecewise linear doesn't mean that the function is continuous. Indeed, it isn't because there is a jump in the solution at each cell interface. The slopes used to specify this approximate solution are to some extent arbitrary. However, their values must provide as good as possible approximations for the partial derivative  $\partial/\partial x$. We have observed that choosing certain slopes introduces nonlinearities in the discretized equation. Care must be taken to chose these slopes in such a way that no spurious non linear terms be introduced. Such terms are avoided by the Lax-Wendroff, the Beam-warming and the Fromm schemes. Beam-warming and Fromm give the best results.

 

 

Table B.1: Eigenfrequencies and growth rates for several modes. The values on the left of "/'' correspond to the simulations and the values on the right of "/'' correspond to the linear analysis.

Mode m	Eigenvalue $\omega$
	Rotation speed $\Omega_m$	Growth rate
	$m~\Omega_m = Re(\omega)$	$\gamma=Im(\omega)$
2	2.731 / 2.692	0.992 / 1.016
5	7.140 / 7.040	2.318 / 2.447
10	14.200 / 14.244	3.181 / 3.681

Taking care of these subtleties, we have run some simulations of the linear approximation to our system and estimated the eigenfrequencies and the growth rates by examination of the eigenfunctions. A sample is shown in Fig. B.1 for the modes 2, 5 and 10. The eigenvalues found by this method agree within a few percent with those obtained in our previous linear analysis (Table B.1). The growth of the amplitude is perfectly exponential, as can be seen in lower right panels in B.1. We note however small discrepancies between the growth rates obtained by these two methods, the most significant one arising for m=10 and the difference reaching 14%.

   
Appendix C: Numerical solution of the quasilinear equations

The quasilinear evolution of the disk can be followed by numerically solving the linear integral Eq. (29) and the nonlinear diffusion Eq. (30). The numerical procedures employed are briefly described in this section.

C.1 Diffusion equation

Choosing a numerical scheme for solving the diffusion equation is in general no problem. An implicit method, of first order in time, or the Crank-Nicholson scheme, of second order in time, works well in most cases where the diffusion coefficient explicitly depends on space and time. Here, however, this coefficient is part of the solution. The problem becomes nonlinear. Symbolically, we can write it down in the form:

$$% \kappa(x)~\frac{\partial u}{\partial t} = \frac{\partial}{... ...t( D(u,x,t) ~ \frac{\partial u}{\partial x} \right) + s(x,t).$$ (C.1)

We seek for solutions u in the interval I=[a,b] satisfying the initial condition  u(x,0)=u₀(x) and the boundary conditions  $u(a,\cdot)=u_G$ and  $u(b,\cdot)=u_D$. The functions $\kappa$, D and s are known and continuous on I. $\kappa$ has been introduced to easily adapt the procedure to different coordinates systems, cartesian or cylindrical.

An explicit method would impose a very small time step, because the stability criterion demands that  $\Delta t<(\Delta x^2/D)$. This maximum time step is much smaller than the characteristic time of evolution over which the system is to be solved. An implicit method, combined with an iterative process, has to be used. Suppose we know the solution u_iⁿ and the diffusion coefficient D_iⁿ at time t. The scheme to advance one time step from t to   $t+\Delta t$ is as follows:

1.

Make an initial guess for the diffusion coefficient at step n+1 by  $\tilde{D}_i^{n+1}=D_i^n$;

2.

Find the solution to the diffusion equation with this guess. This gives  $\tilde{u}_i^{n+1}$;

3.

Update the diffusion coefficient  $\tilde{D}^{n+1}=D(\tilde{u}^{n+1},x,t)$;

4.

Go back to step 2 until the convergence of $\tilde{u}$ is reached within a given accuracy.

Testing this scheme, we have met some difficulties in implementing the Crank-Nicholson method to solve the diffusion equation, the convergence, if any, being poor. We have therefore prefered a fully implicit method despite its first order accuracy in time.

C.2 Fredholm integral equation

In addition to the diffusion equation, we have to solve the Fredholm integral Eq. (29):

$$% \omega~\delta\sigma_m = m~<\Omega>~\delta\sigma_m + m ~ F(r) 2~\pi~\int_{R_1}^{R_2} G_m(r,r')~\delta\sigma_m(r')~r'~{\rm d}r'$$ (C.2)

where  $F(r)=\frac{1}{r}~\frac{{\rm d}}{{\rm d}r}\left(\frac{<\sigma>}{B_0}\right)$is a continuous function of the variable r. We have solved this equation when dealing with the linear analysis of the diocotron instability (Pétri et al. 2002b). The method was however very expensive in time, which would create difficulty here because the equation must be solved several times per time step. We therefore developed another method. We approximate  $\delta\sigma_m$with a piecewise constant function  $\delta\sigma_m(r)\approx\delta\sigma_m^i$for  $r\in[r_i,r_{i+1}]$. The corresponding quadrature formula is:

$$% 2~\pi~\int_{R_1}^{R_2} G_m(r_{j+1/2},r')~\delta\sigma_m(r'... ...rm d}r' \approx \sum_{i=0}^{n-1} \delta\sigma_m^i~RG_{m,j}^i.$$ (C.3)

The term  $RG_{m,j}^i= 2~\pi~\int_{r_i}^{r_{i+1}} G_m(r,r')~\delta\sigma(r')~r'~{\rm d}r'$ represent the residual functions evaluated at r_i+1/2. The integration then reduces to a matrix product between  RG_m,jⁱ and  $\delta\sigma_m^i$. Our numerical evaluation of the  RG_m,jⁱ is explained in B.2. Setting  r=r_j+1/2 in the Fredholm equation, the latter is approximated by:

$$% \omega~\delta\sigma_m^j = (m<\Omega>+m~F(r_{j+1/2})~RG_{m,j}^j) ~\delta\sigma_m^j + m~F(r_{j+1/2})~\sum_{i\ne j} RG_{m,j}^i ~\delta\sigma_m^i.$$ (C.4)

We find the eigenvalues and eigenvectors of the matrix A, defined by:

$$% A_{ij} = m <\Omega> \delta_{ij} + m~F(r_{j+1/2})~RG_{m,j}^i.$$ (C.5)

The eigenvectors are defined up to a multiplicative factor. We choose this factor by noting that, according to the wave kinetic Eq. (32), the squared norm N_m²(t) of an eigenvector must evolve according to:

$$% N_m^2(t)=N_m^2(0)~{\rm e}^{\int_0^t 2~\gamma_m(t'){\rm d}t'}.$$ (C.6)

The norm of any initial perturbation  $\delta\sigma_m(r,t=0)$ is known. At a later time, the norm of the eigenvector is chosen in accordance with Eq. (C.6) so that  $\vert\delta\sigma_m(r,t)\vert^2=N_m^2(t)$.