Contents

• Radio pulsar drifting sub-pulses and diocotron instability
• 1 Introduction
• 2 Pulsar magnetosphere
• 3 Diocotron instability
• 4 Numerical results
• 5 Pulsar application
• 6 Discussion and conclusion
• References
• 7 Online Material
• Appendix A: Full form of the perturbed fluid-Maxwell's equation

A&A 445, 779-794 (2006)
DOI: 10.1051/0004-6361:20053040

Radio pulsar drifting sub-pulses and diocotron instability

P. K. Fung¹ - D. Khechinashvili^2,3 - J. Kuijpers²

1 - Astronomical Institute, Utrecht University, PO Box 80000, 3508 TA Utrecht, The Netherlands
2 - Department of Astrophysics, IMAPP, Radboud University Nijmegen, PO Box 9010, 6500 GL Nijmegen, The Netherlands
3 - Centre for Plasma Astrophysics, K. U. Leuven, Celestijnenlaan 200B, 3001 Leuven, Belgium; on leave for Institute of Astronomy, University of Zielena Góra, Lubuska 2, 65-265, Zielena Góra, Poland / Astrophysical Observatory, Al. Kazbegi ave. 2a, Tbilisi 0160, Georgia

Received 10 March 2005 / Accepted 12 July 2005

Abstract
The potential role of a diocotron instability in causing drifting sub-pulses in radio pulsar emission is investigated for aligned magnetic rotators. It is assumed that the out-flowing plasma above a pulsar polar cap consists of an initially axially symmetric, hollow beam of relativistic electron positron pair plasma which carries an electric charge as well as a current. The occurrence of instability depends on shear in the angular velocity distribution of the beam as a function of axial distance. Instability occurs under typical pulsar conditions at mode numbers $\leq$40. It destroys the symmetry of the equilibrium configuration and leads to a carousel of density columns which rotates at fixed angular pattern speed. The process is applied to two pulsars with observed carousels of drifting sub-pulses, and the diocotron instability at corresponding mode number and axial distance is used as a diagnostic for the charge and current density of the polar flow.

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

1 Introduction

Although radio pulsars have been observed for more than thirty-five years, our understanding of the emission process remains poor. The radiation is observed to have a high brightness temperature (Melrose 1995) and is believed to be emitted tangentially to open magnetic field lines above the neutron star's polar caps, somewhere between the stellar surface and the light cylinder. As the star rotates, the magnetic pole sweeps past our line of sight and we receive a radio pulse. The temporal structure of the observed pulse corresponds to the spatiotemporal emission along the strip cut out by our line of sight. The radio emission is probably generated by relativistic particles streaming out along the open field lines, but the precise physical conditions under which the radiation is produced, and the precise location of the emission region above the rotating neutron star, remain obscure. The most remarkable feature of the pulsar is the stability of the pulse profile when averaged over about a hundred subsequent pulses, whereas individual pulses are complex and highly variable. An individual pulse is built out of sub-pulses and sometimes micro-pulses of ${\sim} 100~\mu$s duration. In a number of "ordinary'' (i.e., excluding millisecond) pulsars drifting sub-pulses have been observed which show a remarkable drift stability. It now appears that the occurrence of drifting sub-pulses is not a rare phenomenon but can be detected in half of the radio pulsars (Weltevrede et al. 2005; Ashworth 1981). Drifting sub-pulses are sub-pulses that drift in longitude from pulse to pulse at a nearly constant rate. For some pulsars, a "carousel'' model can be constructed, which consists of emission columns arranged along a circle centered at the stellar magnetic axis, which sub-rotate with respect to the star (Deshpande & Rankin 1999,2001). Understanding drifting sub-pulses is probably of critical importance for our understanding of the pulsar electrodynamics and pulsar emission mechanism.

To date, several models have been proposed to explain this phenomenon. In the model by Ruderman & Sutherland (1975, hereafter RS75), drifting sub-pulses are interpreted in terms of spark-associated plasma columns rotating around the magnetic axis due to the ${\vec E} \times {\vec B}$-drift, where ${\vec B}$ is the background magnetic field and ${\vec E}$ is the electric field in the frame corotating with the star. The assumption of the infinitely large binding energy of ions in the stellar surface, crucial for the spark model, has been severely criticized by several authors (see, e.g., Neuhauser et al. 1987,1986). Moreover, there is strong disagreement between the observed drift speed and the one predicted by the spark model (as shown by van Leeuwen et al. 2003, in the case of PSR B0809+74). In the model by Gil et al. (2003) the vacuum gap is partially screened due to thermionic emission of ions/electrons, and sparks still occur but at a smaller drift rate.

Another model (Clemens & Rosen 2004) ascribes the drifting sub-pulses to non-radial pulsations of the neutron star. A difficulty with this model is the exclusive selection of one particular high-order spherical eigenmode (with a spherical harmonic numbers in between $l \approx 500 {-} 700$ and m=0) at the stellar surface. Also, there is no clear link of the internal oscillations of the star to the processes leading to generation of radio emission higher above the stellar surface. Alternatively, Kazbegi et al. (1991) explain the sub-pulse drift in terms of modulation of the emission region by large-scale "drift waves'', which alter the local curvature of the magnetic field lines and thus affect strongly the occurrence of the Cherenkov-drift resonance - a main mechanism of generation of pulsar radio emission, according to Kazbegi et al. (1991). A recent revised version of this model, by Gogoberidze et al. (2005), favours modulation either through modification of the distribution function of particles by the drift wave electric field, or through changes in the emission angle. The stability of a distinct drift pattern is achieved by the accumulation of the drift modes in a specific azimuthal eigenmode, resulting from their non-linear evolution. The model requires relatively weak magnetic field values, which prevail at distances corresponding to a significant fraction of the light cylinder radius while observations suggest that the emission processes take place at much closer distances to the pulsar surface, typically 50-100 stellar radii (see, e.g., Kijak & Gil 2003). Finally, a recent empirical model constructed by Wright (2003) attributes the formation of a drifting regular pattern of emission nodes to constructive interaction between electron and positron beams traveling, respectively, up and down in the magnetosphere, between inner and outer acceleration gaps. The model, though, does not elaborate on the detailed electrodynamics.

Here, the proposed explanation of drifting sub-pulses is the occurrence of a diocotron instability (also known as "rotational shear'' or "slipping stream'' instability) in the pair plasma on the open pulsar field lines. The diocotron instability is an instability of a sheared flow, much like the familiar Kelvin-Helmholtz instability in a quasi-neutral plasma, now, however, in a non-neutral plasma, which is why it is sometimes also called the "electrostatic'' Kelvin-Helmholtz instability (Hasegawa 1975). Specifically, the diocotron instability occurs in a charged plasma beam along a uniform background magnetic field and a vanishing background electric field component parallel to the background magnetic field (e.g., Davidson 1974, Sect. 2.10). The electric self-field of the charged beam would lead to the familiar ${\vec E} \times {\vec B}$ drift around the cylinder axis which is now, however, modified by the effect of a magnetic self-field from the beam which, in general, carries an electric current as well as a charge. When the resulting differential rotation of the beam around the cylinder axis satisfies a certain condition, a non-axially symmetric perturbation in density, velocity and electric potential can grow over the shear layer. To linear order, the unstable surface modes are characterized by a radial eigenfunction, an azimuthal mode number l, an axial wave-number $k_\parallel$, and an angular frequency $\omega$, whose imaginary part gives the growth rate ${\rm Im}~{\omega} = \Gamma$ and whose real part is related to the angular pattern speed through $\omega_{\rm pat} = {\rm Re}~{\omega}/l$. It is this angular pattern speed which we identify with the rotation speed of the drifting sub-pulses when viewed from the lab frame.

In this paper we model the plasma outflow above the pair-creation front as an infinitely long, non-neutral annular beam composed of electrons and positrons propagating with relativistic speeds along the pulsar's magnetic field. Stability analysis reveals a large domain of realistic parameters, where the - linear stage of - diocotron instability may indeed occur under pulsar magnetospheric conditions. We assume that l and $\omega_{\rm pat}$ of the unstable mode can be linked to the observed number of sub-beams and pattern circulation time $\tau_{\rm circ} = \Omega_\ast/(\Omega_\ast - \omega_{\rm pat})$, where $\Omega_\ast = 2 \pi/P_1$ is the pulsar rotation frequency and $\tau_{\rm circ}$ is measured in units of pulsar periods P₁. Using these two assumptions, we then apply our model to the two pulsars with the best studied carousel patterns, viz. PSRs B0943+10 and B0826-34, and demonstrate that the the model is able to account in a natural way for the basic properties of the drift patterns observed in these pulsars. Furthermore, the model can explain drift direction reversals as in PSR B0826-34. For the relativistic beam to nearly corotate (only slightly sub-rotate) with the star, the model requires the charge difference in the beam to exceed the Goldreich-Julian charge density (Goldreich & Julian 1969; see also Eq. (1) below) many times, which is no problem for a typical pair plasma.

