A&A 476, 17-30 (2007)
DOI: 10.1051/0004-6361:20077988

Nonlinear mode coupling in pair plasmas

C. R. Stark¹ - D. A. Diver¹ - A. A. da Costa² - E. W. Laing¹

1 - Dept. of Physics and Astronomy, University of Glasgow, Glasgow G12 8QQ, Scotland, UK
2 - Secção de Telecomunicações, DEEC/SPR, Instituto Superior Técnico, 1049-001 Lisboa, Portugal

Received 31 May 2007 / Accepted 30 September 2007

Abstract
Pulsar magnetospheres are composed of electron-positron plasmas characterised by broadband electromagnetic emission, the source of which remains elusive. This paper investigates one possible emission mechanism in which electrostatic oscillations are coupled to propagating electromagnetic waves by the magnetic field inhomogeneity, thus creating a source of radiation in the pulsar magnetosphere. The full nonlinear equations in cylindrical geometry for a streaming cold electron-positron plasma are solved numerically, together with Maxwell's equations, using a Finite-Difference Time Domain method. Electrostatic oscillations are induced in a streaming plasma in the presence of a non-uniform magnetic field, and the resulting electromagnetic waves are modelled self-consistently. Also presented is the linear perturbation analysis of these model equations perturbed from a dynamical equilibrium in order to probe the fundamental modes present in the system. These simulations successfully exhibit the coupling mechanism and the nonlinear interaction between electromagnetic waves and independent plasma oscillations, confirming the importance of coherent plasma effects and collective plasma processes in the pulsar magnetosphere.

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

1 Introduction

After forty years of study the pulsar emission mechanism still eludes researchers; the theory of pulsar radiation is still in its infancy and there is lack of understanding about the energetic processes involved. The recent observation of Rotating Radio Transients (McLaughlin et al. 2006), periodically active pulsars (Kramer et al. 2006) and pulsar binary systems (Graham-Smith & McLaughlin 2005) show there is a wealth of emission phenomena.

Pulsar magnetospheres are composed of magnetised electron-positron (pair) plasmas characterised by ultra strong magnetic fields and broadband electromagnetic emission. Conventional modelling uses the fact that the dominant Lorentz force produces a host of relativistic charged particles, each of which radiates strongly and stochastically, producing $\gamma$-ray photons. Such single-particle models have been explored as possible radiation sources (da Costa & Kahn 1997), but have not been able to recover the highest energy radiation observed.

However this strict single particle approach, used in many pulsar models, doesn't exploit the properties of the pair plasma; it has been assumed that the magnitude of the pulsar electromagnetic fields dominate the influence of plasma effects. This is not the case since in the plasma rest frame local collective electric fields become significant when transformed into the rotating pulsar laboratory frame. Such fields are comparable with that of the dipolar low-frequency electromagnetic wave driving the magnetosphere (da Costa et al. 2001, hereafter referred as Paper I). Therefore the study of collective plasma processes is of the greatest importance for the dynamics of the magnetosphere, and the radiation mechanisms of pulsars.

Pair plasmas have gained considerable attention over the recent years due to their existence in a number of astrophysical environments such as pulsar magnetospheres. The various wave modes that exist in a pair plasma when the background is homogeneous has been extensively covered in the cold plasma regime (Stewart 1992,1993). Generation of magnetic fluctuations by field-aligned flows in plasmas (Shukla & Shukla 2007) shows a new range of instabilities relevant to strongly flowing constant-density plasmas; further instability studies in homogeneous pair plasmas have been undertaken by Marklund (Marklund et al. 2006), relevant to bursty systems; under certain simplifying conditions, pair-plasma waves can be described by the KdV equation (Verheest & Cattaert 2004). There are also studies of the plasma electromagnetic wave coupling from actual pair production (Bulanov et al. 2005). Our treatment will deal with a cold, free-streaming (but not field-aligned) anisotropic and inhomogeneous pair plasma, using the full nonlinear equations.

In cold non-relativistic plasma theory non-linear electrostatic oscillations in electron-positron plasmas develop a density instability in which the density of both species grows sharply at the edges of the oscillation site (Paper I). Folding thermodynamics into the system provides a possible mechanism for avoiding the onset of the instability, since pressure effects would oppose the density build up. Coupling the oscillation to an electromagnetic mode via an inhomogeneous background magnetic field would provide another means to avoid the onset of the instability. This would allow energy to be radiated away from the oscillation, quenching the density instability and giving a source of radiation in the pulsar magnetosphere. Ultimately thermodynamics will have to be incorporated into a more complete model but it would be advantageous to explore the coupling before relaxing the cold plasma approximation. Previous work has studied this mechanism in the quasi-linear regime (Diver et al. 2002, hereafter referred as Paper II) in which the background magnetic field was inhomogeneous but the plasma density was uniform. A fully consistent nonlinear treatment requires that both the density and magnetic field be inhomogeneous simultaneously.

This paper addresses the nonlinear generation of electromagnetic waves by large amplitude electrostatic oscillations, laying out the details of the equilibrium. The basic mathematical formulation of the electron-positron cold magnetoplasma and the mode coupling mechanism is described in Sect. 2. In Sect. 3 linear analysis of the equations perturbed from dynamical equilibrium is performed to gain an understanding of the modes present in the plasma. Section 4 describes the nonlinear numerical simulations, presenting and discussing the results.

2 Model equations

In cylindrical polar coordinates $\left(r,\theta,z\right)$ with an axial magnetic field $\vec{B}=\vec{\hat{z}}B_{z}$ and an electric field in the $r,\theta$ plane ${E}_{r,\theta}$ the full non-linear model equations for a cold non-relativistic electron-positron plasma are

     

$$r\dot{n}_{+} +(rn_{+}u_{r})'=0$$ (1)

$$r\dot{n}_{-} +(rn_{-}v_{r})'=0$$ (2)

$$\dot{u}_r+u_ru_r'-u_{\theta}^2/r=(e/m)(E_r+u_{\theta}B_z)$$ (3)

$$\dot{u}_{\theta}+u_ru_r'+u_ru_{\theta}/r=(e/m)(E_{\theta}-u_rB_z)$$ (4)

$$\dot{v}_r+v_rv_r'-v_{\theta}^2/r=-(e/m)(E_r+v_{\theta}B_z)$$ (5)

$$\dot{v}_{\theta}+v_rv_r'+v_rv_{\theta}/r=-(e/m)(E_{\theta}-v_rB_z)$$ (6)

$$(rE_r)'=(e/\epsilon_0)r(n_+-n_-)$$ (7)

$$0=-\dot{E}_r/c^2-\mu_0e(n_+u_r-n_-v_r)$$ (8)

$$(rE_{\theta})'=-r\dot{B_z}$$ (9)

$$B_z'=-\dot{E}_{\theta}/c^2-\mu_0e(n_+u_{\theta}-n_-v_{\theta})$$ (10)

where n₊, n_- are the positron and electron number densities; u, v are the positron and electron velocities; ^. denotes $\partial/\partial t$ and ' denotes $\partial /\partial r$. Equation (7) is Poisson's equation for the electric field in the plasma; (9) is the single z component of the $\nabla \times\vec{E}$ equation describing Faraday's law of induction; and (8), (10) are the r and $\theta$ components of the $\nabla \times \vec{B}$ (Ampère-Maxwell) equation. It is assumed that the electrostatic oscillations are radial and the electromagnetic behaviour consists of an axial magnetic perturbation and an azimuthal electric field. The equations can be recast into a form that highlights the symmetry of the electron-positron plasma, this can be done via the following relations

            

$$\Sigma=\frac{1}{2}\left(n_{+}+n_{-}\right)\xi/n_{0}$$ (11)

$$\Delta=\frac{1}{2}\left(n_{+}-n_{-}\right)\xi/n_{0}$$ (12)

$$\sigma=\frac{1}{2}\left(u_{r}+v_{r}\right)/(\omega_{0}L)$$ (13)

$$\delta=\frac{1}{2}\left(u_{r}-v_{r}\right)/(\omega_{0}L)$$ (14)

