A&A 387, 339-346 (2002)
DOI: 10.1051/0004-6361:20020302

Coupling between electrostatic and electromagnetic behaviour in a cylindrical equal-mass plasma

D. A. Diver1 - A. A. da Costa 2 - E. W. Laing 1


1 - Dept of Physics and Astronomy, University of Glasgow, Glasgow G12 8QQ, Scotland, UK
2 - Centro de Electrodinâmica, Instituto Superior Técnico, 1049-001 Lisboa, Portugal

Received 12 December 2001 / Accepted 27 February 2002

Abstract
The electron-positron plasma present in pulsar magnetospheres may have non-linear electrostatic oscillations, which can couple to electromagnetic waves if the magnetic field present is inhomogeneous. Using a cylindrical geometry the basic equations are derived, and using a simplified linear approach characteristic linear modes are identified and analysed. A possible coupling mechanism is discussed, and possible consequences for pulsar radiation mechanisms are inferred.

Key words: plasmas - radiation mechanisms: non-thermal - pulsars: general


1 Introduction

In a previous paper (da Costa et al. 2001, hereafter referred as Paper I) it was shown that the motion of equal-mass plasma in pulsar magnetospheres required the consideration of non-linear effects in the rest frame of the plasma. Within the context of cold plasmas, the non-linear electron-positron oscillation exhibits a density instability, in which the number density of each species grows sharply at the edges of the confined oscillation. Such a phenomenon is unphysical since ultimately annihilation will result.

There are several possible mechanisms which may quench this instability, such as pressure effects which would inhibit density growth. In this paper an investigation into the coupling of electrostatic and electromagnetic modes has been initiated, in which electrostatic perturbations trigger electromagnetic radiation. Such a coupling is desirable for two reasons: not only would it limit the energy in the oscillation by radiating it away, it could also act as a seat of radiation in a pulsar environment.

This mechanism is consistent with other forms of indirect radiation mechanisms, where primary effects generate secondary electromagnetic effects (Melrose 1993; Asseo & Melikidze 1998; Melikidze et al. 2000; Asseo & Riazuelo 2000). There are however two main differences.

First, this is no longer a pure dynamical process in the laboratory frame of the pulsar, but one in the rest frame of the plasma translated in the laboratory frame, as suggested by Melrose (2000). Due to the quasi-static nature of the phenomena in pulsar magnetospheres where the azimuthal dependence is given by $\psi\equiv\varphi-\omega t$, this is a real lighthouse effect which when sitting outside the speed of light cylinder will have superluminal characteristics (da Costa & Kahn 1985; Ginzburg 1989).

Second, the primary process is not strictly thermal, as seen in Paper I, but it might be a stochastic perturbation induced in the rest frame either directly by the quasi-classical process as described by da Costa & Kahn (1997), or other possible phenomena, as the presence of complex non-thermal plasma density distribution functions which will act to produce non-thermal phenomena similar to the one described by Asseo & Melikidze (1998).

The present paper addresses the basic mathematical formulation of the electron-positron (E-P) cold magnetoplasma connection between electrostatic and electromagnetic modes in Sect. 2, in the rest frame of the initially unperturbed magnetoplasma. It then proceeds in Sect. 3 to formulate the numerical algorithm designed specifically to solve the problem. Results are also presented here, and are discussed in Sect. 4. Finally, the implications for future research are analysed in Sect. 5.

2 Model equations

In Paper I the non-linear oscillations of an unmagnetized E-P plasma were studied. However in order to couple the electrostatic mode to propagating electromagnetic radiation a non-uniform magnetic field is required, and therefore oscillations in an E-P magnetoplasma have to be considered.

Given that the propagation of an electromagnetic wave from a confined electrostatic oscillation is intrinsically a two-dimensional phenomenon, the one-dimensional treatment of the e+e- oscillation as considered in Paper I has to be extended. The simplest possible geometry is cylindrical, in which the oscillations are assumed radial, and the electromagnetic behaviour involves an axial magnetic perturbation and an azimuthal electric field.