A similar instability, namely an electrostatic shear flow instability, has been studied before in the context of pulsar magnetospheres by Arons & Smith (1979). Although the authors state that their basic equations contain also the physics of other instabilities, including the diocotron instability, they do not treat the latter process, but rather concentrate on implications of the shear in parallel relativistic momentum. The modes excited due to such a shear are space charge waves, contrasting with the diocotron surface modes, treated in the present paper and leading to azimuthal fragmentation of the flow. It is also important to notice that Arons & Smith (1979) assume a small perturbation wavelength along the external magnetic field, whereas our approximation is just the opposite to that (i.e., we have $k_\parallel=0$). In addition, Arons & Smith (1979) only found the electrostatic shear flow instability to appear at enormously large azimuthal wave-numbers (l > 10⁷). We find the same to be true for the diocotron instability if one assumes that the flow is charge-separated, as in Arons & Smith (1979). However, here we will show that the diocotron instability can operate in a relativistic beam under pulsar conditions if the density distribution is not charge-separated, but instead is formed by a pair plasma.

In Sect. 2 we provide the physical background for the inner magnetosphere and describe a beam model with a minimum amount of detail but which does contain the essence of the various existing theories on pulsar electrodynamics. The diocotron instability of the beam is then investigated in a fluid approach and the analytical results are presented in Sect. 3 while the derivation is given in the Appendix. Numerical results on the instability, the pattern speed, and the growth rate for various choices in parameter space are presented in Sect. 4. The implications for pulsars and the observed drifting sub-pulses are given in Sect. 5. We end with a discussion and conclusions in Sect. 6. We use SI units throughout the paper.

   
2 Pulsar magnetosphere

In order to study the physics of the diocotron instability and construct a correct model of the non-neutral pair plasma beam, one should first understand the electrodynamics of the pulsar magnetosphere. Here, after brief review of existing magnetospheric theories, we discuss the nature of the "rotational'' drift in pulsars, in case of both aligned and inclined pulsars. As we will show by simple reasoning, the assumption of independence of the current and charge densities from each other is essential for "near-corotation'' of the open magnetic flux tube with the pulsar.

2.1 Pulsar electrodynamics

Most pulsar theories have in common the existence of a corotating closed magnetosphere and two bundles of "open'' field lines defining the polar caps on the stellar surface and passing through the light cylinder. Everywhere in the corotating magnetosphere the plasma is charge separated - it consists of one kind of charge only - and the charge density is given by the Goldreich-Julian density

$$n_{\rm GJ}({\vec r}) \equiv - \frac{2 \epsilon_0}{e} \frac{{\... ...x - 2 \epsilon_0 {\bf\Omega}_\ast \cdot {\vec B}_0({\vec r})/e$$ (1)

where e>0 is the absolute value of the electron charge, $\mbox{\boldmath$\Omega$ }_\ast$ is the stellar rotational vector, ${\vec B}_0 ({\vec r})$ is the background magnetic field at position ${\vec r}$and c is the speed of light. Theories largely differ, however, on the plasma and electric field distribution above the polar caps (and of course on the precise description of the relativistic pulsar wind outside the light cylinder which is of no direct concern in this paper). Still, most of them agree on a number of qualitative aspects: the existence of an electric potential drop as one goes out from the star along the open field lines, the creation of a relativistic (Lorentz factor $\gamma$ of a few hundred) pair plasma on part of the open bundles (often still accompanied by the remains of an extremely relativistic primary beam, $\gamma \sim 10^6 {-} 10^7$), and the existence of a return current along or in a sheath around the boundary of the closed magnetosphere. A main difference between existing models is in the magnitude of the potential drop available for acceleration of the primary beam. For comparison, the electric field parallel to the background magnetic field at an altitude $r_{pc}\equiv r_\ast (\Omega_\ast r_\ast /c)^{0.5}$, $r_\ast$ denotes the stellar radius, above the magnetic pole can be conveniently written as

$$E_\parallel(r_\ast+r_{pc}) \approx - p c B_\ast \left(\frac{\Omega_\ast r_\ast}{c} \right)^q,$$ (2)

where $B_\ast$ is now the surface magnetic field strength at the magnetic pole, {p, q} are respectively {1, 1} for the aligned vacuum model (Deutsch 1955), {-2, 1.5} for RS75 which is derived for an anti-aligned rotator upon assumption of infinite binding energy for the ions, {0.5, 2.5} for Arons (1983) which has freely streaming electrons from the stellar surface, and {1, 2} for the general relativistic model of Muslimov & Tsygan (1992). The other main difference between models is in the field lines on which pair creation does or does not take place. For example, Arons (1983) and Mestel et al. (1985) obtain different results on the locations of pair creation because of a different assumption about the initial speed of up-flowing primary electrons from the star. The geometry of the various models is summarized schematically in Fig. 1.

Figure 1: Sketch of various models for pair creation (hatched) on the open field lines ${\vec B}$ in the pulsar magnetosphere with level curves of electrostatic potential $\mathcal{V}$ and corresponding electric field ${\vec E}$: a) axially symmetric vacuum model à la Deutsch (1955); b) RS75, anti-aligned pulsar and infinite ion binding energy; c) Arons (1983); d) Muslimov & Tsygan (1992), general relativistic model. The magnetic moment ${\vec m}$ is oriented along the vertical direction, and the rotation vector is ${\bf\Omega }$.

2.2 General characteristics of the pulsar pair plasma

2.2.1 The electric circuit

In the above models the rotating magnetized star acts as a voltage source which drives an electric circuit in the out-streaming plasma flow which is created by the star itself! To guide the formulation of our problem we will use the simplified picture for the electric current as proposed by Goldreich & Julian (1969). In a steady state the voltage generator will draw a current from the star. In the aligned rotator the current comes in through a cylinder around the magnetic axis, is then diverted through the star across the polar cap, where it brakes the stellar rotation by the Lorentz force, goes outward again in a cylindrical sheath, and closes by crossing the magnetic field lines far away from the star, thereby accelerating the stellar wind and depositing stellar angular momentum. Actually the incoming current consists of an outgoing electron/positron plasma, and it is important to realize that - in contrast to a charge-separated flow - the current is not directly coupled to the charge density since not only convection of net charge but also relative drift between electrons and positrons creates electric current. Likewise, the outgoing current again consists of outgoing plasma, either in the form of positive ions from the star and/or an electron/positron plasma (Fig. 2).

The field lines on the magnetic flux surface separating the incoming and the outgoing current are called the critical field lines by Goldreich & Julian (1969). The critical field lines are at the same electric potential as the interstellar medium. A simple estimate of their position on the polar cap can be obtained if one assumes that all over the polar cap the current density is just the Goldreich-Julian density at the stellar surface multiplied with - or + the speed of light depending on whether the current comes in or goes out. The feet of the critical field lines are then positioned at an angle $\theta_{\rm c}$ as seen from the stellar center given by $\theta_{\rm c}=\theta_{pc}/\sqrt{2}= (\Omega_\ast r_\ast /2c)^{0.5}$ which differs only by a factor 0.96 from the feet of the field lines which pass through the cross-section of both the light cylinder and the null lines, and which are located at an angle $\theta_{\rm n}=\theta_{pc}(2/3)^{1.5}$ at the stellar surface. It is interesting to note that the latter field lines define the outer emission cone in Wright (2003).

2.2.2 Rotational drift in oblique pulsars

To obtain a clear picture of the "rotational'' drift of plasma on the open field lines in the general case of an oblique magnetic rotator, it is easiest to split the electric field in the lab frame into two parts: a first part which would lead to corotation of an otherwise static quasi-neutral plasma, and a remaining part, which in view of the assumed steady rotation can be written in terms of an electrostatic potential $\Psi$:

$${\vec E} = -(\mbox{\boldmath$\Omega$ }_\ast \times {\vec r}) \times {\vec B} - \mbox{\boldmath$\nabla$ }\Psi.$$ (3)

Note that such a decomposition can always be made, whether the pulsar is aligned or oblique. Since we assume that both the star and the closed part of the magnetosphere are in solid body rotation at the rate $\bf\Omega_\ast$, by implication both the boundary of the open field bundle and the polar cap surface are at a constant potential.

Figure 2: The current system I in the axially symmetric case; magnetic field ${\vec B}$, magnetic moment ${\vec m}$, rotation vector ${\bf\Omega }$, null line n, light cylinder lc, critical field line c separating incoming from outgoing current.

As a result electric fields deriving from arbitrary $\Psi$ distributions only cause the plasma on the open bundle to circulate inside the open bundle (see Fig. 3).

Figure 3: Level curves of the electrostatic potential $\Psi$ on the open field lines in the oblique case; magnetic field ${\vec B}$, magnetic moment ${\vec m}$, rotation vector $\mbox{\boldmath$\Omega$ }$. The plasma circulation is on surfaces of constant $\Psi$.

As long as we stay well inside the light-cylinder we can approximate the open flux bundle as a cylinder which is invariant in the vertical direction which we call ${\vec s}$. Each charged particle satisfies force equilibrium which - neglecting centrifugal forces - is approximately given by ${\vec E} + {\vec v}_\alpha \times {\vec B}=0$. Then, substituting ${\vec E}$ from Eq. (3), it follows that the particle motion is described - independent of species $\alpha$ - by

 

$${\vec v} = \mbox{\boldmath$\Omega$ }_\ast \times {\vec r} + {\vec u},$$ (4)

$${\vec u} = - \frac{\mbox{\boldmath$\nabla$ }\Psi \times {\vec B}}{B^2}\cdot$$ (5)