$$\chi=\frac{1}{2}\left(u_{\theta}+v_{\theta}\right)/(\omega_{0}L)$$ (15)

$$\zeta=\frac{1}{2}\left(u_{\theta}-v_{\theta}\right)/(\omega_{0}L)$$ (16)

$$\rho=\frac{eE_{r}}{m\omega_{0}^{2}L}$$ (17)

$$\theta=\frac{eE_{\theta}}{m\omega_{0}^{2}L}$$ (18)

$$\beta_{0}+\beta=\frac{eB_{z}}{m\omega_{0}}$$ (19)

$$r=\xi L$$ (20)

$$t=\tau/\omega_{0}$$ (21)

$$\omega_{0}^{2}=n_{0}e^2/(\epsilon_{0}m)$$ (22)

where $L,1/\omega_{0}$ are appropriate characteristic length, time scales and $\beta_{0}$, $\beta$ are the equilibrium and perturbed magnetic field respectively. The system of governing equations then becomes

     

$$\dot{\Sigma}=-\left(\Sigma \sigma + \Delta \delta \right)'$$ (23)

$$\dot{\Delta}=-\left(\Delta \sigma + \Sigma \delta \right)'$$ (24)

$$\dot{\sigma}=-\frac{1}{2}(\sigma^{2}+\delta^{2} )'+(\chi^{2}+\zeta^{2})/\xi+\zeta\left(\beta_{0}+\beta\right)$$ (25)

$$\dot{\delta}=-\left( \sigma \delta \right)'+2\chi\zeta/\xi+\rho+\chi\left(\beta_{0}+\beta\right)$$ (26)

$$\dot{\chi}=-\chi'\sigma-\zeta'\delta-\left(\chi\sigma+\zeta\delta\right)/\xi-\delta \left(\beta_{0}+\beta\right)$$ (27)

$$\dot{\zeta}=-\sigma\zeta'-\delta\chi'-\left(\sigma\zeta+\delta\chi\right)/\xi+\theta-\sigma \left(\beta_{0}+\beta\right)$$ (28)

$$(\xi\rho)'=2\Delta$$ (29)

$$\dot{\rho}=-\frac{2}{\xi}\left(\Delta\sigma+\Sigma\delta\right)$$ (30)

$$\dot{\theta}=-p\left(\beta_{0}+\beta\right)'-\frac{2}{\xi}\left(\Sigma\zeta+\Delta\chi\right)$$ (31)

$$\dot{\beta}=-\theta'-\theta/\xi$$ (32)

where $p=c^{2}/(\omega_{0}^{2}L^{2})$ is the dimensionless speed of light squared in the plasma. Note Poisson's equation, Eq. (29), and the r component of the Ampère-Maxwell equation, Eq. (30), are not both required and have been included for completeness.

2.1 The coupling mechanism

Consider a homogeneous background magnetic field. When the density of the plasma is perturbed, resulting in a charge imbalance, a radial electric field is created that accelerates the electrons and positrons in opposite directions. The magnetic field causes the particle trajectories of both species to participate in partial Larmor orbits with the same azimuthal velocity. The particles overshoot their initial positions, due to their acquired kinetic energy, and produce a new charge imbalance. The plasma oscillates at the hybrid frequency $\omega_{\rm H}$ given by

$$\omega_{\rm H}^{2}= 2\omega_{\rm p}^2+\omega_{\rm c}^2$$ (33)

where $\omega_{\rm p}^{2}$ is the plasma frequency and $\omega_{\rm c}^{2}$is the cyclotron frequency. If the magnetic field is inhomogeneous there is a $\vec{B}\times\nabla B$ drift which causes a net current density in the azimuthal motion of the plasma during the electrostatic oscillation. The resulting current density induces axial magnetic field fluctuations and an accompanying time-varying azimuthal electric field. These together constitute a plasma electromagnetic wave that propagates radially away from the electrostatic oscillation site. Hence an electrostatic oscillation in the presence of a non-uniform magnetic field couples to a propagating electromagnetic wave.

Having an inhomogeneous magnetic field permeating the plasma requires the equilibrium to be non-uniform for the system to be self-consistent.

2.2 Dynamical equilibrium

In equilibrium there is no time evolution, we shall assume that the plasma is overall electrically neutral and that there is no electric field present. Hence $\partial/\partial \tau=0$, and $\Delta=\rho=\theta=\beta=0$ implying $\delta=\chi=0$ from Eqs. (30), (26) respectively. The system of governing equations become:

    

$$-(\Sigma\sigma)'$$                     0 = (34)

$$-\frac{1}{2}(\sigma^{2})'+\zeta^{2}/\xi+\beta_{0}\zeta$$ 0 = (35)

$$-\sigma(\zeta'+\zeta/\xi+\beta_{0})$$ 0 = (36)

$$-p\beta_{0}'-2\Sigma\zeta/\xi.$$ 0 = (37)

The resulting equations show that the plasma is in a dynamical equilibrium described by

    

$$\Sigma_{0} \kappa_{0}/\sigma_{0}$$ = (38)

$$\sigma_{0} \sqrt{\kappa_{1}^{2}-\zeta^{2}_{0}}$$ = (39)

$$\zeta_{0}' -\frac{\zeta_{0}}{\xi}-\beta_{0}$$ = (40)

$$\beta_{0}' -\frac{2\kappa_{0}}{p \xi}\frac{\zeta_{0}}{\sqrt{\kappa_{1}^{2}-\zeta_{0}^{2}}}$$ = (41)

where 0 subscripts denote equilibrium values and $\kappa_{0}, \kappa_{1}$ are naturally occurring non-dimensional constants. The equilibrium equations describe the self-consistent response of the plasma to a prescribed inhomogeneous background magnetic field, $\beta_{0}$. $\kappa_{0}$ via Eq. (38) defines the conservation of the total number density flux; $\kappa_{1}$ via Eq. (39) defines the kinetic energy conservation of the equilibrium flow; and Eqs. (40), (41) describe the magnetic field generation via the $\zeta _0$, $\beta_{0}$ coupling. Note that this set of equations assumes $\sigma_{0}\neq 0$, since the equilibrium velocity field is linked to the form of the magnetic inhomogeneity.

3 Linear analysis

The nonlinear system of Eqs. (23)-(32) does not have a closed form analytical solution and must be solved numerically. Insight into the possible dynamical responses of the magnetised pair plasma can be obtained via linear analysis.

Linearise the full set of governing equations and look at large values of $\xi$. In this regime $\sigma_{0}\gg\zeta_{0}$ and $\beta_{0}, \Sigma_{0} \rightarrow$ const., as $\sigma_{0}\rightarrow \kappa_{1}$ corresponding to a net motion of the plasma in the radial direction only, $\zeta_{0}=\chi_{0}=0$. To be consistent with Eq. (40) requires that in the scale length, R, of interest $\zeta_{0}'R\approx\beta_{0}R\ll\kappa_{1}$. The governing equations then become:

         

$$\dot{\Sigma} -\left(\Sigma_{0}\sigma+\Sigma\sigma_{0}\right)'$$ = (42)

$$\dot{\Delta} -\left(\sigma_{0}\Delta+\delta\Sigma_{0}\right)'$$ = (43)

$$\dot{\sigma} -\sigma_{0}\sigma'+\beta_{0}\zeta$$ = (44)

$$\dot{\delta} -\sigma_{0}\delta'+\rho+\beta_{0}\chi$$ = (45)

$$\dot{\chi} -\sigma_{0}\chi'-\sigma_{0}\chi/\xi-\delta\beta_{0}$$ = (46)

$$\dot{\zeta} -\sigma_{0}\zeta'-\sigma_{0}\zeta/\xi-\sigma_{0}\beta-\beta_{0}\sigma+\theta$$ = (47)

$$\dot{\rho} -\frac{2}{\xi}\left(\Delta\sigma_{0}+\Sigma_{0}\delta\right)$$ = (48)

$$\dot{\theta} -p\beta'-2\Sigma_{0}\zeta/\xi$$ = (49)

$$\dot{\beta} -\theta'-\theta/\xi.$$ = (50)

Upon inspection of the linearised governing equations it is evident that they can be split into two independent, self-consistent sets namely Eqs. (43), (45), (46), (48) forming one set and Eq. (42), (44), (47), (49) and (50) forming the other. This simplification was exploited to obtain the following solutions.

3.1 Electrostatic solution

3.1.1 Non-zero streaming

The electrostatic solution is characterised by $\beta=\beta_{0}$, $\theta=0$ and $\zeta=0$. Substituting these conditions into the governing equations yields,

    

$$\dot{\Delta} -(\sigma_{0}\Delta+\delta\Sigma_{0})'$$ = (51)

$$\dot{\delta} -\sigma_{0}\delta'+\rho+\beta_{0}\chi$$ = (52)

$$\dot{\chi} -\sigma_{0}\chi'-\sigma_{0}\chi/\xi-\beta_{0}\delta$$ = (53)

$$\dot{\rho} -2(\Delta\sigma_{0}+\Sigma_{0}\delta)/\xi$$ = (54)

with,

 

$$\dot{\Sigma} -\sigma_{0}\Sigma'$$ (55)

$$\sigma \sigma_{0}$$ (56)

Eq. (55) has solution $\Sigma(\xi,\tau)=f(\tau-\xi/\sigma_{0})$ where f is an arbitrary function. Combining Eqs. (52) and (53) to eliminate $\chi$ as well as Eqs. (54) and (51) to eliminate $\Delta$ produces two differential equations both involving $\delta$ and $\rho$. Substituting one differential expression into the other to eliminate $\rho$ yields the partial differential equation

 

$$\ddot{\delta}+2\sigma_{0}\dot{\delta}'+\sigma_{0}^{2}\delta''+\si... ...\delta'/\xi +(\beta_{0}^{2}+2\Sigma_{0}/\xi)\delta -2\sigma_{0}C_{1}(\xi)/\xi=0$$ (57)

where C₁ is an arbitrary function of $\xi$. Using Eq. (29) instead of (54) in the derivation implies C₁=0, $\forall$ $\xi$. Setting $\delta=y(\xi)\exp{(-{\rm i}\omega\tau)}$ yields

$$\sigma_{0}^{2}y''\!+\!(\sigma_{0}^{2}/\xi\!-\!2{\rm i}\omega\... ...}/\xi\!-\!\omega^{2} \!-\!{\rm i}\omega\sigma_{0}/\xi)y\!=\!0.$$ (58)