Using the cylindrical geometry ( $r, \theta, z$) in which the magnetic field $\vec{B}=
\hat{\vec{z}}~B_z$ is only axial, the electric field $\vec{E}$ lies in the $r,\theta$-plane and all plasma variables depend only on r and t. The full non-linear model equations are then (for a cold electron-positron plasma)

 
$\displaystyle {r\dot{n}_+ +(rn_+u_r)'=0}$     (1)
$\displaystyle {r\dot{n}_- +(rn_-v_r)'=0}$     (2)
$\displaystyle {\dot{u}_r+u_ru_r'-u_{\theta}^2/r=(e/m)(E_r+u_{\theta}B_z)}$     (3)
$\displaystyle {\dot{u}_{\theta}+u_ru_r'+u_ru_{\theta}/r=(e/m)(E_{\theta}-u_rB_z)}$     (4)
$\displaystyle {\dot{v}_r+v_rv_r'-v_{\theta}^2/r=-(e/m)(E_r+v_{\theta}B_z)}$     (5)
$\displaystyle {\dot{v}_{\theta}+v_rv_r'+v_rv_{\theta}/r=-(e/m)(E_{\theta}-v_rB_z)}$     (6)

and
    
$\displaystyle {(rE_r)'=(e/\epsilon_0)r(n_+-n_-)}$     (7)
$\displaystyle {(rE_{\theta})'=-r\dot{B_z}}$     (8)
$\displaystyle {B_z'=-\dot{E}_{\theta}/c^2-\mu_0e(n_+u_{\theta}-n_-v_{\theta})}$     (9)
$\displaystyle {0=-\dot{E}_r/c^2-\mu_0e(n_+u_r-n_-v_r)}$     (10)

where n+, n- are the positron and electron number densities; u, v are the positron and electron velocities; . denotes $\partial/\partial t$; ' denotes $\partial /\partial r$; (7) is Poisson's equation; (8) is the single z component of the induction equation; and (9), (10) are the $\theta$and r components of the $\nabla \times \vec{B}$ equation.

This non-linear system of equations does not have a closed form analytical solution, and the electrostatic and electromagnetic modes are mixed. In this preliminary investigation only the normal modes arising from the linearised equations are considered in order to identify a simple mechanism for cross-coupling between electrostatic and electromagnetic behaviour. The linearised equations can be cast in non-dimensional form as

     
$\displaystyle \dot{f}=-a'$     (11)
$\displaystyle \dot{g}=-c'$     (12)
$\displaystyle \dot{a}=\rho+p_1b$     (13)
$\displaystyle \dot{b}=\theta-p_1a$     (14)
$\displaystyle \dot{c}=-\rho-p_1d$     (15)
$\displaystyle \dot{d}=-\theta+p_1c$     (16)
$\displaystyle \dot{\rho}=-p_2(a-c)$     (17)
$\displaystyle \dot{\theta}=-\xi p_3\beta'-p_2(b-d)$     (18)
$\displaystyle \dot{\beta}=-\theta'/\xi$     (19)

where
$\displaystyle r=R\xi$   (20)
$\displaystyle t=T\tau$   (21)


$\displaystyle f=\xi n_+/n_0$     (22)
$\displaystyle g=\xi n_-/n_0$     (23)
$\displaystyle a=\xi T u_r/R$     (24)
$\displaystyle b=\xi T u_{\theta}/R$     (25)
$\displaystyle c=\xi T v_r/R$     (26)
$\displaystyle d=\xi T
v_{\theta}/R$     (27)
$\displaystyle \rho=\xi eT^2E_r/(mR)$     (28)
$\displaystyle \theta = \xi
eT^2E_{\theta}/(mR)$     (29)
$\displaystyle \beta = eTB_z/m$     (30)

and
$\displaystyle p_1=\omega_{\rm c}T$     (31)
$\displaystyle p_2=\omega_{\rm p}^2T^2$     (32)
p3=c2T2/R2     (33)
$\displaystyle \omega_{\rm p}^2=n_0e^2/(\epsilon_0m)$     (34)
$\displaystyle \omega_{\rm c} =eB_0/m$     (35)

with n0 and B0 the equilibrium electron (or positron) number density and equilibrium magnetic field respectively. In the normal mode analysis, B0 and n0 are constant and non-zero, and only in this case there are closed form solutions. Note that . now stands for $\partial /\partial \tau$, and ' denotes $\partial /\partial \xi$.

  
2.1 Electrostatic oscillation

The electrostatic oscillation is characterised by Bz=B0, a constant, requiring $E_{\theta}=0$, the latter condition removing the possibility of an undesirable non-zero time-invariant azimuthal electric field. Hence $\beta=0$, $\theta=0$, implying b=d and a=-c. On substitution into Eqs. (11)-(19) the equation for $\rho$ is obtained:

\begin{displaymath}\dot{\rho}=-2p_2a \end{displaymath} (36)

so that

\begin{displaymath}\ddot{\rho}+2p_2\rho=-2p_1p_2b \end{displaymath} (37)

and thus

\begin{displaymath}\partial_{\tau}(\ddot{\rho}+\Omega^2\rho)=0.
\end{displaymath} (38)

The characteristic oscillation frequency, $\Omega=(2p_2+p_1^2)^{1/2}$, when translated from non-dimensional parameters, yields

\begin{displaymath}\omega^2=2\omega_{\rm p}^2+\omega_{\rm c}^2
\end{displaymath} (39)

where $\omega$ is the oscillation frequency. Note that the magnetic field enters via $\omega_{\rm c}$.

The physical mechanism behind the oscillation is as follows. A positive radial electric field (for example) causes electrons to move radially inward, and positrons radially outward, from their rest positions. Hence a=-c. Each particle however is deflected by the axial magnetic field, executing a partial Larmor orbit with the same azimuthal velocity; hence b=d. The resulting charge separation produces a restoring force causing the particles to return to their initial positions, but overshoot. However, the characteristic time taken to retrace their trajectories depends on the strength of the accelerating electric field (via the number density) and also on the magnetic curvature. The electrostatic oscillation here then is a radial expansion plus torsional twist, similar to the dynamics of a watch-spring.

  
2.2 Electromagnetic wave

When Bz=Bz(r,t), an electromagnetic wave is present with an azimuthal electric field component $\theta=\theta (r,t)$. This gives a radially radiating mode. Hence

                   $\displaystyle \ddot{\theta}$ = $\displaystyle -\xi p_3\dot{\beta}'-p_2(\dot{b}-\dot{d})$ (40)
  = $\displaystyle -\xi p_3(-\theta '/\xi)'-p_2[2\theta-p_1(a+c)].$ (41)

Defining
$\displaystyle \mathcal{L}=\ddot{\theta}-\xi p_3(\theta ' /\xi)'+2p_2\theta$     (42)

then yields
$\displaystyle \dot{\mathcal{L}}=p_1^2p_2(b-d)$     (43)
$\displaystyle \ddot{\mathcal{L}}=2p_1^2p_2\theta-p_1^2\mathcal{L}.$     (44)

Assuming $\theta=y(\xi)\exp(-i\bar{\omega} \tau)$ in order to identify the normal modes, irrespective of initial conditions or transients, yields

 \begin{displaymath}(y'/\xi)'+\kappa^2y/\xi =0
\end{displaymath} (45)

with

 \begin{displaymath}\kappa^2=\frac{\bar{\omega}^2}{p_3}\frac{\bar{\omega}^2-p_1^2-2p_2}{\bar{\omega}^2-p_1^2}
\end{displaymath} (46)

showing that, as a function of $\xi$, the solution for y is

 \begin{displaymath}y \sim\left\{\begin{array}{ll}
\xi H^{(1)}_1(\kappa \xi) & \...
...xi) & \textrm{along\ $\xi\rightarrow 0$ }.
\end{array}\right.
\end{displaymath} (47)

taking into account the boundary conditions at $\xi=0$, to avoid singularities of $E_\theta$, and the boundary conditions at $\xi=\infty$, since we are assuming that the oscillation is radiating into an unbounded constant density plasma. The wave nature of the phenomenon is kept since J1(x)=[H(1)1 (x)+ H(2)1(x)]/2 and therefore for $\xi\rightarrow 0$ the solution is a superposition of two waves: the first described by H(2)1, generated at the source region and moving towards $\xi=0$; the second given by H(1)1 generated at the origin after reflection of the first and propagating towards the source region. Note that only latter form persists for $\xi$ between the source region and infinity, since only retarded solutions are acceptable here. Hence for large values of $\xi$, using the asymptotic expansion of the Hankel function H(1)1, the non-dimensional azimuthal electric field takes the form

\begin{displaymath}\theta\propto\sqrt{\xi}\exp\left[i\left(\kappa\xi-\bar{\omega}\tau\right)\right]
\end{displaymath} (48)

and therefore

\begin{displaymath}E_\theta\propto\theta/\xi\propto\xi^{-1/2}\exp\left[i\left(\kappa\xi-\bar{\omega}t\right)\right].
\end{displaymath} (49)

Since $\kappa$ is to be real, in dimensional terms, the wave frequency $\omega$ is restricted in that

 \begin{displaymath}\kappa^2>0\Rightarrow \omega^2>2\omega_{\rm p}^2+\omega_{\rm c}^2.
\end{displaymath} (50)

This result is not valid if a=-c, since then b=d and $\dot{\mathcal{L}}=0$, eventually yielding a contradiction. It is also worth noting that (45) did not depend on restricting $\rho$; in fact, setting $\rho=0$ yields the same result more readily.

Summarizing, assuming the medium to be an unbounded, constant density plasma with a uniform equilibrium magnetic field, the electromagnetic mode in its simplest state propagates radially and has an azimuthal electric field, and an axial magnetic field. The "wavelength'' of the mode is determined by the behaviour of the Hankel function H(1)1, and its dispersion relation is given by (46).

2.3 The mode coupling mechanism

There are two normal modes: a pure electrostatic oscillation, in which $\theta=0$, a=-c, b=d, $\beta=0$ and $\omega^2=2\omega_{\rm p}^2+\omega_{\rm c}^2$; and a pure electromagnetic mode, for which $\rho=0$, $a \neq -c$, $b \neq d$, $\beta \neq 0$and $\omega^2>2\omega_{\rm p}^2+\omega_{\rm c}^2$. These modes are totally independent, because they have mutually exclusive conditions.

The key to unlocking the coupling mechanism and resolving the conflicting conditions between these modes lies in the velocity fields. If the equilibrium magnetic field is spatially non-uniform, the $\vec{B}\times \nabla B$ drift breaks the azimuthal symmetry in the velocity fields of the positrons and electrons undergoing an electrostatic oscillation, relaxing also the strict conditions of the pure electrostatic mode. The resulting azimuthal current density triggers magnetic fluctuations, which propagate away from the oscillation site provided the magnetic gradient makes the plasma underdense to any such e-m disturbance, at least in this simple linear limit.

This can be seen combining Eqs. (18) and (19), and assuming as before $\theta=\eta(\xi)\exp(-i\bar{\omega} \tau)$, but this time leaving the azimuthal current as a general driving term:

 \begin{displaymath}\left(\frac{\eta'}{\xi}\right)'+\bar{\omega}^2\frac{\eta}{\xi}=\mathcal{P(\xi)}\equiv
-\frac{p_2}{p_3\xi}(b-d)
\end{displaymath} (51)

with the most general solution (Zwillinger 1989)
 
$\displaystyle \eta =\xi H^{(1)}_1( \bar{\omega}~\xi) \int^{\xi}_{\xi_1} {\rm d}...
...playstyle\frac{\zeta
H^{(2)}_1(\bar{\omega}~\zeta)\mathcal{P}(\zeta)}{W(\zeta)}$     (52)

where $W(\xi)=-2~{\rm i}/(\pi \xi)$ is the Wronskian of the Bessel function $J_1(\bar{\omega}\xi)$ and Hankel function $H^{(1)}_1(\bar{\omega}\xi)$. This solution generalises Eq. (47) because taking into account that $\mathcal{P(\xi)}=0$, both for $\xi<\xi_1$ and $\xi>\xi_2$
$\displaystyle {\int^{\xi}_{\xi_1} {\rm d}\zeta
\displaystyle \frac{\zeta H^{(1)...
...omega}~\zeta)
\mathcal{P}(\zeta)}{W(\zeta)}=0\quad \textrm{for\
$\xi<\xi_1$ },}$      
$\displaystyle {\int^{\xi}_{\xi_2} {\rm d}\zeta
\displaystyle \frac{\zeta J_1( \...
...omega}~\zeta)
\mathcal{P}(\zeta)}{W(\zeta)}=0\quad \textrm{for\
$\xi>\xi_2$ },}$     (53)

and therefore

\begin{displaymath}\eta=
\left\{\begin{array}{ll}
\xi H^{(1)}_1( \bar{\omega}~...
...)}{W(\zeta)} &
\textrm{for\ $\xi<\xi_1$ },
\end{array}\right.
\end{displaymath} (54)

emphasising the dependence on the driving azimuthal current densities.

A note of caution is needed here, since the drift is actually a particle phenomenon, and needs careful handling to translate to the fluid context. However, the cold plasma model is unusual in that the linearised equations for a species fluid are identical in form to those for an individual particle, except that the electric and magnetic fields are functions of the total particle dynamics via the current and charge density equations.

Moreover, the azimuthal drift current produced by the magnetic inhomogeneity must remain bounded, in order to be consistent with the linearised modelling elsewhere.

3 Numerical technique

Analytical solutions for steady small amplitude wave motion are presented in Sects. 2.1 and 2.2 for a homogeneous equilibrium plasma. (Note that the solutions in Sect. 2.2 do not describe transients.)

However for the case of an inhomogeneous magnetic field the set of Eqs. (11)-(19) although still linear does not have a closed-form analytical solution, and so we must solve the problem numerically.

The form of Eqs. (11)-(19) is suitable for relatively straightforward numerical integration, with perturbed velocity fields as the initial conditions (see later discussion). The equations for $\dot{f},~\dot{g},~\dot{\theta},$ and $\dot{\beta}$ are amenable to solution via a Lax-Wendroff type method, since they are fundamentally hyperbolic in nature. The remaining equations are ordinary differential equations, and therefore suitable for a Runge-Kutta approach.

A vital aspect of any such numerical modelling lies in ensuring that the $\nabla B$ drift is accounted for, in order to provide the requisite azimuthal current density at the grid values appropriate for the Lax-Wendroff computation. This differential azimuthal velocity is essentially a second-order effect, and if the system is stimulated in the electrostatic mode, it remains the only such effect (since to first order, b-d =0). Given the underlying symmetry, the drift is actually evaluated very close to the grid points by the Runge-Kutta algorithm. However it is clear that this drift must remain bounded in order to maintain acceptable limits in accuracy between the two numerical techniques, meaning that the magnetic gradient cannot be too large.

At each time step it is necessary therefore to combine the Runge-Kutta and the Lax-Wendroff finite difference approaches. However it is possible to discuss the methods independently.

The next sections deal with the numerical approach.

3.1 Runge-Kutta

The equations to be solved here are:

   
$\displaystyle {\dot{r_{\rm p}}=a}$     (55)
$\displaystyle {\dot{a}=\rho+p_1(r_{\rm p})b}$     (56)
$\displaystyle {\dot{q}=2\theta-p_1(r_{\rm p})a-p_1(r_{\rm e})c}$     (57)
$\displaystyle {\dot{r_{\rm e}}=c}$     (58)
$\displaystyle {\dot{c}=-\rho-p_1(r_{\rm e})d}$     (59)
$\displaystyle {\dot{d}=-\theta+p_1(r_{\rm e})c}$     (60)
$\displaystyle \dot{\rho}=-p_2(a-c)$     (61)

where $r_{\rm p}$ and $r_{\rm e}$ are the instantaneous radial positions of the electrons and positrons, respectively, and q=b-d, is a direct and more accurate calculation of the azimuthal current density responsible for the generation of the electromagnetic wave. This calculation takes into account the displacement of the electrons and positrons around the mesh point. In fact bearing in mind the caveat in the previous section, the treatment here is a hybrid one, in that the positron fluid "particle'' can be integrated along its trajectory into a region of different magnetic field strength, whilst retaining a global fluid interpretation in terms of the current density. The ambiguity in such an approach is intrinsic to the construction of the cold plasma model itself, meaning that the values of a, b, c, dand $\rho$ are all provided at the perturbed positions.

Finally, since the magnetic field gradient lies in the radial direction, we anticipate that the drift velocity will lie in the azimuthal direction. Strictly, the updating of the local magnetic field strength corresponding to the particle's perturbed position is a second-order correction, and therefore is really a non-linear correction. However, the magnetic field that appears on the right-hand side as the parameter p1 is the equilibrium one, and its gradient is prescribed. Note that the $\nabla B$ must be significant over a Larmor orbit, otherwise there can be no drift, and so it could be argued that the size of this gradient makes this second-order contribution technically of a similar magnitude to the first order ones, and certainly bigger than any $\vec{v}\cdot \nabla \vec{v}$ advective non-linearity. Hence some preliminary justification for updating the field strength along the particle trajectory is possible, particularly if these corrections are seen as a quasi-linear coupling between the linear, homogeneous modes.

3.2 Lax-Wendroff

The standard 1-spatial dimension explicit Lax-Wendroff finite difference solver is applied to the remaining equations, cast as a vector system:

\begin{displaymath}\dot{\vec{y}}+ \vec{A}\vec{y}'=\vec{b}
\end{displaymath} (62)

where
$\displaystyle \vec{y}$ = $\displaystyle \left(\begin{array}{c}
\theta \\
\beta
\end{array}\right)$ (63)
$\displaystyle \vec{A}$ = $\displaystyle \left(\begin{array}{cc}
0 & \xi p_3 \\
1/\xi & 0
\end{array}\right)$ (64)
$\displaystyle \vec{b}$ = $\displaystyle \left(\begin{array}{c}
-p_2 q \\
0
\end{array}\right).$ (65)

Since
           $\displaystyle \by_m^{n+1}$ = $\displaystyle \exp(k \partial_t)\by_m^n$ (66)
  $\textstyle \approx$ $\displaystyle (1+k\partial_t+\frac{1}{2}k^2\partial_t^2)\vec{y}$ (67)

we can use
                    $\displaystyle \partial_t \vec{y}$ = $\displaystyle \vec{b}-\mathbf{A}\partial_x\vec{y}$ (68)
$\displaystyle \partial_t^2\vec{y}$ = $\displaystyle \partial_t \vec{b}
-\mathbf{A}\partial_x[\vec{b}-\mathbf{A}\partial_x\vec{y}]$ (69)

to produce the algorithm for updating $\beta$ and $\theta$ provided

\begin{displaymath}\sigma_m = p^2p_3 \le 1
\end{displaymath} (70)

in order to satisfy the CFL criterion for stability, namely $\vert p\lambda_i\vert\le 1$ where the $\lambda_i$ are the eigenvalues of $\vec{A}$.

The time derivatives $\dot{\vec{b}}$, are supplied via the R-K Eqs. (55)-(61), which are solved simultaneously with the L-W equations.

3.3 Initial values and boundary conditions

To simplify matters, the velocity field is perturbed as the initial condition, driving the plasma evolution. As in Paper I the perturbations are confined to a sub-range of the full spatial mesh: in other words, if the spatial range is [1, mm], then the initial perturbation ranges over $[m_i, m_f]\equiv[m_2-m_r,
m_2+m_r]$, with m2=mm/2. Moreover, in order to avoid sharp gradients at the edges of the initial perturbation, the form of the initial disturbance across its whole range 2mr is multiplied by a Gaussian envelope.

The perturbation is calculated as follows:

$\displaystyle {f(m)=\sin\left(\pi\frac{m-m_i}{m_r}\right)\exp\left[-\left(\lambda\frac{m-m_2}{m_r}\right)^2\right]}$     (71)
$\displaystyle {r_{\rm p}(mh)=mh+r_{p_0} f(m)}$     (72)
a(mh)=a0 f(m)     (73)
q(mh)=(b0-d0) f(m)     (74)
$\displaystyle {r_{\rm e}(mh)=mh+r_{e_0}f(m)}$     (75)
c(mh)=c0 f(m)     (76)
d(mh)=d0 f(m).     (77)

The initial conditions allow the excitation of both the electrostatic and electromagnetic modes separately when $\nabla
B=0$. The electrostatic mode is excited through a0=0.1, c0=-0.1 over a fixed range of radial mesh points, $b, d, \rho,\theta, \beta=0$. The electromagnetic mode in these circumstances is excited through $r_{p_0}\neq 0$, $r_{e_0}\neq 0$, $b_0\neq 0$ and $d_0\neq 0$. Of course when $\nabla
B\neq 0$ then rp0=re0=b0=d0=0.

The values of p2 and p3 are chosen to be consistent with non-relativistic velocities as in Paper I. Moreover, all perturbation coefficients are chosen such that the results do not violate the linear nature of the calculation.

When calculating the mode coupling due to the magnetic gradient, the initial conditions here are identical to those for the electrostatic mode, but this time the magnetic field has a constant gradient, increasing linearly with increasing spatial mesh point. Hence the magnetic field takes the form

p1=p10+p11 x (78)

with x being the position coordinate.

The only boundary condition at finite distance is at $\xi=0$. Since the region of excitation is far from the origin, two waves must propagate away from the source region in opposite directions. There is therefore a delay before the inward propagating wave reaches this boundary and is reflected, meaning that it is possible to stop the calculation before such a reflection takes place, simplifying the numerics. All the calculations reported here employ this tactic.


 
Table 1: Parameters for velocity perturbation computation.
mm 500 number of spatial points
n 5000 number of time steps
h 0.7 spatial mesh increment
p 0.3 ration spatial/temporal meshes
$\lambda$ 10 Gaussian coefficient
p10 1.0 magnetic field amplitude
p11 3.0 e-3 magnetic gradient
p2 4 mesh plasma frequency, squared
p3 1.2e-2 normalized c2
mr 25 half-width of initial perturbation


4 Numerical results

Using the stated algorithms, it was possible to produce a set of numerical results, separating the normal modes from the coupling of them through the magnetic gradient. In all figures only 1% of the time steps was plotted for reasons of clarity.

4.1 Normal modes

The electrostatic mode, with radial electron and positron motion yielding a radial electric field, is shown in Figs. 1 and 2 for the case of a uniform magnetic field. The simulations show a confined and regular oscillation in the radial components; all other quantities are identically zero. As expected no electromagnetic wave is present.


  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{rhohomr.eps} \end{figure} Figure 1: Computed radial electric field for a pure electrostatic oscillation in a uniform magnetic field.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{velrehomr.eps}\end{figure} Figure 2: Computed radial velocity component for electrons for a pure electrostatic oscillation in a uniform magnetic field.
Open with DEXTER

The electromagnetic mode for a uniform magnetic field might be initiated by the boundary condition b0=0.1, d0=-0.1 applied at time $\tau=0$ over the same restricted range of radial mesh points as the electrostatic case, with $a,c,\rho,\theta,\beta=0$. Since the purpose of this paper is to show the connection between electrostatic and electromagnetic modes, and since the latter do not create the former, this case will not be shown.

4.2 Mode coupling via a magnetic gradient

When a magnetic field gradient is introduced, the results are shown in Figs. 3-10. These first results show the persistence of the electrostatic-type behaviour, but with a transmitted electromagnetic wave escaping from the locality of the original disturbance. This agrees with the discussion in Sect. 2, where it was stated that the magnetic field gradient forces an azimuthal current in the "electrostatic'' mode via the magnetic drift, and an electromagnetic wave is generated. In this case given the initial conditions, and our tactic of stopping the numerical solution before the transient wave reaches the origin ($\xi=0$), there are two waves propagating in opposite directions. Note there is no contradiction with solution (52), where the boundary condition at $\xi=0$ is taken into account; the numerical solution is simply at sufficiently early times.

Note that Eq. (51) suggests that the propagating wave acquires its characteristic frequency from the driving term, namely the drift current. This latter frequency is double that of the uncoupled electrostatic mode, since it is produced by the higher-order quadratic coupling inherent in q, as seen in the last two terms on the right-hand side of Eq. (57). Hence the generated electromagnetic wave automatically satisfies Eq. (50), and so is able to propagate. This is evident in Figs. 6 and 8, when compared with Fig. 4, and is presented more clearly in Fig. 10, where the temporal behaviour at a point in the centre of the oscillation is compared with that in the propagation region; the frequency in the propagation zone is nearly double that in the central oscillation.

The radiating fields obtained are very small simply because the coupling is a second-order effect, and therefore limited in magnitude for consistency. In a fully non-linear treatment this restriction does not apply.


  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{velregradr.eps}\end{figure} Figure 3: Computed radial velocity component for electrons in a magnetic field gradient.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{rhogradr.eps}\end{figure} Figure 4: Computed radial electric field component in a magnetic field gradient.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.6cm,clip]{rhog550r.eps}\end{figure} Figure 5: Snapshot of evolution of radial electric field component in a magnetic field gradient, for time step 500 (solid line) and 5000 (dotted line).
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{thegradr.eps}\end{figure} Figure 6: Computed azimuthal electric field component in a magnetic field gradient.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.6cm,clip]{theg550r.eps}\end{figure} Figure 7: Snapshot of evolution of azimuthal electric field component in a magnetic field gradient, for time step 500 (solid line) and 5000 (dotted line).
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{betgradr.eps}\end{figure} Figure 8: Computed perturbed axial magnetic field component, in a magnetic field gradient.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.6cm,clip]{betg550r.eps}\end{figure} Figure 9: Snapshot of evolution of axial magnetic field component in a magnetic field gradient, for time step 500 (solid line) and 5000 (dotted line).
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.6cm,clip]{rhother.eps}\end{figure} Figure 10: Comparison of the time evolution of the radial electric field at the centre point of the initial disturbance (solid line, left-hand axis) with the azimuthal electric field at a point in the propagation zone (dotted line, right-hand axis). Note the frequency doubling in the latter, due to the higher-order drift current interaction.
Open with DEXTER

5 Concluding comments

This paper shows that provided a cold e-p magnetoplasma has an inhomogeneous magnetic field, which are conditions present in pulsar magnetospheres and in particular the magnetospheric plasma rest frame, any electrostatic perturbation will be able to generate an electromagnetic wave propagating away in all directions by virtue of the drift currents. This was demonstrated in the context of a linear approximation, but we expect that a general full non-linear calculation will confirm these findings. Moreover, this mechanism will undoubtedly play a significant part in quenching the density instability found in Paper I, by radiating away the total energy content of the large amplitude oscillation before the onset of the large velocities that are characteristic of unstable density evolution. However, given that this linearised calculation is necessarily limited in the radiation amplitude that it can generate, only a full non-linear calculation show instability quenching in this way.

In fact, the full, non-linear evolution of the plasma will automatically produce large inhomogeneities in all the fields, so that there is the possibility that even with an initially homogeneous magnetic field, electromagnetic radiation might still be created. This simulation is currently underway, and will be reported here soon.

It is appropriate also to note that in the cold, fluid context the existence of a strong magnetic field gradient will mean that the equilibrium density distribution of particles in practice will not be uniform, since the plasma will tend to avoid regions of strong magnetic pressure; however, the equilibrium density is assumed constant here, though the numerical algorithms can be extended to incorporate the effect.

Finally, it is expected that in a pulsar magnetosphere these plasma rest frame phenomena will assume great significance in the laboratory frame as stated in Paper I. Since we might still expect the quasi-static condition to hold, and therefore this mechanism to rotate around with the angular velocity of the pulsar, the radiated frequency in the rest frame has to be Doppler shifted to the laboratory frame taking into account the Lorentz factor of the moving plasma. The linear calculation shows frequency doubling arising from the drift current; a full non-linear calculation will yield a characteristic spectrum of frequencies reflecting local ambient plasma conditions at the oscillation site. If a broad space region is implicated in such radiating oscillations, then it seems clear that the resulting radiation will be broadband even in the plasma rest frame. If this source is rotating outside the speed of light cylinder (da Costa & Kahn 1985), radiation generation in the rest frame might be the basis for a real superluminal lighthouse effect.

Acknowledgements
The present work was developed in two visits, one (DAD) to Lisbon, May 2001, supported by both Centro de Electrodinâmica, and the Astronomy & Astrophysics Group of Department of Physics and Astronomy, University of Glasgow; the other (AAC) to Glasgow, October 2001 supported by the PPARC Short Term Visitors Grant PPA/V/S/2000/00572.

References

 


Copyright ESO 2002