Because of the assumption of invariance in the ${\vec s}$-direction there is no electric field along the vertical, and it follows that ${\vec u}$ is everywhere tangent to the local magnetic flux surface which moreover coincides with a surface of constant electrostatic potential $\Psi$!

Figure 4: Axially symmetric case: Sketch of the level curves of the electrostatic potential $\Psi$ in the frame corotating with the star on the left, and of the potential V in the lab frame on the right. In the middle are the level curves due to pure corotation with the star. Darker level of grey corresponds to a higher potential value. Hollow arrows mark the corresponding electric field E, with the arrow size being proportional to the field strength. The sum of the potentials on the left equals to the potential on the right. Magnetic field B (thin arrows), magnetic moment m and rotation vector $\Omega$ are displayed.

2.2.3 Drift of relativistic beam in aligned case

In the rest of the paper, we will study axially symmetric systems and aligned magnetic dipolar rotators, and adopt cylindrical coordinates $\{R, \phi, z\}$. Faraday's law in a steady state implies that the electric field in the lab frame can be entirely written in terms of an electrostatic potential $\mathcal{V}$. Note that we do not split off a part due to stellar rotation as in Eq. (3). As a prelude to the exact derivation in the next section, it is instructive to discuss some general properties of the potential $\mathcal{V}$. The potential $\mathcal{V}$ in the lab frame is related to the potential $\Psi$ in the stellar frame by adding the potential due to pure corotation as is shown in Fig. 4. As a result, the electric field direction at large altitude changes from a direction away from the magnetic axis to a direction towards the magnetic axis when going from stellar to lab frame. Obviously, an electric potential drop across the magnetic field is always associated with a drop of the same magnitude along the field at the foot-points for a steady potential. This applies both to the $\Psi$ and the $\mathcal{V}$ potentials in the aligned, axially symmetric case. A simple relation exists between the perpendicular field higher up and parallel electric fields below. Let us define a potential difference $\mathcal{V}_{ji}$ between points i and j by $\mathcal{V}_{ji}= \int_i^j {\vec{E}} \cdot {\rm d}{\vec{s}}$ then it is immediately obvious from the loop 1,2,3,4,1 in Fig. 4 that the perpendicular electric field high up is given by

$$E_\perp (R)= - \frac{\partial }{\partial R}{\mathcal{V}_\parallel(R)}$$ (6)

where $\mathcal{V}_\parallel(R) = \int_1^2 E{\rm d}l$ is the integrated electric field along a magnetic field line from the stellar surface to point 2. In regions of pair creation and above, we will assume that the plasma is ideal so that the magnetic field can be considered to be "frozen-in'' into the plasma. In this domain, the electric field has no component along the magnetic field, and the situation is as sketched in Fig. 4. The conventional ${\vec E} \times {\vec B}$ motion of plasma and of frozen-in field across the magnetic field at large altitudes is now matched by slippage of the upper part of the flux tube over the lower part which is frozen-in and corotates with the star at the foot-points. The continuous rupture of field lines is made possible by parallel electric fields. This continuous reconnection of field lines is called generalized magnetic reconnection as it is always associated with dissipation of magnetic - i.e. electric circuit - energy.

The field in the rupture zone near the foot-points is of course not frozen into the plasma. Far above the acceleration region, the situation becomes practically independent of altitude and the electric field is entirely due to space charges inside the beam and not to any surface-charges as is the case near the stellar surface. In the RS75 model, such would apply already at a distance of a few polar cap radii above the surface while in the general relativistic model of Muslimov & Tsygan (1992) this would apply above a few stellar radii. Since we are interested in the drift motion of a relativistic plasma beam - which creates its own transverse magnetic field $B_\phi^{\rm S}$ - the drift motion of the plasma as obtained from the approximate force equilibrium condition ${\vec E} + {\vec v} \times {\vec B}=0$ is given by:

$$v_{\alpha\phi}^{0}= -\frac{E_R^0}{B_0} + v_{\alphaz}^{0} \frac{B_\phi^{\rm S}}{B_0},$$ (7)

and is independent of species $\alpha$ for the same outflow velocity. Here the drift is with respect to the lab frame and positive drift is in the sense of stellar rotation. For an aligned rotator, E_R<0 and $B_\phi^{\rm S}<0$, so that the first term on the r.h.s. is positive while the second is negative. It follows that the same perpendicular electric field induces a smaller rotational drift in a charged and out-flowing beam than in a static plasma. We stress that the rotational drift we are talking about is with respect to the lab frame, that is absence of drift means complete lack of corotation. For a pulsar observer who translates these results into the comoving star frame this would mean counter rotation at the stellar rotation rate. In particular, for a relativistic beam of Lorentz factor $\gamma$consisting of one kind of particle only, the repulsive Coulomb electric self-field is largely - by a factor $\gamma^{-2}$ - reduced by the pinching magnetic self-field, a result well-known from lab plasma physics. Note, however, that in a relativistic beam consisting of both positive and negative charges, charge density and current density do not any more directly depend on each other because of the additional freedom of velocities of the two species. Below we will see that this is important for radio pulsars.

As to whether the magnetic field is strong enough to enforce corotation of any plasma with the stellar rotation, the answer depends on the global distribution of electric fields. If there is no potential drop between the stellar surface and a relativistic plasma beam higher up then in a Poynting-flux dominated situation, as is the case for the pulsar, the beam would corotate with the star. The required relatively strong electric field would simply be given by a strong super Goldreich-Julian density. But in a realistic pulsar, the relativistic nature of the beam is caused by a parallel voltage drop that allows for rupture of the magnetic field lines below the beam and may cause a strong sub-rotation of the beam with respect to the stellar rotation.

   
3 Diocotron instability

To study the diocotron instability, the relativistic macroscopic fluid description is used. Each species $\alpha$ is described by the coupled fluid-Maxwell's equations:

      

$$\frac{\partial n_{\alpha}^{}}{\partial t} + \mbox{\boldmath$\nabla$ }\cdot(n_{\alpha}^{}{\vec v}_{\alpha}^{}) = 0$$ (8)

$$\frac{\partial }{\partial t}{\vec P}_{\alpha}+ ({\vec v}_{\alpha}... ...c P}_{\alpha}= e_{\alpha}({\vec E}^{}+ {\vec v}_{\alpha}^{} \times {\vec B}^{})$$ (9)

$$\mbox{\boldmath$\nabla$ }\times {\vec E}^{}= -\frac{\partial }{\partial t}{\vec B}^{}$$ (10)

$$\mbox{\boldmath$\nabla$ }\times {\vec B}^{}= \mu_0{\vec J} + \frac{1}{c^2}\frac{\partial }{\partial t}{\vec E}^{}$$ (11)

$$\mbox{\boldmath$\nabla$ }\cdot{\vec E}^{}=\frac{\rho}{\epsilon_0}$$ (12)

$$\mbox{\boldmath$\nabla$ }\cdot{\vec B}^{}= 0$$ (13)

where ${\vec P}_{\alpha}= m_{\alpha}{\vec v}_{\alpha}^{} (1 - {\vec v}_{\alpha}^{2}/c^2)^{-1/2}$ is the relativistic particle momentum. For simplicity, the dependence on $(t, {\vec x})$ is omitted. In Maxwell's Eqs. (11) and (12), ${\vec J}$ and $\rho$are the current density and the charge density, given by:

$${\vec J} = \sum_{\alpha}e_{\alpha}n_{\alpha}^{} {\vec v}_{\alpha}^{}$$ (14)

$$\rho = \sum_{\alpha}e_{\alpha}n_{\alpha}^{}.$$ (15)

Finite temperature effects are not taken into account in this description.

Figure 5: Top: model for the hollow beam. The density difference between positrons and electrons is $\Delta n = n_{+} - n_{-} = h~n_{\rm GJ}$ in region II, and zero in regions I and III. At R > R3, the plasma is assumed to be in solid-body rotation with the star. Bottom: the axial velocity of the beam particles across the beam is taken to be constant.

The stability of the equilibrium under electrostatic perturbations is investigated where each quantity is perturbed as follows:

$$\psi(t, {\vec x}) = \psi^0({\vec x}) + \psi^1(t, {\vec x})$$ (16)

where superscripts "0'' and "1'' denote the equilibrium and perturbed quantity respectively.

The surface waves are the solution of the linearized fluid-Maxwell equations. Before we analyze the equilibrium situation and the solution of the perturbed system, we first describe the geometry of the charged beam, because the assumptions we make enable us to simplify the set of equations.

   
3.1 Beam model: hollow, nearly hollow and return current

The outgoing plasma above the pair-creation zone is modelled as a differentially rotating cylinder of out-flowing electrons and positrons (denoted with "-''and "+'' respectively). The background magnetic field is assumed to be uniform, ${\vec B}_0 ({\vec r})= B_0 ~\hat{z}$. The plasma is assumed to be cold, and its quantities independent of the z-coordinate. Also, the equilibrium radial density  $n_{\pm}^{0}(R)$, and the equilibrium velocity components  $v_{\pmz}^{0}(R),v_{\pm\phi}^{0}(R)$ are all assumed to be axially symmetric about the magnetic axis. Finally, we assume that electric fields are perpendicular to the magnetic field.

We consider three possible geometries: a hollow beam, a hollow beam including a core component, and a hollow beam in which the current reverses sign.