For the particular case in which the plasma density varies inversely with $\xi$ so that $\Sigma _0$ = const., we can write the full solution to Eq. (57) as:

 

$$\delta\left(\xi,\tau\right)=\xi^{-1/2}{\rm e}^{{\rm i}\omega(\xi/... ...}}{\sigma_{0}\beta_{0}},0,\frac{2{\rm i}\beta_{0}\xi}{\sigma_{0}}\right)\right]$$ (59)

where M,W are Whittaker functions of the first and second kind respectively (Abramowitz & Stegun 1972) and C₂, C₃ are constants. Note that once again we assume $\sigma_{0}\neq 0$ for Eq. (59) to be valid.

M and W are related to the confluent hypergeometric functions $\mathcal{M}$ and $\mathcal{U}$ as follows

$$M(\kappa,\mu,z) {\rm e}^{-z/2}z^{1/2+\mu}\mathcal{M}(1/2+\mu-\kappa,1+2\mu,z)$$ = (60)

$$W(\kappa,\mu,z) {\rm e}^{-z/2}z^{1/2+\mu}\mathcal{U}(1/2+\mu-\kappa,1+2\mu,z).$$ = (61)

Rewriting the solution in terms of the confluent hypergeometric functions and using the asymptotic expansion of $\mathcal{U}$ and $\mathcal{M}$, the solution takes the form

$$\delta(\xi,\tau) C_{4}\xi^{-1/2}\exp{[{\rm i}\omega(\xi/\sigma_{0}-\tau)-i\beta_{0}\xi/\sigma_{0}-\varphi_{1}]}$$ =

$$+C_{5}\xi^{-1/2}\exp{[{\rm i}\omega(\xi/\sigma_{0}-\tau)+i\beta_{0}\xi/\sigma_{0}+\varphi_{1}]}$$ (62)

where C₄, C₅ are constants and $\varphi_{1}$ is a phase term such that $\vert\varphi_{1}\vert\lesssim \Sigma_{0}\ln{\xi}/(\sigma_{0}\beta_{0})$. In dynamical equilibrium the electrostatic oscillation is being convected at the plasma flow speed $\sigma _0$. The ${\pm} {\rm i}\beta_{0}\xi/\sigma_{0}$ term describes the upstream and downstream movement of the plasma due to particle gyrations around the magnetic field, this effect decreases with increasing $\sigma _0$. The special case where $\Sigma_{0}/\xi$ = const. is dealt with in Appendix B.

3.2 Electromagnetic solution

The equations governing the decoupled electromagnetic response of the plasma are

     

$$\dot{\Sigma} -\left(\Sigma_{0}\sigma+\Sigma\sigma_{0}\right)'$$ = (63)

$$\dot{\sigma} -\sigma_{0}\sigma'+\beta_{0}\zeta$$ = (64)

$$\dot{\zeta} -\sigma_{0}\zeta'-\sigma_{0}\zeta/\xi-\sigma_{0}\beta-\beta_{0}\sigma+\theta$$ = (65)

$$\dot{\theta} -p\beta'-2\Sigma_{0}\zeta/\xi$$ = (66)

$$\dot{\beta} -\theta'-\theta/\xi.$$ = (67)

Differentiating Eq. (66) with respect to $\tau$ and eliminating $\dot{\beta}'$ using (67), one can define the operator $\mathcal{L}_{0}$ as,