As a first simple case, we model the beam in this paper as a hollow beam of electrons and positrons with nearly equal axial velocities (Fig. 5). The densities of electrons and positrons are assumed to be constant in different regions of the beam. Then, only the density difference $\Delta n(R) = n_{+}^{}(R) - n_{-}^{}(R)$ is relevant and is taken to be different only in region II. We normalize $\Delta \tilde{n} = \Delta n /~n_{\rm GJ}$ where $n_{\rm GJ} \equiv -2 \epsilon_0 { \Omega}_\ast {B_0}/e$ is not the GJ-density Eq. (1) but an abbreviation involving the absolute value of the local magnetic field strength B₀:

$$\Delta\tilde{n}(R) = \left\{\begin{array}{ll} 0 &\quad 0< R< R_1 ... ...region II)} \\ 0 &\quad R_2 < R < R_3\:\mbox{(region III)}. \end{array}\right.$$ (17)

The closed magnetosphere occupies the region $R \geq R_3$, has a Goldreich-Julian density everywhere, and corotates with the star (this also applies for the resting two cases).

For the motion of the two species, we take $v_{+z} \approx v_{-z} = \beta^0_z c$ independent of R. Additionally, particles are only relativistic in the axial direction, so that $\gamma \approx (1 - \beta_z^{02})^{-1/2}$. We express the axial current density of the beam formally as

$$J_z (R) = f(R)\beta^0_z(R)~e~n_{\rm GJ}~c.$$ (18)

At the polar cap, the density at which the particles stream out is approximately $n_{\rm GJ}$, therefore, we expect $0 \leq f(R) \leq h$ if the current occupies the entire open flux tube, but of course can be much larger if the current is confined to a thin hollow cylinder.

Two additional configurations are examined: the influence of a nonzero core and when part of the beam moves in the opposite direction (return current).

In Fig. 6 the model including a nonzero core is presented. The charge density difference in region I is $h_{\rm I}$. We expect that in this region no pair creation occurs because of the large radius of curvature of the magnetic field, and, therefore, that the beam is charge-separated. The axial velocity remains unchanged, i.e. both in region I and in region II, approximately the same axial velocity applies to all particles.

As for the return current, the axial velocity of the hollow beam changes sign at the radial location

 

R_p = R₁ + p (R₁ - R₂) (19)

which separates region IIa from IIb (see Fig. 7). Note that when the net ring current is zero, we have

 

R_p² = (R₂² + R₁²)/2. (20)

Figure 6: Similar to Fig. 5; now, a core component is included. The dimensionless charge difference in Region I is $h_{\rm I}$.

Figure 7: Similar to Fig. 5, but now, part of the beam moves in the opposite direction (Rp is defined by Eq. (19)).

  
3.2 Equilibrium

We will proceed by analyzing the equilibrium situation. Due to the assumptions of z-invariance and axial symmetry, the equilibrium electric field is determined from Poisson's equation (Eq. (12)) and is given by

 

$$\frac{1}{\epsilon_0 R}\sum_{\alpha}e_{\alpha}\int_0^R \mbox{d}R'~R'n_{\alpha}^{0}(R')$$                       E⁰_r(R) =

$$\frac{e n_{\rm GJ}}{\epsilon_0 R}\int_0^R\mbox{d}R'~\Delta\tilde{n}(R')$$ = (21)

where $\alpha$ denotes the species of the charge in the beam. We use the definition of $\Delta\tilde{n}$ to obtain the latter equality. The equilibrium magnetic field is given by, ${\vec B}^{0} = B_0 \hat{z} + B^{\rm S}_\phi(R)\hat{\phi} + B^{\rm S}_z(R)\hat{z}$, where B₀ is the constant background magnetic field and ${\vec B}^{\rm S}$ is the field generated by the equilibrium axial and azimuthal currents J⁰_z, $J^0_\phi$. Since the background magnetic field is extremely large, we neglect the field $B^{\rm S}_z$ which is generated by $J^0_\phi$. Furthermore, we use Eq. (18) for the equilibrium axial current density, so that:

$$B^{\rm S}_\phi(R) = \frac{1}{\epsilon_0 c R} \int_0^R \mbox{d}R'~R'J_z(R')/c.$$ (22)

Using the assumptions of Sect. 3.1, the zeroth order of both the continuity equation and Maxwell's equations are satisfied, while the azimuthal and axial component of the momentum equation are trivial. Thus, we are left with the radial component of the momentum equation (Eq. (9)), which relates the density distribution $n_{\alpha}^{0}(R)$, the axial velocity $\beta^0_z$, the axial current (in terms of $e n_{\rm GJ} c$) f and the angular speed $\omega_{\alpha}(R) = v_{\alpha\phi}^{0}(R)/R$ of species $\alpha$ to each other:

$$\omega_{\alpha}^2(R) + \epsilon_\alpha\omega_{\alpha}(R) \Ome... ...alpha}R}\left[ E_R^0 - v_{\alphaz}^{0}(R) B_\phi^S\right] = 0$$ (23)

where $\Omega_{B\alpha}\equiv e B_0/\gamma m_{\alpha}$ is the - absolute value of - relativistic gyro-frequency, and $\epsilon_\alpha \equiv \mathop{\mbox{sgn}} e_\alpha$.

If one neglects the first term one immediately derives the earlier result ${\vec E} \times {\vec B}$ (cf. Eq. (7)). Actually, the first term represents the inertial acceleration which can be neglected even for large Lorentz factors because of the strength of the background magnetic field. Defining $Q \equiv \beta_z^{02} f/h$, $Q_{\rm I} \equiv \beta_z^{02} f_{\rm I}/h_{\rm I}$and $\tilde{h}_{\rm I} = h_{\rm I}/h$, the angular rotational velocity $\tilde{\omega}_\alpha(R) \equiv \omega_\alpha(R)/\Omega_\ast$ for the hollow beam (Fig. 5) reads:

 

$$\tilde{\omega}_\alpha(R) = h \left\{\begin{array}{ll} 0 & \mbox{r... ...[2.5mm] (1 - Q)\frac{R_2^2 - R_1^2}{R^2}& \mbox{region III}. \end{array}\right.$$ (24)

For the beam including a core component (Fig. 6), we have

 

$$\tilde{\omega}_\alpha(R) = h\left\{ \begin{array}{ll} \tilde{h}_{... ...- Q_{\rm I}) +\frac{R_2^2 - R_1^2}{R^2}(1 - Q) & \mbox{III}. \end{array}\right.$$ (25)

For our model of the beam including a return current (Fig. 7), we obtain

 

$$\tilde{\omega}_\alpha(R) = h \left\{\begin{array}{ll} 0 &\mbox{re... ...2}\right) + 2 Q\frac{R_p^2 - R_1^2}{R^2}& \mbox{region III}. \end{array}\right.$$ (26)

Note that exact co-rotation only occurs when the beam is solid and $h (1 - Q) \rightarrow 1$ (where h and Q are now the solid beam parameters). This can be seen in either the hollow beam case, then R₁ = 0 and for a beam with h = 1, $\beta_z^{02} f$must be close to zero, or in the second case where we included a core component, where co-rotation is only in Region I. Also, a relativistically out-flowing plasma beam with a Goldreich-Julian current density does not corotate at all with the star unless the charge density is of order $\gamma^2$ times the Goldreich-Julian charge density (e.g. $h = \tilde{h}_{\rm I} = 1$ in Region I of Eq. (25)). Put differently, rotational speeds exceeding the stellar rotation ("super-rotation'') are possible if the charge difference density exceeds the Goldreich-Julian density sufficiently. In general, it is clear from Eqs. (23)-(26) that a realistic pulsar outflow will have a sheared rotational velocity distribution across the beam.

For illustrative purposes, the rotational velocities for the cases (24) to (26) with R₁ = 0.3, R₂ = 0.6 and R₃ = 1 are shown in Fig. 8: hollow beam with Q = 0 (solid), hollow beam with Q= 0.5 (long dash), beam with core ( $h_{\rm I} = 0.2, Q_{\rm I} = 0.1$) and Q = 0.5 (short dash) and beam including return current with Q= 0.5 where R_p is determined by zero ring current (Eq. (20), dotted) and another where p = 0.5(dash-dotted). Note that the presence of a return current in a relativistic beam has the effect of pushing up the rotational velocity toward corotation (compare dotted and long-dashed solutions).

Figure 8: Examples of equilibrium rotational velocity distributions across the beam, where R1,R2 are the inner and outer boundary of the beam. See text for the cases and the beam parameters.

3.3 Linearization

Linearizing the fluid-Maxwell's equations is done by assuming that the perturbations are electrostatic, i.e. $\vert{\vec B}^{1}\vert$ is small compared to $\vert{\vec E}^{1}\vert$ in the momentum equation and Faraday's law. With this, the perturbed electric field can be described as the gradient of a scalar potential:

$${\vec E}^1 = - \mbox{\boldmath$\nabla$ }V^1.$$ (27)

Furthermore, each perturbed quantity is expressed as:

$$\psi^1(t, R, \phi, z) = \sum_{l = -\infty}^\infty \sum_{k_z =... ...y \tilde{\psi}(R) \exp[-{\rm i}(\omega t - l \phi - k_z z)].$$ (28)

In this way, the set of equations with five variablies can be rewritten into one perturbed Poisson's equation with $\tilde{V}$as the only variable (see Appendix A). Because we are considering beams whose motion is only relativistic in the axial direction, terms of order $\beta_{\alpha\phi}^{02}$ and higher are neglected. Also, we assume that the beam density is low and the background magnetic field strength is large, so that $\omega_{p\alpha}^2 \ll \Omega_{B\alpha}^2$. As for the perturbations, we are interested in waves with k_z = 0 and $\vert\omega - l \omega_{\alpha}(R)\vert^2 \ll \Omega_{B\alpha}^2$ (both $\omega$ and $\omega_\alpha$ are of the order of $\Omega_\ast$). With these assumptions: $P \approx -(l/R) \epsilon_\alpha \Omega_{B\alpha}$, $\nu_{\alpha}^2(R) \approx -\Omega_{B\alpha}^2$. Then, the perturbed Poission's equation becomes:

$$\frac{1}{R}\frac{\partial }{\partial R}\left[R\frac{\partial ... ...a}_\alpha(R)}\frac{\partial }{\partial R}{\Delta \tilde{n}(R)}$$ (29)

where $\tilde{\omega} = \omega/\Omega_\ast$.

Now, since we are using step-functions for $\Delta n$, Eq. (29) can be solved analytically. This applies even in the case of a "return current'', because in the derivation of Eq. (29) we have neglected the contributions of k_z and $\beta_{\alpha \phi}^2$.

3.3.1 Eigenvalue equation

The general solution of Eq. (29) for our model, is given by:

 

$$\tphi_{\rm I} = \left(C_1 + \frac{C_2}{R_1^{2l}}\right)R^l,\qquad R < R_1,$$

$$\tphi_{\rm II} = C_1R^l + \frac{C_2}{R^l} , \qquad R_1 < R < R_2,$$ (30)

$$\tphi_{\rm III} = (C_1 R_2^{2l} + C_2) \frac{R_3^{2l} - R^{2l}}{R_3^{2l} - R_2^{2l}}\frac{1}{R^l}, \qquad R_2 < R < R_3.$$

The integration constants C₁ and C₂ are related to each other and can be determined by integration across the discontinuities R = R₁ and R = R₂. This results in an eigenvalue equation of the perturbation for $\tilde{\omega}$:

$$\tilde{\omega}^2 - h G_1 \tilde{\omega} + h^2 G_2 = 0$$ (31)

where

 

$$\frac{l}{h}\Big[\tilde{\omega}_\alpha(R_1) + \tilde{\omega}_\alpha(R_2)\Big] + \Big[\left(R_2/R_3\right)^{2l} - \left(R_1/R_3\right)^{2l}\Big]$$                                             G₁ =

$$- \tilde{h}_{\rm I}\Big[1 - \left(R_1/R_3\right)^{2l}\Big]$$ (32)

 

$$(1 - \tilde{h}_{\rm I})\Bigl[\frac{l}{h}\tilde{\omega}_\alpha(R_2)\big\{1 - (R_1/R_3)^{2l}\big\} \Bigr.$$                                  G₂ =

$$- \Bigl. \big\{1 - \left(R_1/R_2\right)^{2l}\big\}\big\{1 - (R_2/R_3)^{2l}\big\}\Bigr]$$

$$- \frac{l}{h}\tilde{\omega}_\alpha(R_1) + \frac{l^2}{h^2}\tilde{\omega}_\alpha(R_1)\tilde{\omega}_\alpha(R_2).$$ (33)

Solutions for $\tilde{\omega}= {\rm Re}~\tilde{\omega}+ i\:{\rm Im}~\tilde{\omega}$with positive imaginary part correspond to growing perturbations and they occur when

 

G₁² - 4 G₂ < 0. (34)

It can be seen from Eqs. (32) and (33) that an instability is completely determined by just two values: $\tilde{\omega}_\alpha(R_1)/h$ and $\tilde{\omega}_\alpha(R_2)/h$. For example, it then follows from Eqs. (24) and (26) (see also Fig. 8) that the same (in)stability occurs for a non-relativistic hollow beam and a non-relativistic return current. For unstable solutions, the normalized (to $\Omega_\ast$) growth rate $\tilde{\Gamma}$ and the normalized (to $\Omega_\ast$) rotational speed $\tilde{\omega}_{\rm pat}$ of the perturbation (or pattern speed) are given by:

  

$$\tilde{\Gamma} {\rm Im}~\tilde{\omega}= (h/2)\sqrt{4 G - F^2}$$ = (35)

$$\tilde{\omega}_{\rm pat} \frac{{\rm Re}~\tilde{\omega}}{l} = \frac{1}{2}\Big[\tilde{\omega}_\alpha(R_1) + \tilde{\omega}_\alpha(R_2)\Big]$$ =

$$+ \frac{h}{2l}\Big[\big\{(R_2/R_3)^{2l} - (R_1/R_3)^{2l}\big\} - \tilde{h}_{\rm I}\big\{1 - (R_1/R_3)^{2l}\big\}\Big].$$ (36)

From the first term on the r.h.s. of Eq. (36), the pattern speed derived turns out to be approximately halfway between $\tilde{\omega}_\alpha(R_1)$ and $\tilde{\omega}_\alpha(R_2)$, just like for the Kelvin-Helmholtz instability in a quasi-neutral plasma. Physical beams are constrained by: $0 \leq R_1 < R_2 \leq R_3$, so it is obvious that for $h_{\rm I}=0$, Q = 1, we have $\omega_\alpha(R_1)=0$ and $\omega_\alpha(R_2)=0$ from Eq. (24), and therefore G₂ < 0 and the solutions are always stable. Also, for $h_{\rm I}=0$, the fundamental mode l=1 is always stable.

3.3.2 Surface waves

The diocotron instability gives rise to surface waves. This is illustrated in Fig. 9 where the equipotential curves are drawn for $Q = 0, \:R_1 = 0.5, \: R_2 = 0.6, \: R_3 = 1$ and l = 5 at t = 0.

Figure 9: Equipotential for annular beam with Q = 0, R1 = 0.5, R2 = 0.6 and l = 5 at t = 0. Black areas denote negatively valued $\tilde{V}$, and white the positive values.

Obviously, extremes are formed at the two surfaces with five azimuthal modes. Moreover, the surface waves are shifted in phase with respect to each other. To demonstrated this, we use the expression for the perturbed potential at the surfaces Eq. (30), but now, C₁/C₂ is given by:

$$\frac{C_1}{C_2} = \frac{\tilde{\omega}/h - 1}{R_1^{2 l}}\cdot$$ (37)

Because for unstable waves, $\tilde{\omega}$ is complex and the ratio of the amplitude on the surfaces relates as (for simplicity we take t = 0):

$$\frac{\tilde{V}(R_2)}{\tilde{V}(R_1)} = (R_1/R_2)^{-l} + \exp[-i\phi']\frac{(R_1/R_2)^{l} - (R_1/R_2)^{-l}}{\sqrt{G_2}}$$ (38)

where $\tan \phi' = {\rm Im}~\tilde{\omega}/{\rm Re}~\tilde{\omega}$.

Figure 10: Left: unstable diocotron modes in the parameter space of R1/R3 (along the horizontal axis) and R2/R3(along the vertical axis). In levels of grey is plotted $l_{\rm max}$, which is the azimuthal wave number l for which the combination of R1/R3 and R2/R3 results in the largest growth rate (Eq. (35)). The parameter $Q = \beta _z^{02} f/h$ is found above each graph. Right: maximum growth rates $\tilde{\Gamma}/h$ (triangles) and pattern velocities $\tilde{\omega}_{\rm pat}/h$ (squares) as functions of $l_{\rm max}$.

Figure 11: The same as in Fig. 10, but now including a charge-separated core in Region I (see Eq. (25) and Fig. 6). The corresponding value of  $\tilde{h}_{\rm I}$ is specified above each graph, and $Q_{\rm II}=0.4$ in all graphs. The growth rate and the pattern speed are displayed as squares and triangles, respectively.

Figure 12: The same as in the left column of Fig. 10, but now including a return current (see Eq. (26) and Fig. 7). Above each graph are seen the corresponding parameters Q and p, the latter defining Rp, i.e. the radius where the axial velocity of the beam particles changes sign (see Eq. (19)).

Figure 13: Maximum growth rate (squares) and pattern velocities (triangles) corresponding to the case in Fig. 12.

Thus, in general, the two waves are shifted in phase, which is determined by the geometry. Due to the shift, the two surface waves are coupled to each other through their perturbed potential. To be more specific, under the condition determined by Eq. (34), perturbations at the outer surface can influence the perturbations at the inner surface, and vice versa, leading to exponential growth.

It is this system that probably will evolve into a non-linear stage where particles are trapped and vortices form.

   
4 Numerical results

We proceed with the analysis by searching the domain of instability for $l_{\rm max}$: the mode number lwith the largest growth rate as a function of ( R₁/R₃, R₂/R₃). These results for a hollow beam, a beam including a core component and a beam including return current are displayed in greyscale in Figs. 10-13. In the contourgraphs, darker areas represent larger $l_{\rm max}$: for $l_{\rm max} \in [2,10]$ we use different level of greyscale for each integer $l_{\rm max}$ (some of which are depicted in the graphs), whereas for the ranges: [11, 20], [21, 30] and [31, 40], one greyscale is used for each range; black represents $l_{\rm max} > 40$, and such values are apparent in the graphs for $Q \geq 0.6$. Complimentary to the contourgraphs we present the maximum of growth rates $\tilde{\Gamma}/h$ (Eq. (35)) and pattern velocity $\tilde{\omega}_{\rm pat}/h$ (Eq. (36)) found in each $l_{\rm max}$-region.

For a hollow beam, in which Q ranges from 0 to 0.8 in steps of 0.2, $l_{\rm max}$ is presented in Fig. 10. In the case of Q = 0 (i.e. either the beam is nonrelativistic or the beam carries no current), the beam displays instability for all modes. The general trend for all Q's is that thinner beams, i.e. $R_1 \simeq R_2$, are unstable for higher l.

An increase in Q leads to a shift of the $l_{\rm max}$regions in parameter space and a decrease in their area, until, for small  $l_{\rm max}$, instability eventually vanishes. In the case of Q = 0.9, almost any hollow beam is stable and therefore, the corresponding parameterspace is not presented. This is in agreement with the results of Arons & Smith (1979) for a charge-separated, relativistic, particle beam. In that case: fequals h, the first term for both G₁ and G₂ becomes negligible when we assume $\gamma_\alpha \gg 1$, and solutions for $\tilde{\omega}$ are very stable unless $l \gg \gamma^2$ (see Eqs. (33) and (34)).

Figure 11 depicts the influence of adding a charge separated core component by increasing $\tilde{h}_{\rm I}$ in steps from 0.1 to 0.9 (in region II: Q = 0.4). Since it is charge separated, we have $1 - Q_{\rm I} = \gamma^{-2}$. Clearly, when $Q \neq Q_{\rm I}$ the shear is still present (see Eq. (25)) and we observe a diocotron instability.

In Figs. 12 and 13 is shown the effect of a return current on the occurence of instability. In the case where R_p is determined by the condition that the net ring current is zero, the rotational velocities at R₁ and R₂are identical to those of Q=0 in the hollow beam case. Since the instability condition depends only on these velocities (Eqs. (31)-(34)) there is no need to locate the unstable domains in the parameter space again.

   
5 Pulsar application

To demonstrate the applicability of the diocotron instability to radio pulsars, we use the hollow beam model as a simple example. For the above investigated cases of the diocotron instability we made the following assumptions: the rotation axis is aligned with the magnetic axis, the magnetic field is uniform in the axial direction, $k_\parallel=0$, and a core component is absent (or at least there is no spacing between the beam components). Since our model is valid for aligned pulsars, we have chosen two pulsars which show drifting sub-pulses and are nearly aligned: PSR B0943+10 and PSR B0826-34. To apply our results, we need to extract the following information from the observations: the azimuthal mode number l, the pattern frequency  $\omega_{\rm pat}$ in the laboratory frame, the inner and outer radius of the subpulse, R₁, R₂, relative to the cone radius R₃ (which we assume to correspond to the last open field line). R₁ and R₂ are retrieved from the observations of the specific pulsar. The best estimate for R₃ (at the height of emission) is obtained by using the width of the average profile together with $\alpha$ and $\beta$ of that pulsar (see, e.g., Eq. (2) in Rankin 1993).

The data available from the literature is shown in Table 1.

 

 

Table 1: Observationally inferred pulsar data. Note that in Esamdin et al. (2005) no value for $\beta$ is quoted. Also, we infer P₃ from their values of P₂ and using the relation for the drift rate $D (^\circ /P_1) = P_2(^\circ )/P_3(P_1)$, which is a good approximation for this system, as it is nearly aligned. Positive and negative values are obtained due to the definition of drift rate.

PSR	$\alpha (^\circ)$	$\beta (^\circ)$	P1(s)	P3 (P1)	l	$\tau_{\rm circ} (P_1)$	$\tilde{\omega}_{\rm circ}$
B0943+10 (Deshpande & Rankin 2001)	11.58	-4.29	1.09	1.87	20	37	0.973
B0826-34 (Gupta et al. 2004)	1-2	1	1.848	0.50	15	7.5	0.833
B0826-34 (Esamdin et al. 2005)	0.5	-	1.848	{+6, -7}	13	+78, -91	{ 0.99, 1.01}

5.1 PSR B0943+10

Drifting sub-pulses of this pulsar have been studied extensively by Deshpande & Rankin (1999,2001). They argue that the alias problem is solved and they found 20 sub-beams which sub-rotate with respect to the pulsar, and complete one circulation within 37 pulsar periods. This corresponds to an angular frequency in the laboratory frame of $\omega_{\rm pat} = 0.973~\Omega_\ast$. From the average pulse profile (Fig. 1 in DR01) we obtain for the full pulse width $W = 34^\circ$. Together with $\alpha = 11\hbox{$.\!\!^\circ$ }58$ and $\beta = -4\hbox{$.\!\!^\circ$ }29$, we obtain for the cone radius $\rho = 5\hbox{$.\!\!^\circ$ }08$. To apply our model of a hollow beam, we estimate $R_2/R_3 \sim 0.86 {-} 0.93$ (Fig. 10 in DR01). If we assume the sub-pulse shape to be circular, we get $R_1/R_3 \approx 0.65$.

For R₁/R₃ = 0.65 and R₂/R₃ = 0.92, we calculate $\tilde{\Gamma}/h$ in the parameterspace (l, Q) (see Fig. 14). Apparently, growth rate is not at optimum for l = 20. Since we do not have more information on Q, we conclude that for l = 20 to develop, we need $Q \approx 0.8$. By combining this with measured $\tilde{\omega}_{\rm pat} = 0.973$, we find from Eq. (36) that $h = \Delta n/ n_{\rm GJ}\approx 20$.

Since the values for R₁/R₃ and R₂/R₃ are just estimations, it is important to see what the influence is if we change the beam configuration. The results are shown in Figs. 15. The values of Q-parameter, providing instability for l = 20, are listed in Table 2. If we apply the measured pattern velocity $\tilde{\omega}_{\rm pat} = 0.973$ to all five cases shown in Table 2 and plug the corresponding values into Eq. (36), we obtain that for the instability to develop, the charge density difference should be $\Delta n \approx 19.5~ n_{\rm GJ}$.

 

 

Table 2: The results of Q and $\tilde{\Gamma}$ for possible beam configurations for PSR B0943+10 by assuming l = 20 and $\tilde{\omega}_{\rm pat} = 0.973$. The charge density difference is found to be $h \approx 19.5$. Initial guess for R₁/R₃ and R₂/R₃ is taken from Fig. 10 in DR01.

R1/R3	R2/R3	Q	$\tilde{\Gamma}$
0.65	0.92	0.80	$1.8 \times 10^{-2}$
0.60	0.92	0.83	$3.6 \times 10^{-3}$
0.80	0.92	0.60	1.14
0.65	0.86	0.77	$7.2 \times 10^{-2}$
0.65	0.93	0.81	$1.4 \times 10^{-2}$

Figure 14: Growth rates as functions of Q, for different values of azimuthal mode number l (indicated in the plot), for a hollow beam with R1/R3 = 0.65 and R2/R3 = 0.92.

5.2 PSR B0826-34

This pulsar was observed by Manchester et al. (1978), Durdin et al. (1979), Biggs et al. (1985), Reyes et al. (1995) at different frequencies. Already in the first observations it was noticed that the radiation spans the full $360^\circ$. The average pulse profile is divided into four regions, labeled I-IV. At low frequencies, e.g. 408 MHz, the main pulse (III) dominates over a weak interpulse (I) in the average pulse profile, whereas at higher frequencies, e.g. 1374 MHz, the relative intensities switch, and the intensity in I is much higher than in III (this is seen in average pulse profiles at different frequencies, e.g., in Fig. 1 in Esamdin et al. 2005, hereafter, E05). Also, this pulsar is remarkable due to its long nulls, e.g. 70$\%$ of the time at 409 MHz (Durdin et al. 1979). Already in the early observations of PSR B0826-34, the drifting sub-pulses were observed with an unusual apparent drift pattern, namely a sign reversal of the drift rate (Biggs et al. 1985).

Figure 15: Similar to Fig. 14, but now for different combinations of the inner and outer radii of the beam. Top: R1/R3 = 0.6 ( left) and R1/R3 = 0.8 ( right). Bottom: R2/R3 = 0.86 ( left) and R2/R3 = 0.93 ( right).

Recently, two detailed observational studies have been made of this pulsar. Gupta et al. (2004) studied the main pulse (III) with the GMRT at 318 MHz. They claim that the apparent sign reversal is merely an effect of the geometry and aliasing. Our line-of-sight is very near the magnetic axis $\beta = 1^\circ$, and this axis in turn is close to the rotation axis, $\alpha \leq 5^\circ$(probably ${\approx} 1^\circ {-} 2^\circ$). In this way, the observations are still compatible with the Ruderman & Sutherland model. They concluded that the carousel has about 15 sparks (at 318 MHz), rotating with an angular frequency of $\tilde{\omega}_{\rm pat} = 0.867$. Another observation was made with Parkes at 1.374 GHz by E05. They studied the whole pulse profile and used a method called "phase tracking'' where (as the name suggests) the phase of the pulse throughout the whole profile is tracked, which resulted in 13 sub-beams. At this frequency, the interpulse is much stronger than the main pulse, therefore, they came to the conclusion that we observe two nested cones, (see Fig. 10 in their paper, note that both cones have 13 sub-beams). They argue that the drift reversal cannot be accounted for in terms of geometry and aliasing (which would be too much of a coincidence), but that it is a true reversal, implying that the beams are rotating $3.6^\circ$ and $3.2^\circ$ per rotation period in either direction, corresponding to pattern velocities of $\tilde{\omega}_{\rm pat} = 1.01$ and 0.99, respectively. Contrary to the model of Rutherman and Sutherland, the diocotron instability can account for a reverse drift. Therefore, we will illustrate the drift reversal in this context.

Here, we assume that the two observed sub-beam rings are formed by an instability of two nested hollow beams. The treatment of the inner hollow beam is unaffected by the outer beam, since the electric field causing the drift is only determined by the enclosed charge. The outer beam, on the other hand, cannot be described correctly within the context of our model without further assumption, since two different hollow beams would lead to an eigenvalue equation for $\tilde{\omega}$ of 4th degree instead of the quadratic Eq. (31). For simplicity, we therefore assume that the outer beam can be treated independently. Although E05 suggested a polar cap geometry, R₁/R₃ and R₂/R₃ are not well enough determined. So, unlike in PSR B0943+10, we cannot first constrain Q by fixing the hollow beam geometry. Instead, we allow different Q and plot the instability region for l = 13 and Q = 0 - 0.7 (Fig. 16; white regions are stable). Inside the domain we separate - by using greyscales and black - the regions according to equal "angular'' velocities $\tilde{\omega}_{\rm pat}/h$ (black represents $\tilde{\omega}_{\rm pat}/h > 0.08$, and the other two have: $0 < \tilde{\omega}_{\rm pat} \leq 0.04$, $0.04 < \tilde{\omega}_{\rm pat} \leq 0.08$). Note that instead of (R₂/R₃,R₁/R₃), we now use the center of the beam, R_b = (1/2)(R₁ + R₂) and the beam half-width, $\Delta R = (1/2)(R_2 - R_1)$.

To have a reverse drift, our model implies that the charge density differences are changing around a critical value. For instance, if we take Q = 0.7, and for the beam radius of the inner cone $R_{\rm b} = 0.56$ (a rough estimation from Fig. 10 in E05), then for l = 13, the maximum of $\tilde{\Gamma}/h$ occurs at $\Delta R = 0.0995$, which gives us $\tilde{\omega}_{\rm pat}/h = 0.077$. For exact rotation with the star, we need $h = 1/0.077 \approx 13$ (which corresponds to $\tilde{\Gamma} \approx 0.12$). Thus, for slightly larger or smaller $\Delta n$, the pattern speed will be either sub- or super-rotational as compared to the stellar rotation. Also, note that the beam thickness $\Delta R$ increases with the beam radius $R_{\rm b}$ for the same (l,  $\omega_{\rm pat}$), just like E05 suggested for the two rings on the basis of their observations (see their Fig. 10). If we choose, e.g., $R_{\rm b} = 0.77$ for the outer cone, then, for this beam to be unstable for l = 13 and to have the same pattern velocity as the inner beam, we find a maximum of $\tilde{\Gamma}/h$, now corresponding to $\Delta R = 0.13$, which is indeed slightly larger than the beam thickness of the inner cone.

Figure 16: Domain of instability for l=13 as applied to PSR B0826-10. The three regions denote different intervals of $\tilde{\omega}_{\rm pat}/h$ (see legend).

Despite the absence of observationally well-determined values of $R_{\rm b}$ and $\Delta R$, our results demonstrate the onset of instability for the required l-value and pattern speed over a broad range of thin beams.

   
6 Discussion and conclusion

We have investigated the diocotron instability (alternatively called a "slipping stream'' or an "electrostatic Kelvin-Helmholtz'' instability) in the context of drifting sub-pulses observed in radio profiles of a number of pulsars. We modelled an out-flowing relativistic pair plasma as a non-neutral, differentially rotating, infinitely long and azimuthally symmetric cylindrical annular beam, propagating along the pulsar's strong magnetic field. In the limit of long axial wavelengths ( $k_\parallel=0$), treated in this paper, we find that the instability occurs for a broad range of beam geometrical parameters.

We find that for charge separated flow, only the modes with very high azimuthal harmonic number, $l\geq \gamma^2$, are unstable. However, as we have demonstrated, when the flow consists of both positrons and electrons, the charge density and the current density become largely independent of each other, and the beam becomes unstable in the relativistic case already for moderate azimuthal numbers, $l\leq 40$.

We apply our model to two pulsars with sufficiently detailed carousel observations, namely PSRs B0943+10 and B0826-34, identifying the azimuthal number l of the unstable surface modes with the number of sub-beams in a particular pulsar profile. In both cases we have demonstrated that the hollow beam model can account for the observed sub-pulse drift. It is important that since the drift direction is determined by the non-dimensional charge density h (i.e., relative excess of the charge of one sign over another), together with the normalized current density Q, direction reversal of the drift is possible within our model. This is observed in PSR B0826-34. Such an application opens up the possibility to determine details of the polar cap plasma, such as charge and current densities of the particles, within the framework of the diocotron instability model.

However, an exact comparison between the model and the observations is far from being complete. The number of sub-beams observed in a particular pulsar depends sensitively on our line-of-sight. Furthermore, a detailed description of carousels underlying observed drifting sub-pulses is only known for very few pulsars. Empirical models of the emission geometry point to radiation coming from a number of nested cones plus a "core'' component. On the other hand, for more than one cone, the number of discontinuities to match becomes more than two. This results in a higher-order polynomial equation, instead of the quadratic Eq. (31), which is not accessible to analytical treatment but rather requires involvement of numerical methods. This applies, in particular, to our treatment of PSR B0826-10 where E05 suggest that the whole carousel drifts with one speed.

In both applications we find that a strong super-GJ charge density difference is needed for the relativistic beam to nearly corotate with the star. There are no difficulties for the net charge of a pair plasma beam to achieve such densities. This is in contrast with a charge-separated relativistic beam which has a low rotation rate ( $h \approx 1$ and $1 - \beta_z^{02}f$ is small in Eq. (24)). This has been noticed before by Fawley et al. (1977), and results from the approximate balance of the repelling electric radial force by the pinching force of the self-magnetic field. Of course, when a relativistic beam is in approximate corotation with the star there is only a small parallel voltage drop deep below the resonance region, close to the stellar surface, and conversely when the beam is not rotating with the star there exists a large parallel voltage drop at low heights which allows for strong slippage of the magnetospheric field with respect to the star.

As mentioned above, we have limited ourselves to purely perpendicular wave modes, which is a textbook approach to diocotron instability of surface waves at the beam boundary (see, e.g., Davidson 1974). On the other hand, if one includes a finite wave number parallel to the vertical direction (as done by Chen & Uhm 1985, in the case of a relativistic pure electron beam), the analysis becomes more complicated, whereas for the application to drifting sub-pulses we expect only a gradual change in the results by extending the geometry of fluted columns to spiralling structures. To justify the use of invariance in the vertical direction in our treatment, we note that at the pulsar surface a feedback exists between the potential structure at and above the stellar surface, the charge density perturbations in the beam and the density perturbations. Even though the particles themselves stream out at approximately the speed of light, the perturbations in potential are left behind, which in turn determine perturbations in the energy of primary particles. The latter might affect the density (multiplicity) of out-flowing pair plasma, causing filamentation of the flow, and therefore modulating the radio emission pattern generated by this plasma.

In the above analysis we have simplified the problem by assuming that the pulsar is aligned, $\mbox{\boldmath$\Omega$ } \vert\vert {\vec m}$. Since it now appears that the majority of pulsars exhibit the drifting sub-pulse phenomenon, the diocotron instability should be examined in oblique rotators, where the equilibria become time-dependent. Future work on this instability should incorporate this, as well as the other points mentioned above.

Although we have investigated the linear stage of the diocotron instability only, practical applications would correspond to solutions in the non-linear regime. In published work on the non-linear stage of the diocotron instability, the surface waves interact with each other through the perturbed potential, and result in stable vortices (Rosenthal et al. 1987; Davidson et al. 1991; Driscoll & Fine 1990). The diocotron instability modulates and fragments a charged cylindrical beam in the azimuthal direction as a result of initial shear of the plasma angular velocity around the axis. We expect that drifting sub-pulses correspond to such vortices from the non-linear phase of the diocotron instability instead of the linear structures.

The radio emission from pulsars presents more puzzles than the drifting sub-pulse phenomenon, the major one being the mechanism of the radio emission itself. Here we note that Mahajan et al. (1997) have suggested that the diocotron instability can couple non-escaping longitudinal Langmuir waves to escaping electromagnetic waves.

Acknowledgements

We would like to thank Axel Jessner and Jérôme Pétri for stimulating discussions. D.Kh. thanks the Netherlands Research School for Astronomy (NOVA) for support, and the Department of Astrophysics at the Radboud University Nijmegen for hospitality during his stay there.

This work has been funded under a joint research project between the Centre for Plasma Physics and Radiation Technology (CPS) and The Netherlands Research School for Astronomy (NOVA).

References

 
• Arons, J. 1983, ApJ, 266, 215
   
• Arons, J., & Smith, D. F. 1979, ApJ, 229, 728
   
• Ashworth, M. 1981, Ph.D. Thesis, The University of Manchester, Manchester, UK
   
• Biggs, J. D., McCulloch, P. M., Hamilton, P. A., Manchester, R. N., & Lyne, A. G. 1985, MNRAS, 215, 281
   
• Chen, H. C., & Uhm, H. S. 1985, Phys. Rev. A, 32, 1657
   
• Clemens, J. C., & Rosen, R. 2004, ApJ, 609, 340
   
• Davidson, R. C. 1974, Theory of Nonneutral Plasmas (W. A. Benjamin, Inc.)
   
• Davidson, R. C., Chan, H., Chen, C., & Lund, S. 1991, Rev. Mod. Phys., 63, 341
   
• Deshpande, A. A., & Rankin, J. M. 1999, ApJ, 524, 1008
   
• Deshpande, A. A., & Rankin, J. M. 2001, MNRAS, 322, 438
   
• Deutsch, A. J. 1955, Annales d'Astrophysique, 18, 1
   
• Driscoll, C. F., & Fine, K. S. 1990, Phys. Fluids B, 2, 1359
   
• Durdin, J. M., Large, M. I., Little, A. G., et al. 1979, MNRAS, 186, 39P
   
• Esamdin, A., Lyne, A. G., Graham-Smith, F., et al. 2005, MNRAS, 356, 59
   
• Fawley, W. M., Arons, J., & Scharlemann, E. T. 1977, ApJ, 217, 227
   
• Gil, J., Melikidze, G. I., & Geppert, U. 2003, A&A, 407, 315
   
• Gogoberidze, G., Machabeli, G. Z., Melrose, D., & Luo, Q. 2005, ArXiv Astrophysics e-prints
   
• Goldreich, P., & Julian, W. H. 1969, ApJ, 157, 869
   
• Gupta, Y., Gil, J., Kijak, J., & Sendyk, M. 2004, A&A, 426, 229
   
• Hasegawa, A. 1975, Plasma instabilities and nonlinear effects, Physics and Chemistry in Space (Berlin: Springer)
   
• Kazbegi, A. Z., Machabeli, G. Z., & Melikidze, G. I. 1991, Aust. J. Phys., 44, 573
   
• Kijak, J., & Gil, J. 2003, A&A, 397, 969
   
• Mahajan, S. M., Machabeli, G. Z., & Rogava, A. D. 1997, ApJ, 479, L129
   
• Manchester, R. N., Lyne, A. G., Taylor, J. H., et al. 1978, MNRAS, 185, 409
   
• Melrose, D. B. 1995, JA&A, 16, 137
   
• Mestel, L., Robertson, J. A., Wang, Y.-M., & Westfold, K. C. 1985, MNRAS, 217, 443
   
• Muslimov, A. G., & Tsygan, A. I. 1992, MNRAS, 255, 61
   
• Neuhauser, D., Koonin, S. E., & Langanke, K. 1987, Phys. Rev. A, 36, 4163
   
• Neuhauser, D., Langanke, K., & Koonin, S. E. 1986, Phys. Rev. A, 33, 2084
   
• Rankin, J. M. 1993, ApJ, 405, 285
   
• Reyes, F., Aparici, J., & Olmos, F. 1995, A&A, 301, 182
   
• Rosenthal, G., Dimonte, G., & Wong, A. Y. 1987, Phys. Fluids, 30, 3257
   
• Ruderman, M. A., & Sutherland, P. G. 1975, ApJ, 196, 51
   
• van Leeuwen, A. G. J., Stappers, B. W., Ramachandran, R., & Rankin, J. M. 2003, A&A, 399, 223
   
• Weltevrede, P., Edwards, R., & Stappers, B. W. 2005, in prep.
   
• Wright, G. A. E. 2003, MNRAS, 344, 1041
  

 

  
7 Online Material

   
Appendix A: Full form of the perturbed fluid-Maxwell's equation

The coupled fluid-Maxwell's equations (8)-(13) are linearized and electrostatic perturbations are applied. Using the assumption (27), the set of first order equations becomes:

     

$$-iA_\alpha \tilde{n}_\alpha + \frac{1}{R}\frac{\partial }{\partia... ...}^{0}\tilde{v}_{\alpha \phi}}{R} + i k_z n_{\alpha}^{0} \tilde{v}_{\alpha z}= 0$$ (A.1)

$$-i A_\alpha \tilde{v}_{\alpha R} - \left[A_2 + \gamma_{\alpha}^{2... ...e_{\alpha}}{\gamma_{\alpha}^{} m_{\alpha}}\frac{\partial \tilde{V}}{\partial R}$$ (A.2)

$$\left[A_3 + \gamma_{\alpha}^{2}\beta_{\alpha\phi}^{02}\frac{\part... ...pha z} = -\frac{il}{R}\frac{e_{\alpha}}{\gamma_{\alpha}^{} m_{\alpha}}\tilde{V}$$ (A.3)

$$\left[(1 +\gamma_{\alpha}^{2}\beta_{\alphaz}^{2})\frac{\partial v... ...e{v}_{\alpha z} = -ik_z\frac{e_{\alpha}}{\gamma_{\alpha}^{}m_{\alpha}}\tilde{V}$$ (A.4)

$$\frac{1}{R}\frac{\partial }{\partial R}\left(R \frac{\partial }{\... ...- k_z^2\tilde{V}= -\frac{1}{\epsilon_0}\sum_{\alpha}e_{\alpha} \tilde{n}_\alpha$$ (A.5)

where $\gamma_{\alpha}^{} = (1 - \beta_{\alphaz}^{02} - \beta_{\alpha\phi}^{02})^{-1/2}$, $A_2 = 2\omega_{\alpha} + \epsilon_\alpha\Omega_{B\alpha}$, $A_3 = \omega_{\alpha}+ \epsilon_\alpha\Omega_{B\alpha}+ \frac{\partial v_{\alpha\phi}^{}}{\partial R}$ and $\Omega_{B\alpha}= e B_0/\gamma_{\alpha}^{}m_{\alpha}$.

Equations (A.1)-(A.5) result in one equation for $\tilde{V}$:

 

$$\frac{1}{R}\frac{\partial }{\partial R}\left[R\left( 1 - \sum_{\a... ...R} + \frac{\partial v_{\alphaz}^{0}}{\partial R}\right) \right)\right]\tilde{V}$$

$$- k_z^2\left[ 1 - \sum_{\alpha}\frac{\omega_{p\alpha}^2}{A_\alpha... ...rtial \gamma_{\alpha}^{}}{\partial R}}{\nu_{\alpha}^2}\right) \right]\tilde{V}=$$

$$\frac{\tilde{V}}{R} \sum_{\alpha}\frac{1}{A_\alpha}\frac{\partial... ...}{R}\sum_{\alpha}\frac{\omega_{p\alpha}^2}{\nu_{\alpha}^2} \frac{1}{A_\alpha^2}$$

$$\times \left[\beta_{\alpha\phi}^{}\beta_{\alphaz}^{}\left(2 A_\al... ...A_1 + 2 \omega_{\alpha}\beta_{\alpha\phi}^{02}\gamma_{\alpha}^{2}\right)\right]$$ (A.6)

where

$$\omega_{p\alpha}^2= \frac{e^2n_{\alpha}^{0}}{\gamma_{\alpha}^{} m_{\alpha}\epsilon_0}$$ (A.7)

$$\Omega_{B\alpha}= \frac{e B_0}{\gamma_{\alpha}^{}m_{\alpha}}$$ (A.8)

$$A_\alpha = \omega - l \omega_{\alpha}- k_z v_{\alphaz}^{0}$$ (A.9)

$$A_1 = \omega_{\alpha}+ \epsilon_\alpha\Omega_{B\alpha}$$ (A.10)

$$G = -\frac{l}{R}A_2 +\left(\frac{l}{R}\beta_{\alpha\phi}^{02} + k_z \beta_{\alphaz}^{0}\beta_{\alpha\phi}^{0}\right)A_1$$ (A.11)

$$\nu_{\alpha}^2(R) = A_\alpha^2 - A_1 A_2 + \omega_{\alpha}^2\frac... ...az}^{02} - \frac{1}{\gamma_{\alpha}^{2}}\right) + \beta_{\alpha\phi}^{02}A_1^2.$$ (A.12)

For our purpose, Eq. (A.6) reduces to:

 

$$\frac{1}{R}\frac{\partial }{\partial R}\left[R\frac{\partial \tilde{V}}{\partial R}\right] - \frac{l^2}{R^2}\tilde{V} -\tilde{V}\frac{l}{R} \sum_{\alpha} \frac{1}{\left[\omega - l \om... ...rtial R}\left(\frac{\omega_{p\alpha}^2}{\epsilon_\alpha\Omega_{B\alpha}}\right)$$ =

$$\tilde{V}\frac{l}{R}\frac{2 \Omega_\ast}{\omega - l \omega_\alpha(R)}\frac{\partial }{\partial R}{\Delta \tilde{n}(R)}$$ = (A.13)

$$\tilde{V}\frac{l}{R}\frac{2}{\tilde{\omega} - l \tilde{\omega}_\alpha(R)}\frac{\partial }{\partial R}{\Delta \tilde{n}(R)}.$$ = (A.14)