$$\mathcal{L}_{0}(\theta)=\ddot{\theta}-p\left[ \frac{(\xi \theta)'}{\xi} \right]'=-2\Sigma_{0}\dot{\zeta}/\xi.$$ (68)

Substituting for $\dot{\zeta}$ from Eq. (65) and collecting all $\theta$ terms together obtain,

$$\mathcal{L}(\theta)=\mathcal{L}_{0}+2\Sigma_{0}\theta/\xi=\fr... ...left[\beta_{0}\sigma+\sigma_{0}(\beta+\zeta'+\zeta/\xi)\right]$$ (69)

where $\mathcal{L}$ is another operator. Defining $\mathcal{K}$ as the operator $\mathcal{K}=\partial_{\tau}+\sigma_{0}\partial_{\xi}$and applying to Eq. (69) yields,

$$\mathcal{K}(\xi\mathcal{L})=2\Sigma_{0}\beta_{0}\mathcal{K}\sigma +2\Sigma_{0}\sigma_{0}\mathcal{K}(\beta+\zeta'+\zeta/\xi)$$ (70)

where $\mathcal{K}\sigma=\dot{\sigma}+\sigma_{0}\sigma'=\beta_{0}\zeta$. The second term in the expression can be expressed as,

$$\mathcal{K}(\beta+\zeta'+\zeta/\xi)=\theta'-\beta_{0}\sigma'+\dot{\beta}+\dot{\zeta}/\xi$$ (71)

where we have substituted from Eq. (65). Returning to Eq. (70) and collecting terms that can be expressed directly in terms of $\theta$ on the left hand side gives,

$$\mathcal{K}(\xi\mathcal{L})-\mathcal{H}=2\Sigma_{0}\beta_{0}(\beta_{0}\zeta-\sigma_{0}\sigma')$$ (72)

where $\mathcal{H}=2\Sigma_{0}\sigma_{0}(\theta'+\dot{\beta}+\dot{\zeta}/\xi)$. Substituting from Eqs. (67) and (69) yields $\mathcal{H}=-\sigma_{0}\mathcal{L}$ hence,

$$\mathcal{K}(\xi\mathcal{L})+\sigma_{0}\mathcal{L}=2\Sigma_{0}\beta_{0}(\beta_{0}\zeta-\sigma_{0}\sigma').$$ (73)

Applying the operator $\mathcal{K}$ to Eq. (73) and using Eqs. (64) and (69) give,

$$\mathcal{K}[\mathcal{K}(\xi\mathcal{L})]+\sigma_{0}\mathcal{K}(\mathcal{L}) 2\Sigma_{0}\beta_{0}[\beta_{0}(\dot{\zeta}+\sigma_{0}\zeta')-\sigma_{0}\beta_{0}\zeta']$$ = (74)

$$-\beta_{0}^{2}\xi(\mathcal{L}-2\Sigma_{0}\theta/\xi).$$ = (75)

Therefore,

$$\mathcal{K}[\mathcal{K}(\xi\mathcal{L})+\sigma_{0}\mathcal{L}]+\beta_{0}^{2}\xi(\mathcal{L}-2\Sigma_{0}\theta/\xi)=0.$$ (76)

Expanding this expression and writing in terms of the parameter $\sigma _0$ yields,

$$\xi\ddot{\mathcal{L}}\!+\!\beta_{0}^{2}\xi(\mathcal{L}-2\Sigm... ...\!+\!\sigma_{0}^{2}[\xi\mathcal{L}''\!+\!2\mathcal{L}']\!=\!0.$$ (77)

3.2.1 Zero streaming

Using a regular perturbation method (Zwillinger 1989) and expressing the solution in terms of a power series in $\sigma _0$: $\theta(\xi,\tau)=g(\xi,\tau)+\sigma_{0}h(\xi,\tau)+\sigma_{0}^{2}l(\xi,\tau)$; and initially including only the zeroth order term in the expansion the solution $\theta=g(\xi,\tau)$ satisfies,

$$\mathcal{D}(\theta)=\ddot{\mathcal{L}}+\beta_{0}^{2}(\mathcal{L}-2\Sigma_{0}\theta/\xi)=0.$$ (78)

Assuming $\theta(\xi,\tau)=y(\xi)\exp{(-{\rm i}\omega\tau)}$ and substituting from Eq. (69) we obtain,

$$\mathcal{D}(y)=\left[ \frac{(\xi y)'}{\xi} \right]'+\frac{\om... ...Sigma_{0}/\xi-\beta_{0}^{2})}{p(\omega^{2}-\beta_{0}^{2})}y=0.$$ (79)

Asserting that $\Sigma _0$ = const., this has solution of the form

 

$$y(\xi) \frac{C_{6}}{\sqrt{\xi}} M\left(\frac{{\rm i}\omega\Sigma_{0}}{\sqrt{p}(\omega^{2}-\beta_{0}^{2})},1,\frac{2{\rm i}\omega\xi}{\sqrt{p}}\right)$$ =

$$+\frac{C_{7}}{\sqrt{\xi}}W\left(\frac{{\rm i}\omega\Sigma_{0}}{\sqrt{p}(\omega^{2}-\beta_{0}^{2})},1,\frac{2{\rm i}\omega\xi}{\sqrt{p}}\right)$$ (80)

where M, W are Whittaker functions of the first and second kind respectively and C₆, C₇ are constants. Hence for large values of $\xi$ the solution takes the form

$$g(\xi,\tau) y(\xi)\exp{(-{\rm i}\omega\tau)}$$ =

$$\sim \frac{C_{8}}{\sqrt{\xi}}\exp{\left[-{\rm i}\omega\left(\xi+\sqrt{p}\tau\right) /\sqrt{p}+\varphi_{2}\right]}$$

$$+\frac{C_{9}}{\sqrt{\xi}}\exp{\left[{\rm i}\omega\left(\xi-\sqrt{p}\tau\right) /\sqrt{p}-\varphi_{2}\right]}$$ (81)

where C_8,9 are complex coefficients and $\varphi_{2}$ is a phase term such that $\vert\varphi_{2}\vert\lesssim \omega\Sigma_{0}\ln{(\xi)}/(\!\sqrt{p}(\omega^{2}-\beta_{0}^{2}))$. The phase velocity of the wave is $\sqrt{p}$, the speed of light in the plasma. The special case where $\Sigma_{0}/\xi$ = const. is dealt with in Appendix B.

3.2.2 Non-zero streaming

Including first order terms in $\sigma _0$ of the perturbation $\theta=g(\xi,\tau)+\sigma_{0}h(\xi,\tau)$ and setting $h(\xi,\tau)=\tilde{y}(\xi)\exp{(-{\rm i}\omega\tau)}$ yields from Eq. (77),

 

$$\mathcal{D}(h) -2\dot{\mathcal{L}}'-3\dot{\mathcal{L}}/\xi$$ = (82)

$${\rm i}\omega(2\mathcal{L}'+3\mathcal{L}/\xi)$$ = (83)

$$\mathcal{D}(\tilde{y}).$$ = (84)

Recalling that $g(\xi,\tau)=y(\xi)\exp{(-{\rm i}\omega\tau)}$, Fourier analysing Eq. (78) and rewriting in terms of $\mathcal{L}$ yields,

$$\mathcal{L} \frac{\mathcal{D}(y)+2\beta_{0}^{2}\Sigma_{0}y/\xi}{(\beta_{0}^{2}-\omega^{2})}$$ = (85)

$$\frac{2\beta_{0}^{2}\Sigma_{0}y/\xi}{(\beta_{0}^{2}-\omega^{2})}$$ = (86)

where $\mathcal{D}(y)=0$ has been used. Substituting into (83) gives,

 

$$\mathcal{D}(\tilde{y}) \left[ \frac{(\xi \tilde{y})'}{\xi} \right]'+\frac{\omega^{2}(\omega^{2}-2\Sigma_{0}/\xi-\beta_{0}^{2})}{p(\omega^{2} -\beta_{0}^{2})}\tilde{y}$$ = (87)

$$\frac{2{\rm i}\omega\beta_{0}^{2}\Sigma_{0}}{(\beta_{0}^{2}-\omega^{2})}[2y'/\xi+y/\xi^{2}]=R(\xi).$$ = (88)

We can now proceed in solving Eq. (88) by obtaining the complementary function from the homogeneous equation and the particular integral from a particular solution of the inhomogeneous equation. Solving $\mathcal{D}(\tilde{y})=0$ yields the solution (80) for $\tilde{y}$. For the particular integral we use the variation of parameters method (Zwillinger 1989),

$$\tilde{y}_{PI}(\xi)=-\tilde{y}_{1}\left(\xi\right)\int{\frac{\til... ...\tilde{y}_{1}\left(\xi\right)R(\xi)}{W(\tilde{y}_{1},\tilde{y}_{2})}}{\rm d}\xi$$ (89)

where $\tilde{y}_{1}$, $\tilde{y}_{2}$ are the two linearly independent solutions to the homogeneous equation ( $y=\tilde{y}_{1}+\tilde{y}_{2}$) and $W(\tilde{y}_{1},\tilde{y}_{2})$ is the Wronskian. Asymptotically expanding the particular integral and integrating yields,

$$\tilde{y}_{PI}\sim C_{10}\xi^{-1/2}\exp{\left[-{\rm i}\omega\xi/\sqrt{p}+\varphi_{2}\right]}$$

$$+C_{11}\xi^{-1/2}\exp{\left[{\rm i}\omega\xi/\sqrt{p}-\varphi_{2}\right]}$$ (90)

where C_10,11 are complex coefficients. Therefore,

$$\theta(\xi,\tau) g(\xi,\tau)+\sigma_{0}h(\xi,\tau)$$ = (91)

$$\approx [(1+\sigma_{0})C_{8}\xi^{-1/2}+\sigma_{0}C_{10}\xi^{-1/2}] {\rm e}^{-{\rm i}\omega(\xi+\sqrt{p}\tau)/\sqrt{p}+\varphi_{2}}$$ (92)

$$+[(1+\sigma_{0})C_{9}\xi^{-1/2}+\sigma_{0}C_{11}\xi^{-1/2}]{\rm e}^{{\rm i}\omega(\xi-\sqrt{p}\tau)/\sqrt{p}-\varphi_{2}}.$$ (93)

In this non-relativistic treatment we must have $\sigma_{0}\ll\sqrt{p}$ implying the radial plasma flow does not affect the propagation of the electromagnetic mode significantly. Hence the zeroth order solution dominates and higher order $\sigma _0$ terms are not required.

3.3 Convective solution

When both electromagnetic and electrostatic modes are suppressed the convective solution is obtained. If $\delta=\theta=\zeta=0$, $\beta=\beta_{0}$ and $\sigma=\sigma_{0}$ this yields

     

$$\dot{\Sigma} -\sigma_{0}\Sigma'$$ = (94)

$$\dot{\Delta} -\sigma_{0}\Delta'$$ = (95)

$$\dot{\chi} -\sigma_{0}(\xi\chi)'/\xi$$ = (96)

$$\dot{\rho} -2\sigma_{0}\Delta/\xi$$ = (97)

$$\rho -\beta_{0}\chi.$$ = (98)

Rearranging Eq. (97) for $\Delta$ and substituting into (95) gives a differential expression for $\rho$. The same result can be achieved by rearranging (98) for $\chi$ and substituting into (96). Therefore,

$$\dot{\rho}=-\sigma_{0}\left(\xi\rho\right)'/\xi.$$ (99)

These have the general solution,

$$F(\xi,\tau) f(\tau-\xi/\sigma_{0})$$ = (100)

$$G(\xi,\tau) \frac{1}{\xi}f(\tau-\xi/\sigma_{0})$$ = (101)

where f is an arbitrary function, $F=\Delta,\Sigma$ and $G=\rho,\chi$. Here F and G are being convected at the streaming velocity of the plasma $\sigma _0$. Assuming $F, G \sim \exp{(-{\rm i}\omega\tau)}$ gives the particular solutions ${\propto}{\exp}{[-{\rm i}\omega(\tau-\xi/\sigma_{0})]}$ describing a perturbation in the plasma variables propagating through the plasma at the streaming velocity $\sigma _0$. Note that the convective solution has no analogous phenomena in the zero streaming case.

4 Nonlinear analysis

4.1 Numerical technique

To solve numerically the nonlinear system of Eqs. (23)-(32) a Finite Difference Time Domain (FDTD) algorithm was employed. The system of equations to be solved consists of a set of first-order hyperbolic partial differential equations (PDEs). FDTD recasts the PDEs as ordinary differential equations (ODEs) with respect to time by replacing the spatial derivatives with finite difference approximations (to fourth-order accuracy here). Defining discrete coordinates $\xi=m$d$\Xi$such that $m\in[0,m_{\max}]$ and using notation $\vec{\Psi}(m$d $\Xi,\tau)=\vec{\Psi}_{m}(\tau)$ we can write for a typical equation:

$$\partial_{\tau} \vec{\Psi}(\xi,\tau) \vec{Q}(\xi,\tau,\vec{\Psi}(\xi,\tau),\partial_{\xi}\vec{\Psi}(\xi,\tau))$$ =

$$\vec{Q}(m\textnormal{d}\Xi,\tau,\vec{\Psi}_{m}(\tau),\partial_{\xi}\vec{\Psi}_{m}(\tau))$$ =

$$\partial_{\tau} \vec{\Psi}_{m}(\tau)$$ = (102)

and

$$\partial_{\xi}\vec{\Psi}_{m}(\tau)=\frac{8(\vec{\Psi}_{m+1}-\... ...1}) -(\vec{\Psi}_{m+2}-\vec{\Psi}_{m-2})}{12\textnormal{d}\Xi}$$ (103)

where $\vec{Q}$ is some function; $m_{\max}$ is the number of spatial grid points; d$\Xi$ is the spatial mesh increment and $\partial_{\xi}$ denotes the partial derivative with respect to $\xi$. The resulting ODEs are then solved via the fourth-order Runge-Kutta method for an integration time $T=n_{\max}$d$\tau$ where $n_{\max}$ is the number of time steps and d$\tau$ is the temporal mesh increment. In comparison to other finite difference methods such as Lax-Wendroff, FDTD is relatively simple to implement.

4.1.1 Initial and boundary conditions

Recall the system of equations describing the dynamical equilibrium of the system, namely Eq. (38)-(41). Defining

$$\zeta_{0} \kappa_{1}\hat\zeta_{0}$$ (104)

$$\beta_{0} \kappa_{1}\hat\beta_{0}$$ (105)

and substituting into (38)-(41) yields

    

$$\Sigma_{0} \frac{\kappa_{0}}{\kappa_{1}}\frac{1}{\sqrt{1-\hat\zeta^{2}_{0}}}$$ = (106)

$$\sigma_{0} \kappa_{1}\sqrt{1-\hat\zeta^{2}_{0}}$$ = (107)

$$\hat\zeta_{0}' -\frac{\hat\zeta_{0}}{\xi}-\hat\beta_{0}$$ = (108)

$$\hat\beta_{0}' -\frac{\kappa_{2}}{ \xi}\frac{\hat\zeta_{0}}{\sqrt{1-\hat\zeta_{0}^{2}}}$$ = (109)

where $\kappa_{2}=2\kappa_{0}/p\kappa_{1}$ controls the strength of the magnetic inhomogeneity via Eq. (109). The system of equations comprise two coupled ordinary differential equations, Eqs. (108), (109), the solution of which gives $\sigma _0$, $\Sigma _0$ via Eqs. (106) and (107). The coupled set were numerically solved using a fourth order Runge-Kutta routine by integrating the equations from a defined starting point $\xi_{i}$ (where $\hat\zeta_{0}$, $\hat\beta_{0}$ are specified) towards $\xi=0$ for a prescribed spatial length, $\Xi\in\xi$ such that $\Xi=m_{\max}$d$\Xi$.

Figure 1: Dynamical equilibrium solutions for $\kappa _{0}=4\times 10^{11}$, $\kappa _{1}=10^{-2}$ and $\kappa _{2}=4\times 10^{9}$. Plots show behaviour as function of position for the following quantities: top left, radial streaming speed $\sigma _0$; top right, total number density $\Sigma _0$; bottom left, differential azimuthal flow $\zeta _0$; bottom right, magnetic field. Solutions are calculated using a bespoke Runge-Kutta routine with the parameters in Table 1. The equilibrium calculated here was used to produce the solutions in Figs. 2-7.

Open with DEXTER

The parameters $\kappa_{0}$, $\kappa_{1}$ and $\kappa_{2}$ characterise the solution and infer the value of p consistent with the system. $\kappa_{2}$ is a key parameter (along with $\Xi$, the initial values of $\hat\beta_{0}(\xi_{i})$, $\hat\zeta_{0}(\xi_{i})$and $\xi_{i}$) since it defines the form of the solution while $\kappa_{0}$ and $\kappa_{1}$ merely scale it. Figure 1 shows a typical example of a calculated equilibrium. In principle $\kappa_{i}$ can be arbitrarily chosen to characterise the system, for instance using a large value of $\kappa_{2}$ to obtain a large magnetic field inhomogeneity. Care has to be taken, however, since the model described here is a non-relativistic one. Recall Eq. (13) for $\sigma$ and define U_r=u_r/c and V_r=v_r/c. In the non-relativistic limit U_r, $V_{r}\ll1$, so $U_{r}+V_{r}\ll1$ hence

$$2\sigma\ll\sqrt{p}.$$ (110)

Setting $\sigma=\kappa_{1}\hat{\sigma}$ and noting in the limit where $\xi \rightarrow \infty$ that $\hat{\sigma}\rightarrow 1$yielding

$$2\kappa^{3}_{1}\kappa_{2}\ll\kappa_{0}.$$ (111)

For smaller values of $\xi$, $\hat{\sigma}<1$. Equation (111) is always true in the non-relativistic limit. Ultimately extending the model to incorporate relativistic effects will allow greater magnetic field strengths and inhomogeneities that are otherwise inaccessible.

Once the equilibrium has been defined the plasma is driven by an initial charge density perturbation,

$$\Delta(m\textnormal{d}\Xi,\tau=0) A_{0} \sin\left[\pi\frac{\left(m-m_{0}+\varsigma\right)}{\varsigma}\right]$$ =

$$\times \exp\left[-\left(\lambda\frac{m-m_{0}}{\varsigma}\right)^{2}\right]$$ (112)

where $\varsigma$ is the half-width of the initial perturbation, A₀ is the amplitude of the perturbation, m₀ is the centre of the computational domain and $\lambda$ is a Gaussian coefficient. To avoid sharp gradients at the edges of the perturbation, the initial disturbance was multiplied by a Gaussian envelope. At $\tau=0$ the initial conditions for the other plasma variables are: $\Sigma_{1}=\sigma_{1}=\delta=\chi=\zeta_{1}=\theta=\beta=0$. The initial condition stimulates the electrostatic mode which couples to an electromagnetic mode via the background magnetic field.

	Figure 2: Spatial and temporal evolution of $\Delta$, showing the plasma density evolution as a function of time $\tau$ and space $\xi$ associated with an electrostatic oscillation. This is the response after a $1\%$ initial charge density perturbation.
Open with DEXTER

	Figure 3: Spatial and temporal evolution of the radial electric field, $\rho$, of the plasma oscillation, consistent with Fig. 2. The nonlinear evolution of the field is visible at the edges of the oscillation site.
Open with DEXTER

	Figure 4: Spatial and temporal evolution of the axial magnetic field, $\beta$, showing clearly the propagation of an electromagnetic signal outwards from the electrostatic oscillation site. Note that there is no magnetic field fluctuation associated with a purely electrostatic phenomenon.
Open with DEXTER

	Figure 5: Spatial and temporal evolution of the azimuthal electric field, $\theta$, of the electromagnetic wave shown in Fig. 4.
Open with DEXTER

	Figure 6: Spatial structure of the axial magnetic field, $\beta$, for time step 800 (dotted line) and 1800 (solid line), corresponding to slices along the $\xi$-axis in Fig. 4, showing that the wave is clearly propagating.
Open with DEXTER

	Figure 7: Spatial structure of the azimuthal electric field, $\theta$, for time step 800 (dotted line) and 1800 (solid line),corresponding to slices along the $\xi$-axis in Fig. 5. Note that the phase of this component, taken with Fig. 6, is consistent with an electromagnetic wave.
Open with DEXTER

Table 1: Table of parameters for numerical simulations Figs. 1-7.

4.2 Numerical results

4.2.1 Mode coupling via inhomogeneous magnetic field

Solving the cold plasma fluid equations yields a continuum description incorporating the coherent motion of the pair plasma. The results shown in Figs. 2-7 were calculated using the equilibrium shown in Fig. 1 and the parameters listed in Table 1, showing the response of the plasma after an initial $1\%$ charge density perturbation: this is in the linear regime. The solutions illustrate the generic form of the coupling in the plasma: when the charge density of the plasma is perturbed, an electrostatic oscillation occurs in which the plasma density fluctuates, Fig. 2, under the influence of the induced radial electric field, Fig. 3. In the presence of an inhomogeneous background magnetic field this oscillation couples to an electromagnetic mode propagating in the radial direction with an axial magnetic field component, Fig. 4, and an azimuthal electric field, Fig. 5. For the $1\%$ charge density perturbation the equilibrium values were: $\Sigma_{0}\sim 4\times10^{13}$; $\sigma_{0}\sim1\times10^{-2}$; $\zeta_{0}\sim 2\times10^{-10}{-}7\times10^{-6}$ and $\beta_{0}\sim5\times10^{-7}{-}2\times10^{-2}$ with $\vert\beta_{0}'\vert\sim1\times10^{-3}{-}6\times10^{1}$. Fully nonlinear simulations are reported later in this section.

The set of parameters $\kappa_{i}$ determine p and hence the quantity $\omega_{0}L$, where the choice of $\omega_{0}$ and Ldefine the physical system to which the results apply. For the solutions in Figs. 1-7 and using Table 1, $p=2\times10^{4}$. Following Paper I and assuming that the required typical electron number densities for $\gamma$-ray radiation in the rest frame is $n_{-}\approx10^{21}$ m^-3, then we have: u_r0, $v_{r0}\sim10^{4}$ m s^-1; $u_{\theta0}$, $v_{\theta0}\sim10^{4}$ m s^-1; $E_{r}\sim10^{9}$ V m^-1; $E_{\theta}\sim10^{-3}$ V m^-1; $B_{z}\sim10^{-3}$ T, where $L\sim10^{-2}$ m. Translated into the laboratory frame in which the pulsar rotates and assuming a typical Lorentz factor $\gamma\sim10^{8}$ gives $B_{z}\sim10^{5}$ T.

The relative emission efficiency of the mode coupling can be characterised by the ratio $\theta/\rho\sim10^{-12}$. This is quite small but it is important to note that the coherent dynamical response of the electrostatic oscillation cannot propagate in the cold plasma context and hence it cannot be observed directly; only the coherent, propagating electromagnetic response can be detected. In contrast to single-particle radiation mechanisms the mode coupling here generates a collective electromagnetic mode, consistent with the background plasma conditions, from the coherent plasma motion.

Figure 8: Dynamical equilibrium solutions for $\kappa _{0}=4\times 10^{7}$, $\kappa _{1}=10^{-2}$ and $\kappa _{2}=4\times 10^{11}$. Plots show behaviour as function of position for the following quantities: top left, radial streaming speed $\sigma _0$; top right, total number density $\Sigma _0$; bottom left, differential azimuthal flow $\zeta _0$; bottom right, magnetic field. Solutions calculated using a bespoke Runge-Kutta routine with the parameters in Table 2. The equilibrium calculated here used to produce solutions in Figs. 9-16.

Open with DEXTER

	Figure 9: Spatial and temporal evolution of $\Delta$ for one period of an electrostatic oscillation after an initial $70\%$ charge density perturbation. Note the nonlinear evolution of the density and the formation of density spikes in contrast to that exhibited in Fig. 2. Note that this plot is of the region of interest between spatial grid points 1150 and 1850.
Open with DEXTER

	Figure 10: Spatial and temporal evolution of the radial electric field, $\rho$, of the plasma oscillation consistent with Fig. 9. Note that this plot is of the region of interest between spatial grid points 1150 and 1850.
Open with DEXTER

	Figure 11: Spatial and temporal evolution of the axial magnetic field, $\beta$, showing clearly the propagation of an electromagnetic signal outwards from the electrostatic oscillation site. Note the hallmark features indicative of the nonlinear response of the plasma. This is more clearly exhibited in Fig. 16.
Open with DEXTER

	Figure 12: Spatial and temporal evolution of the azimuthal electric field, $\theta$, of the electromagnetic wave shown in Fig. 11. Note the hallmark features indicative of the nonlinear response of the plasma. This is more clearly exhibited in Fig. 15.
Open with DEXTER

As the generated mode propagates away from the oscillation site it will encounter regions of differing background density giving a corresponding variation in the plasma frequency, requiring a change in the wave amplitude to conserve energy. This is analogous to waves on a string with an impedance variation, but the plasma context introduces an extra feature: namely the interplay between the driving frequency of the wave and the local plasma frequency.

If the frequency of the wave is very large then the effect of the plasma on the wave propagation becomes less significant. A wave propagating towards smaller values of $\xi$ will encounter a gradual increase in the background plasma density and will eventually be partially reflected and absorbed when the frequency of the wave, $\omega$ falls below that of the local hybrid frequency, $\omega < \omega_{\rm H}(\xi)$ (Paper II). From Figs. 6-7 it is clear that the doppler shift is negligible, consistent with our earlier linear analysis.

Figures 2-7 address small perturbations (1%) to the equilibrium; finite amplitude effects need larger perturbations. The nonlinearity of the mode coupling is specifically exhibited in Figs. 9-16. For these specific simulations the equilibrium, Fig. 8, was calculated using the parameters from Table 2 with an initial $70\%$ charge density perturbation. Over the initial charge density perturbation the equilibrium values were: $\Sigma_{0}\sim 4\times10^{9}$; $\sigma_{0}\sim1\times10^{-2}$; $\zeta_{0}\sim 1\times10^{-6}{-}5\times10^{-5}$ and $\beta_{0}\sim2\times10^{-3}{-}1\times10^{-1}$ with $\vert\beta_{0}'\vert\sim4\times10^{2}{-}2\times10^{4}$. In comparison to the simulations discussed in Figs. 2-7 the ratio $\theta/\rho\sim 10^{-6}$ is much higher. Figures 13 and 14 clearly show the onset of the density instability consistent with Paper I and its consequent effect on the evolution of the electromagnetic wave it generates, Figs. 15-16. For these solutions again setting $n_{-}\approx10^{21}$ m^-3 in the plasma rest frame yields: u_r0, $v_{r0}\sim10^{7}$ m s^-1; $u_{\theta0}$, $v_{\theta0}\sim10^{7}$ m s^-1; $E_{r}\sim10^{9}$ V m^-1; $E_{\theta}\sim10^{3}$ V m^-1; $B_{z}\sim1$ T, where $L\sim5\times10^{-1}$ m.

	Figure 13: Spatial structure of $\Delta$ for the time steps 20(solid line) and 1000 (dotted line). The nonlinear evolution of the plasma oscillation and the onset of the density instability is clearly shown.
Open with DEXTER

	Figure 14: Spatial structure of $\Delta$ for the time steps 1620(solid line) and 2000 (dotted line). The nonlinear evolution of the plasma oscillation and the onset of the density instability is clearly shown. Note the gradual flow of the plasma oscillation downstream by the dynamical background plasma.
Open with DEXTER

	Figure 15: Spatial structure of the azimuthal electric field, $\theta$, for time step 1620 (solid line) and 2000 (dotted line). Note the nonlinear behaviour consistent with Fig. 14.
Open with DEXTER

	Figure 16: Spatial structure of the axial magnetic field, $\beta$, for time step 1620 (solid line) and 2000 (dotted line). Note the nonlinear behaviour consistent with Figs. 14 and 15.
Open with DEXTER

	Figure 17: Variation of the azimuthal electric field $\theta$ of the radiated electromagnetic wave with electrostatic mode amplitude. $\vert\beta _{0}'\vert=3.5\times 10^{-7}$ (arbitrary units) at the centre of the oscillation site.
Open with DEXTER

	Figure 18: Variation of the phase and amplitude of the electric field component, $\theta$, of the EM wave with initial perturbation amplitude in the range 40% to 70% of background density.
Open with DEXTER

Increasing the amplitude (and therefore the energy content) of the initial density perturbation for a given background magnetic field gradient correspondingly increases the amplitude of the generated wave, and makes the wave profile more nonlinear, Figs. 17-19. Increasing the nonlinearity of the oscillation also affects its phase, consistent with the results of Paper I. Note that as the initial charge density perturbation is increased relative to the background plasma density the perturbed fluid velocities become relativistic. Relaxing the non-relativistic restrictions in future code developments will address this.

If the initial perturbation amplitude of the density oscillation is kept constant and the background magnetic field gradient is varied, the amplitude of the resulting EM wave increases nonlinearly as the coupling strength ${\propto}\vec{B}\times\nabla B$ is increased, Fig. 20. Notice that as $\beta _{0}'$ is increased the resulting rate of increase in the EM amplitude begins to slow, hinting perhaps at a maximum fraction of the electrostatic that may be converted in this way.

Table 2: Table of parameters for numerical simulations Figs. 8-16.

	Figure 19: Variation of the phase and amplitude of the electric field component, $\theta$, of the EM wave with initial perturbation amplitude in the range 5% to 30% of background density.
Open with DEXTER

	Figure 20: Variation of the azimuthal electric field $\theta$ of the radiated electromagnetic wave with background magnetic field gradient, $\beta _{0}'$, measured at the centre of the oscillation site, after an initial $10\%$ density perturbation.
Open with DEXTER

4.2.2 Multiple wave interactions

Consider the case where there are a number of individual electrostatic oscillations coupling to propagating electromagnetic modes. There will arise the situation where electromagnetic radiation from one oscillation will interact with another oscillation site and its associated radiation field. In this situation the local magnetic field gradient will be altered momentarily, increasing the wave coupling. To investigate this phenomenon an injected electromagnetic pulse propagated towards an established electrostatic oscillation, Fig. 21. Varying the wavelength of the injected wave and its amplitude will alter the local magnetic gradient hence the resultant coupling. Figure 22 shows the results of 3 different wavelengths of injected wave, and their interaction with the electrostatic oscillation: a reflected component moving upstream; a central oscillating feature; and a transmitted feature moving downstream. Note that Fig. 22 shows not the waves themselves but the difference between the wave solutions as a result of interaction, and the non-interacting case, to assist in the comparison. The resultant is calculated by subtracting the full nonlinear calculation from the linear superposition of the two non-interacting cases (that is, the electrostatic oscillation alone, and the injected electromagnetic wave on its own).

	Figure 21: Plot of axial magnetic field $\beta$ showing injected EM wave (left-hand structure) propagating towards EM wave generated by electrostatic coupling.
Open with DEXTER

Figure 22: Nonlinear response of plasma to EM wave interaction, when the outward directed injected EM wave has passed beyond the edge of the ES oscillation site. Plot shows the difference between the cases where the two EM waves are independently propagating, and when they interact nonlinearly. The solid line is the case where the ratio of the EM wave wavelength to the width of the density perturbation is unity. The dotted and dashed lines correspond to a ratio of 2 and 3 respectively. Note the EM wave amplitude is constant for all ratios. The residual nonlinear response consists of three main features: the reflected feature (left-hand side); the oscillation feature (centre); and a transmitted feature, corresponding to the injected EM wave (right-hand side).

Open with DEXTER

The physical origin of the reflected feature is due to two main phenomena: (1) the spatially and temporally varying refractive index of the electrostatic oscillation as the wave passes through it; and (2) the increased wave coupling, momentarily stimulated as the injected wave traverses the oscillation site. The structure and shape of the reflected feature depends strongly on these two factors. The central feature represents the increased wave coupling at the oscillation site due to the passage of the injected wave. The transmitted feature represents the influence of the interaction on the injected EM wave as it propagates away from the oscillation site.

The results of wavelength variation of the fixed-amplitude injected EM wave encountering a plasma oscillation of varying amplitude are presented in Figs. 23 and 24. These show that the maximum energy in the reflected feature occurs when the EM wavelength is not less than the width of the oscillation for smaller oscillation amplitudes, but this drifts to larger wavelengths as the oscillation becomes more nonlinear in character.

Figure 23: Plot of the magnetic energy of the resultant reflected feature as a function of the wavelength of the injected wave normalised to the width of the ES oscillation site. Note that the resultant is calculated by subtracting the full nonlinear calculation from the linear superposition of the two non-interacting cases (that is, the electrostatic oscillation alone, and the injected electromagnetic wave on its own). Notice that for a given density perturbation, the residual reflected feature is a strong function of wavelength of the injected EM wave peaking close to where the wavelength matches the ES site width. However, this peak response drifts with ES amplitude towards longer wavelengths, reflecting the essential nonlinearity of the coupling.

Open with DEXTER

	Figure 24: Same as for Fig. 23.
Open with DEXTER

Figure 25 shows the results for the central feature, under similar conditions as for Figs. 23 and 24. There does not seem to be a peak in the response here, but instead a monotonic increase.

	Figure 25: Plot of the amplitude of the central oscillating feature as a function of the wavelength of the injected EM wave normalised to the width of the ES oscillation site.
Open with DEXTER

Increasing the amplitude (but maintaining constant wavelength) of the injected wave also increases the local magnetic gradient at the oscillation site, and so we expect to see similar phenomena as in the variable wavelength case. The summary of these results is presented in Figs. 26 and 27, which show broadly similar behaviour to the earlier cases but with one notable difference: the reflected feature shows no maximum amplitude. It is clear that a maximum response is associated only with a resonance in scale-lengths, which is entirely reasonable.

	Figure 26: Plot of the magnetic energy of the residual reflected feature as a function of the injected EM wave amplitude normalised to the to the value of the background magnetic gradient at the centre of the ES oscillation site.
Open with DEXTER

	Figure 27: Plot of the amplitude of the central oscillating feature as a function of the amplitude of the injected EM wave normalised to the to the value of the background magnetic gradient at the centre of the ES oscillation site.
Open with DEXTER

5 Discussion and further developments

This paper has presented a wide-ranging study of the generation of electromagnetic radiation by large amplitude electrostatic oscillations in a streaming, magnetised pair-plasma.

In order to characterise the dynamical responses and small-amplitude modes of such a plasma, we have also presented a linear analysis as a preamble to a fully nonlinear treatment, and also to show continuity with earlier work.

The nonlinear numerical simulations have successfully confirmed the proposed coupling mechanism (Paper II) and displayed the nonlinear nature of the interaction between electromagnetic waves and independent electrostatic oscillations. These results open up a completely new way of describing pulsar radiation mechanisms. They offer the exciting possibility of interpreting sub-pulse structure in terms of local plasma densities and magnetic field gradients in terms of the nonlinear coupling and interactions of the waves and oscillations that are detailed here.

Such phenomena extend the usual pulsar electromagnetic field distribution, that is, the superposition of the underlying dipolar electromagnetic field of the star, plus the self-field of the flowing plasma (da Costa & Kahn 1982; Paper I).

Collective processes in the pulsar rest frame depend very strongly on the local plasma and field conditions, as we have shown, and so present a powerful diagnostic tool if such processes can be identified in the emission structure.

Future considerations include extending the cold plasma treatment to a kinetic one, in order that we can generalise from simple cold electrostatic oscillations to propagating electrostatic waves, such as Bernstein modes (Laing & Diver 2005; Keston et al. 2003). These modes could act as sources of electromagnetic radiation via the coupling mechanism, providing a rich spectrum of electromagnetic waves propagating in the magnetosphere. Such work is currently being considered by the authors.

Acknowledgements

We are gratefully acknowledge the support by PPARC for a studentship for C. R. Stark. The authors are grateful to the anonymous referees for constructive comments that greatly enhanced the presentation of this paper. We are also grateful for travel support from PPARC via Rolling grant PP/C00234/1.

Appendix A: Table of important quantities

Table A.1: Physical meaning of non-dimensionalised plasma quantities.

   
Appendix B: Zero streaming regime

This section describes the connection between the linear perturbation analysis detailed in the main body of this paper and that of Paper II. In Paper II the linear analysis described the modes in the plasma when the background plasma was uniform and stationary; here the linear analysis was conducted in the framework of a uniform, radially streaming plasma.

B.1 Equilibrium

When $\zeta_{0}=\sigma_{0}=0$ the system of governing equations reduces to that of Paper II. Self-consistency requires

       

$$-(\Sigma\delta)'$$                 0 = (B.1)

$$-\frac{1}{2}(\delta^{2})'+\chi^{2}/\xi$$ 0 = (B.2)

$$\chi\beta_{0}$$ 0 = (B.3)

$$-\delta\beta_{0}$$ 0 = (B.4)

$$-\delta(\chi'+\chi/\xi)$$ 0 = (B.5)

$$-2\Sigma\delta\protect/\protect\xi$$ 0 = (B.6)

$$-p\beta_{0}'.$$ 0 = (B.7)

The quantity $\Sigma$ cannot equal zero since this would imply a negative particle number density. As a result this requires that $\delta=0$ in accordance with Eq. (B.6) and is additionally consistent with having a non-zero magnetic field, Eq. (B.4), permeating the plasma. For $\beta \neq 0$ we must also have $\chi=0$ to satisfy Eq. (B.3). These conditions are in agreement with Eqs. (B.1), (B.2) and (B.5). In this regime the background plasma is uniform and therefore the equilibrium magnetic field must be constant as dictated by Eq. (B.7). Note that no constraint is placed on the behaviour of $\Sigma$ in this regime other than $\Sigma \neq 0$.

   
B.2 Electrostatic solution

Recall the differential equation describing the electrostatic response of a radially streaming pair plasma:

$$\ddot{\delta}+2\sigma_{0}\dot{\delta}'+\sigma_{0}^{2}\delta''+\si... ...delta'/\xi +(\beta_{0}^{2}+2\Sigma_{0}/\xi)\delta -2\sigma_{0}C_{1}(\xi)/\xi=0.$$ (B.8)

In the regime where $\sigma_{0}=0$ the differential equation Eq. (57) should reduce to the equivalent expression in Paper II. Setting $\sigma_{0}=0$ the expression becomes,

$$\ddot{\delta}+(\beta_{0}^{2}+2\Sigma_{0}/\xi)\delta=0$$ (B.9)

where the term in the brackets takes the role of the plasma hybrid frequency with the notable exception of the $1/\xi$ dependence. Since $\sigma_{0}=0$ here, there is no information governing the spatial dependence of the equilibrium density. In Paper II the background density n₀ = const. is recovered by setting $\Sigma_{0}/\xi$ = const., using the definition of $\xi$ in Eq. (11). This means that the term in parenthesis in Eq. (B.9) is the square of the constant hybrid frequency $\omega_{\rm H}$, as before.

B.3 Electromagnetic solution

In the regime where $\sigma_{0}=0$ the differential equation governing the electromagnetic response of the plasma, Eq. (79), is:

$$\mathcal{D}(y)=\left[ \frac{(\xi y)'}{\xi} \right]'+\frac{\om... ...Sigma_{0}/\xi-\beta_{0}^{2})}{p(\omega^{2}-\beta_{0}^{2})}y=0.$$ (B.10)

This is equivalent to the expression in Paper II when $\Sigma_{0}/\xi$ = const. Following the same argument described in Sect. B.2 we recover the Paper II result,

$$y(\xi)=C_{12}J_{1}\left(\kappa\xi\right)+C_{13}Y_{1}\left(\kappa \xi\right)$$ (B.11)

where

$$\kappa^{2}=\frac{\omega^{2}(\omega^{2}-2\nu_{0}-\beta_{0}^{2})} {p(\omega^{2}-\beta_{0}^{2})}$$ (B.12)

and $\nu_{0}=\Sigma_{0}/\xi=\rm const$.

References

• Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions (New York, USA: Dover Publications, Inc) (In the text)
• Bulanov, S. S., Fedotov, A. M., & Pegoraro, F. 2005, Phys. Rev. E, 71, 016404 [NASA ADS] [CrossRef] (In the text)
• da Costa, A. A., & Kahn, F. D. 1982, MNRAS, 199, 211 [NASA ADS] (In the text)
• da Costa, A. A., & Kahn, F. D. 1997, MNRAS, 284, 1 [NASA ADS] (In the text)
• da Costa, A. A., Diver, D. A., & Stewart, G. A. 2001, A&A, 366, 129 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Diver, D. A., da Costa, A. A., & Laing, E. W. 2002, A&A, 387, 339 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Graham-Smith, F., & McLaughlin, M. A. 2005, A&G, 46, 1.23 (In the text)
• Keston, D. A., Laing, E. W., & Diver, D. A. 2003, Phys. Rev. E, 67, 036403 [NASA ADS] [CrossRef]
• Kramer, M., Lyne, A. G., O'Brien, J. T., Jordan, C. A., et al. 2006, Science, 312, 549 [NASA ADS] [CrossRef] (In the text)
• Laing, E. W., & Diver, D. A. 2005, Phys. Rev. E, 72, 036409 [NASA ADS] [CrossRef]
• Marklund, M., Eliasson, B., & Shukla, P. K. 2006, Phys. Plasmas, 13, 083102 [CrossRef] (In the text)
• McLaughlin, M. A., lyne, A. G., Lorimer, D. R., Kramer, M., et al. 2006, Nature, 439, 817 [NASA ADS] [CrossRef] (In the text)
• Shukla, N., & Shukla, P. K. 2007, Phys. Lett. A, 362, 221 [NASA ADS] [CrossRef] (In the text)
• Stewart, G. A. 1992, J. Plasma Phys., 47, 295
• Stewart, G. A. 1993, Ph.D. Thesis, University of Glasgow
• Verheest, F., & Cattaert, T. 2004, Phys. Plasmas, 11(6), 3078 (In the text)
• Zwillinger, D. 1989, Handbook of Differential Equations (San Diego, USA: Academic Press) (In the text)