Contents

• A free-electron laser in the pulsar magnetosphere
• 1 Introduction
• 2 The model: A FEL in the pulsar magnetosphere
• 3 Numerical simulations
• 4 Physical parameters for numerical simulations
• 5 Results
• 6 Discussion
• References
• 7 Online Material
• Appendix A: Derivation of Gaussian modes

A&A 422, 817-830 (2004)
DOI: 10.1051/0004-6361:20040295

A free-electron laser in the pulsar magnetosphere

P. K. Fung¹ - J. Kuijpers²

1 - Astronomical Institute, Utrecht University, PO Box 80000, 3508 TA Utrecht, The Netherlands
2 - Department of Astrophysics, University of Nijmegen, PO Box 9010, 6500 GL Nijmegen, The Netherlands

Received 18 February 2004 / Accepted 2 May 2004

Abstract
We have studied systematically the free-electron laser in the context of high brightness pulsar radio emission. In this paper, we have numerically examined the case where a transverse electromagnetic wave is distorting the motion of a relativistic electron beam while travelling over one stellar radius ($\approx$ $10\:\mbox{km}$). For different sets of parameters, coherent emission is generated by bunches of beam electrons in the radio domain, with bandwidths of 3 GHz. Pulse power often reached $10^{13}\:\mbox{W}$, which corresponds with brightness temperature of $10^{30}\:\mbox{K}$. The duration of these pulses is of the order of nanoseconds. In the context of pulsar radio emission, our results indicate that the laser can produce elementary bursts of radiation which build up the observed microstructures of a few tens of microseconds duration. The process is sensitive mostly to the beam particles energy, number density and the background magnetic field, but much less so to the transverse wave parameters. We demonstrate that the operation of a free-electron laser with a transverse electromagnetic wiggler in the pulsar magnetosphere occurs preferably at moderate Lorentz factors $\gamma \geq 100$, high beam density $n \gtrsim 0.1\;n_{{\rm GJ}}(r_\ast)$ where $n_{{\rm GJ}}(r_\ast)$ is the Goldreich-Julian density at a stellar radius $r_\ast$, and finally, at large altitude where the background magnetic field is low $B_0 \leq 10^{-2}\;{\rm T}$.

Key words: masers - radiation mechanisms: non-thermal - plasmas

1 Introduction

Although radio pulsars have been discovered now for more than three decades, the generation of the radio signal which reveals their existence, is still a puzzle. It is believed that a pulsar is a rotating neutron star with a high magnetic field ( $B \sim 10^5 {-} 10^9~\mbox{T}$). The radio pulses span a frequency range from $\sim$ $100\:\mbox{MHz}$ to $\sim$ $30\:{\rm GHz}$ and have the remarkable property that, in a given pulsar, they vary dramatically from pulse to pulse, while the average pulse profile is extremely stable and unique for that pulsar. On one hand, the stability and uniqueness of the averaged pulse shape for every pulsar imply that the radio pulse is generated well inside the light cylinder (defined as the cylinder with radius $R_{{\rm lc}} = c/\Omega_\ast$, where  $\Omega_\ast$ is the rotation frequency of the pulsar and c is the speed of light). On the other hand, the variability of successive pulses suggests that the radio emission process is strongly fluctuating and/or that the acceleration of particles responsible for the emission is highly non-steady. The observations indicate that the radiation is emitted more or less tangential to the open field lines from the magnetic poles and that it is polarised preferentially in the local plane of the magnetic field lines and perpendicular to it. Current belief is that radio emission is related to the development of a cascade of pairs of electrons and positrons on the open field lines anchored at the magnetic poles.

An important property of the pulses is their high brightness temperature  $T_{{\rm b}}$, given by $T_{{\rm b}} \sim 10^{34} \: \mbox{K} \: F_{{\rm Jy}} \: d_{{\rm kpc}}^2\: r_{{\rm m}}^{-2} \: \nu_{{\rm GHz}}^{-2}$, where  $F_{{\rm Jy}}$ is the measured flux in Jansky at the frequency $\nu_{{\rm GHz}}$ for a pulsar at distance $d_{{\rm kpc}}$ and where $r_{\rm m}$is the radius of the emission region (in meters). For characteristic values: $F \sim 1\:\mbox{Jy}$, $d \sim 1\:\mbox{kpc}$ at $\nu \sim 1\:\mbox{GHz}$, this implies that the brightness temperature ranges from $10^{26}\:\mbox{K}$ to $10^{30}\:\mbox{K}$, depending on whether the radius of the emission region is taken to be the whole pulsar ( $r = r_\ast \sim 10^4\:\mbox{m}$) or only the polar cap ( $r = r_{{\rm pc}} = r_\ast (\Omega_\ast r_\ast/c)^{-1/2}\sim 10^2~\mbox{m}$ for $P_\ast = 2\pi/\Omega_\ast = 1\:\mbox{s}$). If the emission process were incoherent this would imply the presence of energetic electrons (positrons) of individual energies ${\cal E} \sim k_{{\rm B}}\:T/m_{{\rm e}} \sim 10^{10} {-} 10^{18}\ \mbox{GeV}$, which are difficult to achieve in view of the maximum available voltage jump in a rotating magnetic star inside the light cylinder. Because of this, the mechanism responsible for the radio emission is assumed to be some form of coherent action. The emission can then be either an antenna or a maser process. If the radiation source is located on the open field lines above the polar cap a maser emission process is most likely powered by a high-energy electron-positron beam. Candidates for the emission process are maser curvature emission (Luo & Melrose 1995), relativistic plasma emission (Melrose & Gedalin 1999), normal and anomalous Doppler instability (Lyutikov 1998), and linear acceleration emission (Melrose 1978).

Despite the efforts in the investigations over the past three decades, the mechanism responsible for the radio emission is still unknown.

The process studied here is a specific type of the free-electron laser (FEL), a laboratory device which produces coherent radiation. In a FEL, a beam of relativistic electrons, with velocity v, passes through a periodic, electromagnetic field (called wiggler) and radiates at the resonance frequency $\nu_{{\rm res}} = \omega_{{\rm res}}/2\pi = \pm(2\gamma^2 (\omega_{{\rm w}} - k_{{\rm w}} v))/2\pi$, where $\gamma={(1-{(v/c)}^2)}^{-1/2}$, and $\omega _{{\rm w}}$ and $k_{{\rm w}}$ are the wave frequency and the wavenumber respectively associated with the wavelength $\lambda _{{\rm w}}$ of the wiggler. Under appropriate conditions the radiation will be enhanced by the particles that follow, due to the bunching of the particles. Previously mentioned linear acceleration emission (Melrose 1978) and coherent inverse Compton scattering (Schopper et al. 2002) are both forms of this mechanism. They differ from the case here in the sense that, in this study, particles will undergo a small transverse displacement, whereas in the other two cases, the distortion is along the path of the particles. Therefore, the solution proposed here is only applicable for the region where the background magnetic field is sufficient small, in comparison to the magnetic field amplitude of the wiggler, i.e. high up in the magnetosphere.

Applying the FEL concept to the pulsar, we will investigate as potential wigglers high-frequency Alfvén waves and other potential periodic structures in the pulsar magnetosphere. Such waves might be generated by a beam instability of the (remnant of the) primary beam in the ambient secondary pair plasma, or by the inhomogeneity of the pair plasma in which faster particles run into dense clumps of pairs (Asseo & Melikidze 1998; Melikidze et al. 2000; Usov 1987). The efficiency of the FEL interaction between these waves and the beam of primary/secondary particles is investigated by doing numerical simulations.

In Sect. 2 the basic concept of the free-electron laser mechanism is presented in detail before it is applied to the pulsar magnetosphere. A description of the simulation code is given is Sect. 3, and a summary of free-electron laser parameters is found in Sect. 4. In Sect. 5, we will present the results of the numerical simulations of this process under pulsar magnetosphere conditions. The conclusions and a discussion are given in Sect. 6.

Throughout the paper, SI units and Cartesian coordinates, in which $\vec{x}$ is the position vector, with local magnetic field parallel to $\hat{z}$ are used.

   
2 The model: A FEL in the pulsar magnetosphere

In this section the basic theory of FEL operation is explained.

Note that the FEL mechanism that we treat here is an antenne process in which bunching occurs in space ("reactive'').

2.1 Free-electron laser theory

2.1.1 1-particle trajectory and resonance frequency

The relativistic equation of motion for a particle in an electromagnetic disturbance (wiggler), with electric field given by  $\vec{E}_{{\rm w}}$ and magnetic field by $\vec{B}_{{\rm w}}$, is given by:

$$\frac{\mbox{d}\gamma m \mbox{\boldmath$\beta$ } c}{\mbox{d}t}... ...rm w}} + \mbox{\boldmath$\beta$ } c\times \vec{B}_{{\rm w}}]$$ (1)

where the electromagnetic field of the wiggler is described by:

 

$$\vec{E}_{{\rm w}}(\vec{x},t)= E_{{\rm w}} \cos(k_{{\rm w}} z - \omega_{{\rm w}} t) \hat{x}$$

$$\vec{B}_{{\rm w}}(\vec{x},t)= B_{{\rm w}} \cos(k_{{\rm w}} z - \omega_{{\rm w}} t) \hat{y}$$

$$\frac{E_{{\rm w}}}{B_{{\rm w}}} = \frac{\omega_{{\rm w}}}{k_{{\rm w}}} = \beta_{{\rm w}} c$$ (2)

and  $\beta _{{\rm w}}$ denotes the wiggler phase velocity normalised to the velocity of light c. For a relativistic particle, with initial Lorentz factor $\gamma_0 = (1 - \beta_{z0})^{-1/2} \gg 1$, the first order approximation of (1), where we neglect the particle's energy change, becomes:

$$\frac{\mbox{d} \beta_x}{\mbox{d} t} = \frac{\sqrt{2}K}{\gamm... ..._{{\rm w}}}\right) \cos~(k_{{\rm w}} z - \omega_{{\rm w}} t).$$ (3)

The solution of this equation is:

$$\beta_x = -\eta\frac{\sqrt{2}K}{\gamma} \sin~ (k_{{\rm w}} z - \omega_{{\rm w}} t)$$ (4)

where $\eta = q/\vert q\vert$ and the strength of the wiggler is re-expressed in a dimensionless parameter K, defined as:

$$K \equiv \left\vert \frac{q B_{{\rm w}} \lambda_{{\rm w}} }{\sqrt{8}\pi mc} \right\vert\cdot$$ (5)

Thus, the effect of the wiggler on the particle orbit is adding a small sinusoidal motion in the direction of the wiggler's electric field, where the periodicity of this motion for an observer in the laboratory frame is given by

$$t_{{\rm eff}} = \frac{\lambda_{{\rm eff}}}{c} = \frac{\lambda_{{\rm w}}}{ (\beta_{z0} - \beta_{{\rm w}}) c}\cdot$$ (6)

Up to the second order, the (angular) frequency of the radiation emitted by a particle moving in this disturbance is derived from $(\omega - \vec{k}\cdot\vec{v}) = \pm (\omega_{{\rm w}} - k_{{\rm w}} v_z)$:

$$\omega_{{\rm res}} = -\frac{2 \gamma^2(\omega_{{\rm w}} - k_{{\rm w}} v_z)}{1 + K^2}$$ (7)

where we have used $\omega - \vec{k}\cdot\vec{v} \approx \omega(1 - \beta_z)$ for radiation beamed forward and $\beta^2_z = \beta^2 - \beta^2_x = 1 - \gamma^{-2} - K^2\gamma^{-2}$.

2.1.2 Collective effect: Bunching

The wiggler field causes a transverse motion of the particles (4). In presence of a radiation field $\vec{E}_{{\rm rad}},\;\vec{B}_{{\rm rad}} \propto \sin(k z - \omega t)$, which is set up by particles radiating at the resonance frequency (7), the particles that follow will move in a beat wave of the wiggler field and the radiation field. This beat is often called the ponderomotive wave, for its spatial energy density is highly non-uniform. The corresponding ponderomotive force is $\vec{F}_{{\rm p}} = q \beta_x c\hat{x} \times {\vec{B}_{{\rm rad}}}$ and therefore, it is proportional to $\sin((k + k_{{\rm w}})z - (\omega + \omega_{{\rm w}})t)$. The beat which is proportional to $\sin((k - k_{{\rm w}}) z - (\omega - \omega_{{\rm w}}) t)$ is superluminal, and, therefore, has no immediate particle-wave resonances.

This force acts on the beam in the z-direction and drives a longitudinal current $\delta J_z = q n v_{{\rm p}}$, where q is the elementary charge of the beam particles, n is the number density of the particles and  $v_{{\rm p}}$ is the velocity which is induced by the motion in the ponderomotive wave. As a result, like $\vec{F}_{{\rm p}}$, $\delta J_z$is proportional to $\sin((k + k_{{\rm w}})z - (\omega + \omega_{{\rm w}})t)$. Related to $\delta J_z$ is a density modulation according to $q \partial (\delta n)/\partial t = - \mbox{\boldmath$\nabla$ }\cdot \delta J_z \hat{z}$. This density modulation is observed as bunching of the particles. Note that the bunches occur at the ponderomotive wavelength

$$\lambda_{{\rm p}} = 2\pi k_{{\rm p}}^{-1}$$ (8)

where $k_{{\rm p}} = k_{{\rm w}} + k$. Due to this bunching, the transverse current is in phase with the ambient field ( $\delta J_\perp = q\delta n \vec{v}_{{\rm w}} \propto \sin(k z - \omega t)$), i.e. particles in the bunch radiate in phase and enhance the ambient radiation field.

2.2 Background magnetic field

The effect of a background magnetic field $\vec{B} = B_0 \hat{z}$ limits the particle's motion perpendicular to the field. The equations of motion for a particle now become:

$$\frac{\mbox{d} \gamma \beta_x}{\mbox{d}t} = \sqrt{2} K \omega_{{\... ...\rm w}}}\right) \sin(k_{{\rm w}} z - \omega_{{\rm w}} t) + \Omega_{B_0} \beta_y$$ (9)

$$\frac{\mbox{d}\gamma \beta_y}{\mbox{d}t} = \Omega_{B_0} \beta_x$$ (10)

$$\frac{\mbox{d}\gamma \beta_z}{\mbox{d}t} = \sqrt{2}K \frac{\beta_x}{\beta_{{\rm w}}} \omega_{{\rm w}}$$ (11)

where K is defined by Eq. (5) and $\Omega_B = \vert q\vert B/m$ is the non-relativistic gyrofrequency for a particle (mass m, charge q) in a magnetic field with strength B. For zero background magnetic field, the equations reduce to (1). In the presence of a uniform background magnetic field, the first order solution for the particle's velocity components (the x- and y-component) are coupled to each other through  $\Omega_{B_0}$. The particle is no longer free to resonate with the wiggler, but is bounded by the magnetic field. This restriction limits the particle's ability to resonate and, therefore, also the FEL laser action.

2.3 Application in the pulsar magnetosphere

Our aim is to apply the mechanism to the pulsar magnetosphere, while making as few assumptions as possible. As is common to most pulsar models, we assume that a large electric field is set up along the open field lines above the polar caps due to the fast rotation of the magnetised neutron star. This electric field pulls out and accelerates electrons, which, in the presence of a strong magnetic field, are in the lowest Landau level and thus move along the magnetic field lines. At some altitude above the pulsar surface - between a few polar cap radii and a few stellar radii - a dense plasma of electron-positron pairs is produced, either from curvature radiation of the primaries in the strong magnetic field or from their inverse Compton radiation. Very likely, the entire process is only stationary in an average sense but highly variable in space and time on small scales with inhomogeneous distributions of pair plasma intermingled with primary beams.

In our model, we investigate radiation from a mono-energetic beam of electrons in the presence of a subluminal transverse electromagnetic wiggler.

In theory, radiation of any frequency can be obtained by tuning  $\omega _{{\rm w}}$, $\vec{k}_{{\rm w}}$ and $\gamma$. For radiation at frequencies $\nu_{{\rm res}} = \omega_{{\rm res}}/2\pi$ in the range of 10⁸ and $10^{10}\;\mbox{Hz}$, and parallel wave vectors  $\vec{k}_{{\rm w}} \parallel \vec{v}$, the required Lorentz factors range from unity up to 10⁷ (Fig. 1).

Figure 1: Values for which $\omega _{{\rm w}}$, $k_{{\rm w}}$ and $\gamma$ give rise to radio emission between $100\ \mbox{MHz}$ (lower solid line) and $10\ \mbox{GHz}$ (upper) according to the resonance condition (7) with K = 10. On the vertical axis is plotted $k_{{\rm w}} v - \omega _{{\rm w}}$, because we assume the particle's velocity to be larger than the wiggler's phase velocity.

Near the stellar surface, particles move one-dimensionally along the magnetic field lines due to extremely fast synchrotron losses in the strong magnetic field. The quantity $p_\perp^2/B_0$ is not invariant anymore when the timescale at which the magnetic field changes, is small compared to one gyration period in the background magnetic field, i.e.:

$$\frac{B_0}{\mbox{d}B_{{\rm w}}/\mbox{d}t} \ll \frac{2 \pi \gamma}{\Omega_{B_0}}$$ (12)

and in our case, using (2):

$$B_0 \ll B_{\rm w}2^{1/4}\pi^{1/2}\sqrt{\frac{(\beta_{z0} - \beta_{{\rm w}})\gamma}{K}}\cdot$$ (13)

So, if the last inequality is satisfied, a particle moving in the background magnetic field acquires momentum transverse to the magnetic field.

Next to consider is the synchrotron loss time, which needs to be larger than the FEL timescale when the latter process causes the pulsar radio emission. The synchrotron loss time  $\tau_{1/2}$, the time within which the particle looses half of its initial energy $E = \gamma m c^2$, is defined as $\tau_{1/2} \equiv E/(2P)$, where $P = \gamma^2 \beta_\perp^2 \Omega_{B_0}^2/(6 \pi \epsilon_0 c)$ is the power of synchrotron radiation emitted by a charged particle moving in a background magnetic field B₀ and $\beta_\perp$ is the component of the normalised velocity perpendicular to B₀. By using the upper limit for B₀ (13), we find for  $\tau_{1/2}$

$$\tau_{1/2} \gg \frac{3 \epsilon_0 mc}{8\sqrt{2} \pi^2 q^2}\frac{\lambda^2_{{\rm w}}}{(\beta_{z0} - \beta_{{\rm w}})K^3}$$ (14)

where we used $\beta_\perp \approx K/\gamma$ from (4).

On timescales much smaller than $\tau_{1/2}$, a change of the particle energy due to synchrotron radiation can be neglected.

In our calculations, we used parameters such that firstly, the resonance is in the radio regime (7), secondly, the inequality (13) is satisfied so that particles can (and do) acquire transverse momentum, and finally, the synchrotron emission can be neglected (14).

3 Numerical simulations

This section describes how the ingredients for the FEL in the pair plasma of a pulsar magnetosphere are represented in our numerical simulations. The Coulomb interactions between the particles in the beam are neglected, due to the large Lorentz factors. This is because, although the beam particles generate an electric field radially outward E_r, due to their relativistic speeds, the generated magnetic field (in the azimuthal direction) reduces the electric field by a factor of $1 - \beta^2$, where $\beta$ is the particles' speed. This results in a reduced radial Lorentz force $F_r = q E_r/\gamma^2$.

Furthermore, the role of the pairplasma is mainly in providing the wiggler. Therefore, we model the presence of the wiggler by Eqs. (2) rather than generating it in a consistent way.

First, we give a brief introduction to the code. Then, the parameters for each ingredient in the case of N-particles simulations are given.

3.1 Code General Particle Tracer (GPT)

We will investigate radiation from a relativistic electron beam travelling through a wiggler, and the formation of bunches is tested by doing numerical simulations with a code called General Particle Tracer (GPT).

This code solves the equation of motion of each (macro)particle in the time domain numerically. For each macroparticle, labelled i, the differential equations

 

$$\frac{\mbox{d}\gamma_i \mbox{\boldmath$\beta$ }_i}{\mbox{d}t} = \frac{q_i}{m_ic}[\vec{E}(\vec{x},t)+ \vec{v}_i(t) \times \vec{B}(\vec{x},t)]$$

$$\frac{\mbox{d}\vec{x}_i}{\mbox{d}t} = \vec{v}_i = \frac{\gamma_i\mbox{\boldmath$\beta$ }_i c}{\sqrt{\gamma_i^2\mbox{\boldmath$\beta$ }_i^2 + 1}}$$ (15)

are solved numerically with a fifth-order embedded Runge-Kutta ODE solver (Press et al. 1992). Here $\gamma = (1 - \mbox{\boldmath$\beta$ }^2)^{-1/2}$ is the Lorentz factor, $\mbox{\boldmath$\beta$ } = \vec{v}/c$, and $\vec{E}(\vec{x},t)$, $\vec{B}(\vec{x},t)$ are the electromagnetic fields in which the particle is moving at position $\vec{r}$ and time t. Because a macroparticle represents N_q particles of elementary mass $m_{{\rm e}}$ and charge -e, the fraction q/m in the equation of motion is the same as in the case of single electrons. The electric field consist of two parts, $\vec{E}(\vec{x},t)= \vec{E}_{{\rm w}}(\vec{x},t)+ \vec{E}_{{\rm rad}}(\vec{x},t)$, where $\vec{E}_{{\rm w}}(\vec{x},t)$ is the wiggler field and $\vec{E}_{{\rm rad}}(\vec{x},t)$ is the radiation field, of which the development is studied in time; the same holds for $\vec{B}(\vec{x},t)$. The space charge effects of the particles in the beam are neglected. The initial conditions of the particles (i.e. number of charged particles, mass, charge, initial distributions in coordinate and velocity space) as well as the spatio-temporal behaviour (including initial amplitude and phase) of the wiggler, are user-specified.

To find the radiation in the simulations, a set of differential equations for the radiated energy are solved, additional to the differential Eqs. (15). To derive this set of differential equations, we first note that there are two essential differences between a FEL in the lab and our astrophysical application: firstly, there are no reflective mirrors in the neutron star magnetosphere. For a FEL to operate under such conditions, it should be a high-gain, single-pass process. Secondly, there are no side walls bounding the FEL cavity in the magnetosphere. This is accommodated in our pulsar study by using Gaussian modes of which the field decreases in a gaussian manner to zero away from the axis.

The radiation electric field is splitted into a set of Gaussian modes (see Appendix A):

 

$$\vec{E}_{{\rm rad}}(\vec{x},t)= \sum_j E_j(\vec{x},t)\hat{x}$$ (16)

$$E_j(\vec{x},t)= A_j T_j \cos~(\theta_j)$$ (17)

where

$$T_j = \frac{w_0}{w_j(z)}{\rm e}^{-r_\perp^2/w_j^2(z)}$$ (18)

$$\theta_j = \omega_j t - k_j z - \frac{k_j r_\perp^2 z}{2(z^2 + z^2_{0,j})} - \arctan\left(\frac{z}{z_{0,j}}\right) + \xi_j,$$ (19)

where we have chosen cylindrical coordinates (z, $r_\perp = (x^2 + y^2)^{1/2}$), and the polarisation is chosen to lie in the x-direction, in agreement with the polarisation of the wiggler which is expected to generate radiation polarised in the x,z-plane. Note that, there should be a small component of electric field in the z-direction as well. But since this is negligible, as compared to the x-component of the Gaussian mode, its interaction with the particles is neglected in the calculations.

Further, $w_j(z) = w_0 \sqrt{1 + z^2/z_{0,j}^2}$ is the waist, w₀ is the waist at z=0, and z_0,j = (1/2)k_jw₀² is the characteristic distance where the wave starts to diverge, A_j is the amplitude of this wave at wavevector k_j, $\omega_j = k_jc$ is the corresponding angular wavefrequency, and $\xi_i$ is an arbitrary phase.

Note that Gaussian modes propagate initially only in one direction and, therefore, they are suitable to describe narrow beams of pulsar radio emission. Figure 2 shows the wavefront of a Gaussian mode (dotted line); the solid lines are level curves at $\mbox{e}^{-1,-2,..}$ of the maximum amplitude (left: near-field, right: far-field).

Figure 2: Wavefronts (dotted lines) and level curves at e-n of the maximum amplitude of the Gaussian mode (solid lines, n = 1, 2, ...); top: near-field, bottom: far-field. z is in units of z0 and $r_\perp$ in units of w0.

As the figure shows, the wavefront of a Gaussian mode changes from planar to spherical. The asymptotic opening angle $\phi$ of this wave is given by:

$$\tan{\phi} = \sqrt{\frac{\lambda}{2 \pi z_0}}\cdot$$ (20)

In the GPT, the user specifies the minimum and the maximum frequency of the Gaussian modes and the number of modes in this range. The frequency intervals are equally spaced (given by  $\Delta N_m$) and each is represented by a Gaussian mode with wavevector k. The spectrum in this range is represented by the sum of these modes (Eq. (17)). When the radiation of the pulsar is assumed to be due to a FEL, then this set of Gaussian modes describes the radiation properly.

To continue the derivation of additional differential equations, we use conservation of energy: the energy gain/loss per unit time of the ith macro-particle with charge Q_i = N_q q_i due to the interaction with the jth electric field component is given by: $N_q q_i\:\vec{v}_i \cdot \vec{E}_j$, which is balanced by the rate of change of the electromagnetic fields. The total power radiated $P = \mbox{d}W_{{\rm light}}/\mbox{d}t$ is obtained from:

$$\sum_{i,j}N_q q_i\:\vec{v}_i \cdot \vec{E}_j + \frac{\mbox{d}W_{{\rm light}}}{\mbox{d}t} =0$$ (21)

where the summation is over all N particles and N_m Gaussian modes.

Instead of calculating the change in amplitude and phase of each mode directly, GPT rewrites the amplitudes and phases into two other independent variables m_j and n_j:

 

$$\sqrt{m_j^2 + n_j^2}/\Delta N_m$$                   A_j =

$$\phi_j \arctan(n_j/m_j).$$ = (22)

Inserting Eqs. (22) into the conservation of energy gives two differential equations for m_j and n_j (j = 1.. N_m):

 

$$\frac{\mbox{d}{m_j}}{\mbox{d}t} -\sum_i \frac{N_q q_i \Delta N_m}{\pi w_0^2 \epsilon_0 L} v_{x,i} T_j \cos~(\theta_j)$$ =

$$\frac{\mbox{d}{n_j}}{\mbox{d}t} \sum_i \frac{N_q q_i \Delta N_m}{\pi w_0^2 \epsilon_0 L} v_{x,i} T_j \cos~(\theta_j)$$ = (23)

where L is the cavity length.

A factor $\Delta N_m$ enters Eqs. (22) and (23) because one mode represents a frequency interval. This gives the correct calculation of m_j and n_j, which is dependent on $\Delta N_m$, and ensures that the resulting spectrum, represented by A_j's, is independent on the number of modes used.

 

 

Table 1: Input parameters and results for all simulation runs. Input parameters are: electron initial Lorentz factor $\gamma _0$, number density in the beam n in terms of Goldreich-Julian density at the pulsar surface $n_{{\rm GJ}\ast }$, dimensionless wiggler strength K, wiggler wavelength $\lambda _{{\rm w}}$, dimensionless wiggler phase speed $\beta _{{\rm w}}$ and (uniform) background magnetic field B₀. Results from the numerical simulations are: central frequency $\nu _{{\rm cen}}$, FWHM bandwidth of the radiation $\Delta \nu$, peak power of the light pulse  $P_{{\rm max}}$, and brightness temperature corresponding to the peak power  $T_{{\rm b}}$. The last four entries are the energy density of each ingredient in the simulation: beam $\epsilon _{{\rm b}}$, wiggler $\epsilon _{{\rm w}}$, background magnetic field $\epsilon _{{\rm b0}}$ and pulse $\epsilon _{{\rm p}}$. Note that every time "-'' appears, it means the entry has the value as in run 1. In run 5, n appears for $\nu _{{\rm res}}$, $\Delta \nu$ and $T_{{\rm b}}$, which means, that for those cases, the spectra are flat, no peaks are seen, and also the brightness temperatures which depend on the bandwidths are not calculated.

Run	Input parameters						Output				Energy densities
	$\gamma$	n	K	$\lambda _{{\rm w}}$	$\beta _{{\rm w}}$	B0	$\nu _{{\rm cen}}$	$\Delta \nu$	$P_{{\rm max}}$	$T_{{\rm b}}$	$\epsilon _{{\rm b}}$	$\epsilon _{{\rm w}}$	$\epsilon _{{\rm b0}}$	$\epsilon _{{\rm p}}$
		( $n_{{\rm GJ}\ast }$)		(m)		(T)	(GHz)	(GHz)	( $10^{13}\:\mbox{W}$)	( $10^{30}\:\mbox{K}$)	(106 Jm-3)	(Jm-3)	(Jm-3)	(102 Jm-3)
1	1000	0.23	10	50	0.9	0	11	2.42	9.09	1.02	1.30	8.15	0	8.25
2	-	-	5	-	-	-	42.8	9.72	36.6	1.04	-	2.04	-	34.5
	-	-	10	-	-	-	11	2.42	9.09	1.02	-	8.15	-	8.25
	-	-	15	-	-	-	4.96	1.02	3.67	1.16	-	8.35	-	3.46
	-	-	20	-	-	-	2.84	0.53	1.61	0.99	-	32.62	-	1.50
3	-	-	-	25	-	-	22.2	3.8	9.52	1.44	-	32.6	-	9.00
	-	-	-	40	-	-	13.8	2.86	9.42	1.02	-	10.07	-	8.89
	-	-	-	50	-	-	11	2.42	9.09	1.34	-	8.15	-	8.25
	-	-	-	100	-	-	5.57	1.28	7.02	1.70	-	2.03	-	6.63
4	-	-	-	-	0.2	-	88.5	14.1	9.59	0.215	-	94.9	-	9.00
	-	-	-	-	0.5	-	55.7	8.12	9.54	0.37	-	18.2	-	9.00
	-	-	-	-	0.8	-	22.2	3.79	9.39	0.79	-	9.35	-	8.86
	-	-	-	-	0.9	-	11	2.42	9.09	1.02	-	8.15	-	8.25
5	-	-	-	-	-	10-3	10.9	2.43	9.46	7.4	-	-	0.4	8.85
	-	-	-	-	-	0.01	2.18	0.61	2.7	1.36	-	-	40	2.55
	-	-	-	-	-	0.025	n	n	0.455	n	-	-	$2.48 \times 10^2$	$42.8 \times 10^{-2}$
	-	-	-	-	-	0.05	n	n	$6.36 \times 10^{-4}$	n	-	-	$9.95 \times 10^{2}$	$4.0 \times 10^{-4}$
6	-	0.1	-	-	-	-	11.37	1.53	1.79	0.38	2.30	-	-	1.7
	-	0.2	-	-	-	-	11.1	2.22	7.2	1.07	1.76	-	-	6.8
	-	0.3	-	-	-	-	10.9	3.01	16.5	1.73	1.15	-	-	15.5
	-	0.4	-	-	-	-	10.7	3.40	29.9	2.8	0.575	-	-	28.2
7	1000	-	-	-	-	-	11	2.42	9.09	1.24	1.3	-	-	8.24
	500	-	-	-	-	-	2.72	0.80	1.77	0.72	0.16	-	-	1.66
	250	-	-	-	-	-	0.69	0.16	0.112	0.23	0.02	-	-	0.10
	100	-	-	-	-	-	0.09	0.006	$2.7\times 10^{-4}$	$4.6\times 10^{-3}$	$1.36\times 10^{-3}$	-	-	$2.5 \times 10^{-4}$

The wiggler part of the electromagnetic field is assumed to have a time-independent amplitude during passage of the beam, and is not included here, for we assume that the change in the wigglers' energy is small on the considered timescales.

Summarising, the differential equations solved at time t by GPT consist of the 2N equations:

$$\frac{\mbox{d}\gamma_i \mbox{\boldmath$\beta$ }_i}{\mbox{d}t} = \frac{q_i}{m_ic}[\vec{E}(\vec{x},t)+ \vec{v}_i(t) \times \vec{B}(\vec{x},t)]$$

$$\frac{\mbox{d}\vec{x}_i}{\mbox{d}t} = \vec{v}_i = \frac{\gamma_i\mbox{\boldmath$\beta$ }_i c}{\sqrt{\gamma_i^2\mbox{\boldmath$\beta$ }_i^2 + 1}}$$

for the N charged particles' phase space positions, and the set of 2N_m equations:

$$\frac{\mbox{d}{m_j}}{\mbox{d}t} = -\sum_i \frac{N_q q_i \Delta N_m}{\pi w_0^2 \epsilon_0 L} v_{x,i} T_k \cos~(\theta_j)$$

$$\frac{\mbox{d}{n_j}}{\mbox{d}t} = \sum_i \frac{N_q q_i \Delta N_m}{\pi w_0^2 \epsilon_0 L} v_{x,i} T_k \cos~(\theta_j)$$

for the total number of modes N_m included in the simulation, to calculate the amplitude and phase at each wavenumber in the specified range during the interaction with the particles. As a result, we obtain the spectrum of the radiation in time.

To retrieve the pulse in the time domain, both the electric field and the power are Fourier transformed from the frequency to the time domain, and Parseval's theorem is applied to the results to ensure a correct transformation between time and frequency domain.

Note that the results obtained with this code are independent of the number of time outputs, the number of macroparticles and the number of modes.

4 Physical parameters for numerical simulations

Except in the second last run of the simulations, we used a beam of 200 macroparticles with a total charge of $-6 \times 10^{-4}\:\mbox{C}$. These are uniformly distributed in a cylinder with radius of one metre and length of three ponderomotive wavelengths (8), which varies from run to run. The corresponding electron number density is $1.59 \times 10^{16}\ \mbox{m}^{-3}$. This is 0.23 times the Goldreich-Julian density at the stellar surface, given by $n_{{\rm GJ}}(r_\ast) = 2 \epsilon_0 \mbox{\boldmath$\Omega$ }\cdot\vec{B}_\ast/q$, for a reference pulsar with $\Omega = 2 \pi/ P$, where $P = 1\:\mbox{s}$ and $\vec{B}_\ast = 10^8~\mbox{T}$. Except in the last run, each macro-particle's initial Lorentz factor is $\gamma = 1000$.

The wiggler is an electromagnetic disturbance propagating in the z-direction. In the pulsar magnetosphere, the scale at which the physical parameters change, e.g. background magnetic field, is of the order of $z/r_\ast$ ($\geq$ $10\:\mbox{km}$). Therefore, the interaction between the beam and the wiggler is simulated until the beam particles reach this distance, so $t_{{\rm end}} \approx 3.3\times 10^{-5}\:\mbox{s}$. The parameters for the wiggler are given in Table 1 (all labelled with "$\mbox{w}$''). The values for the magnetic induction and angular frequency are obtained from (2) and (5).

In the first run, the wiggler's phase speed is 0.9c. The amplitude of the magnetic induction $B_{{\rm w}}$ is $3 \times 10^{-3}\:\mbox{T}$. The wavelength $\lambda _{{\rm w}}$ is 50 metres. These values are chosen such that the resonance is at radio frequencies, and the strength of the magnetic induction is chosen such that the induced transverse velocity is small compared to the initial axial velocity ( $K/\gamma \ll 1$ in Eq. (4)).

In the subsequent runs we investigate the effects of varying each of these parameters as compared to the results of the first run. In run 2, the wiggler strength K is increased from five to twenty, in steps of K=5 (correspondingly the magnetic induction is increased in steps of $B_{{\rm w}} = 1.5 \times 10^{-3}~\mbox{T}$). In the third run, the wavelength of the wiggler takes the values 25, 40, 50 and 100 m. In run 4, the velocity of the wave is changed from 0.2 to 0.9.

In run 5, we study the effect of a non-zero background magnetic field. We took the most simple case is which the field in uniform (in z) and has a background magnetic field strength of B₀ = 10^-3, 10^-2, $2.5 \times 10^{-2}$, $5 \times 10^{-2}$ T.

In run 6 as we change the number density n of the electron beam. In terms of the Goldreich-Julian density, the fractional density ranges from 0.1 to 0.4.

Finally, in run 7, we lower the Lorentz factor of the beam particles (the resonance frequency depends on this as $\omega_{{\rm res}} \propto \gamma^2$). The Lorentz factors are chosen: $\gamma = 1000$, 500, 250, 100.

5 Results

Figures 3a(Run 1) to 9c(Run 7) and Table 1 show the end results of all runs, i.e. at $t = 3.3 \times 10^{-5}\:\mbox{s}$. Rather than discussing the runs sequentially, we present the systematic trends of our computations, using the figures as illustrations. Run 5, where a uniform background field is included, is presented separately.

Figure 3a: Run 1: (z, x)-behaviour of the particles during their passage in the cavity (from left to rigth; from top to bottom). On the horizontal axis is plotted the z-coordinate of all particles, subsequently for z = 2, 4, 6, 8 and 10 km. Except for $t = 3.34 \times 10^{-5}\:\mbox{s}$, the horizontal spacing between the long ticks is 0.02 m (e.g. $\Delta z \approx 0.08 \:\mbox{m}$ for $t = 0\:\mbox{s}$). In the last plot, the spacing between the long ticks is 0.05 m. Bunching clearly occurs at $t = 1.33\times 10^{-5}\;\mbox{s}$, "de-bunching'' at $t = 2.0 \times 10^{-5}\:\mbox{s}$ and "re-bunching'' at $t = 2.67 \times 10^{-5}\:\mbox{s}$.

Figure 3b: Run 1: Power of the pulse at the end of the simulation $t = 3.3 \times 10^{-5}\:\mbox{s}$ ( top: number of macroparticles N=200; bottom: N=2). On the horizontal axis is also plotted t, this is not the simulation time, but the duration of the pulse as a distant observer would measure. Negative t corresponds to earlier arrival. The maximum power for N=2 ( bottom) is $P_{2,{\rm max}} = 9.75 \times 10^{9}\:\mbox{W}$. Since $P_{200} \propto N^2 P_{2}$, this implies the radiation pulse at the top is coherent.

Figure 3c: Run 1: spectrum at $t = 3.3 \times 10^{-5}\:\mbox{s}$. The central frequency is at $\nu = 11\:\mbox{GHz}$, with FWHM of $\Delta \nu = 2.42\:\mbox{GHz}$. The peak on the right, centered on $\nu_{{\rm res}} = 11.9\:\mbox{GHz}$, is caused by the first bunching at $t = 1.33 \times 10^{-5}\:\mbox{s}$. And the peak on the left corresponds to the second buching, where $\gamma$ dropped $5\%$, resulting in a lower $\nu _{{\rm res}}$.

Figure 3d: Run 1: time development of average Lorentz factor (i.e. per macro-particle) in two simulations, one with N=2 ( top) and the other with N=200 ( bottom). The N=2 case reflects the behaviour of particles radiating incoherently, with constant energy loss over time as compared with the bunching case, where there is a steep drop in energy ( $\Delta \bar{\gamma} = 50$) of the beam particles between 1.0 and $1.8 \times 10^{-3}\:\mbox{s}$. Between $2.4 \times 10^{-5}\:\mbox{s}$ and $2.8 \times 10^{-5}\:\mbox{s}$ another drop in beam energy occurs, which corresponds to a second bunching of the particles. This time, the beam energy loss is less than for the first time, and the duration is much shorter. The pulse power (Fig. 3b) reflects these characterics where the first pulse is more powerful and has a longer duration than the second pulse.

Figure 4a: Spectra for Run 2 plotted over each other; we have K = 5 (solid), 10 (dotted), 15 (short-dashed) and 20 (long-dashed). The resonance frequency shifts according to Eq. (7). The bandwidth increases with increasing $\nu _{{\rm res}}$.

Figure 4b: Pulse power for run 2 plotted for different K. The maximum scales roughly as K-2.

Figure 4c: Run 2: ( $\bar{\gamma}$, t)-plot for K = 5, 10, 15, 20 (labelling as in Fig. 4a). Obviously, for all K, the beam particles bunch twice (also, see Fig. 4b). Corresponding to the maximum power, the beam energy loss is largest here for K = 5.

5.1 Bunching

We found bunching of particles during their interaction with the wiggler in a number of cases. Run 1 clearly shows this bunching (Fig. 3a). Here, the particles are plotted in (z, x)-projection at different stages of their interaction: bunch formation ( $t = 1.33 \times 10^{-5}\:\mbox{s}$), "de-bunching'' ( $t = 2.0 \times 10^{-5}\:\mbox{s}$) and "re-bunching'' ( $t = 2.67 \times 10^{-5}\:\mbox{s}$). Note that the distance between bunches is $\lambda_{{\rm p}}$ as expected.

An alternative illustration of the formation of bunches is found in Fig. 3d, where the average Lorentz factor per macro-particle  $\bar{\gamma}$ is plotted. Evidently, the formation of bunches corresponds to a steep drop of  $\bar{\gamma}$ (e.g. Fig. 3d, the first bunch formation starts from $t \simeq 1.0\times 10^{-5}\:\mbox{s}$ to $1.8 \times 10^{-5}\:\mbox{s}$, and the second formation from $2.4 \times 10^{-5}\:\mbox{s}$ to $2.8 \times 10^{-5}\:\mbox{s}$).

We notice that both the bunch duration and the beam energy loss is larger for the first bunch formation. Saturation is reached due to an increased velocity spread. Therefore, for the second bunch formation, the beam starts with a certain spread in velocity, and it reaches saturation quicker than in the initial bunching.

For the beam number density of $n = 0.23\:n_{{\rm GJ}\ast}$, we find a beam energy loss to radiation of $\sim$$5\%$, after the first time bunching. This number is only affected by n. The number of particles in each bunch scales with n, as does $\Delta \bar{\gamma}$ (Run 6, Fig. 8c), which agrees with coherent losses.

As to the starting time of the first bunching, we find the following: for the same beam and different wiggler parameters ( $\beta _{{\rm w}}$ and $\lambda _{{\rm w}}$), about ten oscillations of the particles in the wave as seen by the observer are needed before they start to bunch (Run 3, Fig. 5c and Run 4, Fig. 6c). For the same wiggler, increasing the beam density or decreasing the Lorentz factor give rise to earlier occurence of bunching (Run 6, Fig. 8c and Run 7, Fig. 9c). This suggests that the bunching sets in when the (incoherent) radiation reaches a certain level, because the ponderomotive force is dependent on the radiation field.

As for the wiggler strength K, it determines the maximum beam energy loss during bunching. In Fig. 4c(Run 2) are shown ( $\bar{\gamma}$, t)-plots for K = 5, 10, 15, 20. We have also run cases for $K \leq 1$. The maximum energy loss occurs when K = 1.

5.2 Coherent radiation

That the radiation is coherent is demonstrated by the power of the radiation pulse (Fig. 3b). At the top is presented the simulation run for N = 200macro-particles and at the bottom for N = 2 macro-particles. The maximum pulse power for N = 200 is $P_{200} = 9.09 \times 10^{13}\:\mbox{W}$ and for N = 2, $P_{2} = 9.75 \times 10^9\:\mbox{W}$. If the radiation were incoherent emission from 200 macro-particles, we would expect $P \propto N$, where N is the number of particles and $(P_{200}/P_2)_{{\rm inc}} = 100$. However, the simulation shows that $(P_{200}/P_2) \simeq 9 \times 10^3$, and thus $P \propto N^2$, demonstrating that the radiation is coherent.

Note that only the maximum of the first peak satisfies $P_{200} \propto N^2$. This is not the case for the second peak in Fig. 3b, where now $P_{200} \approx 2 \times 10^{13}\:\mbox{W}$, which is less than the coherent case but much more than if it were incoherent.

Figure 5a: Run 3: spectrum at the end of the simulation for different $\lambda _{{\rm w}} = 100$ m (solid), 50 m (dotted), 40 m (short-dashed), 25 m (long-dashed). Spectral broadening is due to more bunchings occuring during the simulation.

Figure 5b: Pulse power for Run 3 plotted for different $\lambda _{{\rm w}}$ values. Maximum powers are of about the same order of magnitude. Only the number of pulses differs for each run, due to number of bunching times.

Figure 5c: Run 3: the average Lorentz factor of the macro-particles during their passage in the cavity for $\lambda _{{\rm w}} = 100$, 50, 40, 25 m (labelling as in Fig. 5a). The starting time of the first bunching for each run is about $10 \:t_{{\rm eff}}$ (6). The energy beam loss after the first bunching is about $5\%$.

5.2.1 Spectrum

The radiation frequency in our simulations agrees well with the expected resonance condition (compare Eq. (7) and Table 1). Therefore, by construction, the frequency of the emission is in the radio regime, between $1 \sim 10 \:\mbox{GHz}$. In some cases, the central frequency is lower than the resonance condition (Runs 2-4: Figs. 4a, 5a, 6a). The shift in these cases is caused by multiple bunching where the average Lorentz factor decreases after each bunching, and therefore, also the frequency as determined by the resonance condition. As a result of multiple bunchings, the total bandwidth increases (Runs 3, 4, 6 in Figs. 5a, 6a, 8a), whereas the relative bandwidth $\Delta \nu/\nu$ decreases.

5.2.2 Pulse shape, power, brightness temperature and radiation cone

The pulse shapes of the coherent radiation can be found in Figs. 3b, 4b, 5b, 6b, 8b and 9b.

The number of temporal peaks corresponds to the number of bunching times (not the number of bunches that is formed); e.g. in Run 1 (N=200) the beam particles bunched twice (Fig. 3d), and the pulse has two peaks (Fig. 3b). The maximum pulse power always occurs at the first time bunching. The following bunch formations result in less powerful pulses (though not completely incoherent). The maxima of the pulses have the same order of magnitude for most runs (Table 1), except for $\gamma = 100$, where the coherent emission becomes less efficient (and also for $B_0 = 0.05\:\mbox{T}$, but this is discussed in the next section). The optimum efficiency of the emission power is reached for $\gamma \simeq 300$ as follows from the ratio $\epsilon_{{\rm p}}/\epsilon_{{\rm b}}$ (Table 1, Run 7: Fig. 9b)

The structure within each pulse reflects the wiggling motion of the particles (e.g. Fig. 5b, where $\lambda_{{\rm w}} = 100\:\mbox{m}$).

The typical duration of the pulse is a few nanoseconds.

The brightness temperature $T_{{\rm b}}$ is derived using

$$I = \frac{k_{{\rm B}}T_{{\rm b}}\nu^2}{c^2} = \frac{W/\Delta t}{A \Delta\nu \Delta\Omega}$$ (24)

where I is the radiation intensity, $k_{{\rm B}}$ is Boltzmann's constant, $W/\Delta t$ is the power, which we take as $P_{{\rm max}}$, A is the emitting surface, which corresponds to the beam cross-section, $\Delta \nu$ is the bandwidth, and $\Delta \Omega = \pi \tan^2\phi$ the solid angle into which the radiaton is emitted (20).

Figure 6a: Run 4: spectrum for $\beta _{{\rm w}} = 0.2$ (solid), 0.5 (dotted), 0.8 (short-dashed), 0.9 (long-dashed).

Figure 6b: Pulse power for run 4 plotted for different $\beta _{{\rm w}}$.

Figure 6c: Run 4: the average Lorentz factor versus the average distance z of the particles in simulations with changing $\beta _{{\rm w}} = 0.2$ , 0.5, 0.8, 0.9 (labelling same as Fig. 6a). The starting time of the first bunching is  $10 \:t_{{\rm eff}}$.

Converting the maximum power into the brightness temperature, we find that, again, for most cases the brightness temperature is approximately $10^{30}\:\rm K$ (Table 1).

5.3 Background magnetic field

An ambient magnetic field in the z-direction (Run 5) can have a drastic effect on the bunching (Fig. 7c). For a background field just smaller than the wiggler field, $B_0 = 10^{-3}\:\mbox{T}$, the particles' behaviour shows no difference with the simulations where the background field is absent. For $B_0 = 10^{-2}\:\mbox{T}$, we still observe bunching. For larger values of B₀, however, bunching disappears completely. Actually, such behaviour is expected from Eq. (13): since the motion of the particles is no longer free in the transverse direction, the ponderomotive force is less effective and the resulting bunching is less pronounced. For all runs, the angle between the particles' velocity and the total magnetic field $\vec{B} = B_0 + B_{{\rm w}}$, which is proportional to $\vec{v} \cdot \vec{B}$ is indeed found to be nonzero.

Since only the cases where $B_0 \leq 10^{-2}\:\mbox{T}$ show bunching, these are discussed further. The spectrum for this run is plotted in Fig. 7a. Only for $B_0 = 10^{-3}\:\mbox{T}$, the characteristics of the spectrum are similar to the cases where $B_0 = 0 \:\mbox{T}$. For $B_0 = 10^{-2}\:\mbox{T}$, the resonance frequency shifts to $\nu_{\rm res} = 2.18 \:\mbox{GHz}$ (as compared to the expected value of $11.9\:\mbox{GHz}$). The bandwidth is $0.61\:\mbox{GHz}$, and is much smaller than before.

Figure 7a: Run 5: spectra for (from top to bottom along the left vertical axis) B0 = 10-3 T (solid), 10-2 T (dotted), 0.0025 T (short-dashed), 0.050 T (long-dashed). For $B_0 = 10^{-3}\:\mbox{T}$, the spectrum shows the same properties (central frequency and bandwidth) as in the case where the background magnetic field is absent. For $B_0 = 10^{-2}\:\mbox{T}$, the central frequency shifts to $2\:\mbox{GHz}$, whereas the bandwidth drops to $0.61\:\mbox{GHz}$.

Figure 7b: Pulse power (logarithmic) for Run 5. The power of the radiation becomes smaller for B0 = 0.025 and 0.05 T.

Figure 7c: Run 5: the average Lorentz factor versus the average distance z of the particles in simulations with background magnetic field for $B_0 = 10^{-3}, 10^{-2}, 2.5 \times 10^{-2}, 5.0 \times 10^{-2}\:\mbox{T}$ (labelling as in Fig. 7a). The last case coincides practically with the horizontal line $\bar{\gamma} = 1000$. Only for $B_0 = 10^{-3}\:\mbox{T}$ and $B_0 = 10^{-2}\:\mbox{T}$ the beam particles show FEL action.

Figure 8a: Run 6: spectrum (from top to bottom) for $n/n_{{\rm GJ}\ast } = 0.1$ (solid), 0.2 (dotted), 0.3 (short-dashed) and 0.4 (long-dashed). The spectrum broadens as the number density increases, as expected from Fig. 8c, which shows single bunching for $n/n_{{\rm GJ}} = 0.1$, but multiple bunching for $n/n_{{\rm GJ}} > 0.1$.

Figure 8b: Pulse power for Run 6 plotted for different number densities relative to $n_{{\rm GJ}\ast }$. As expected, the maximum power scales as n2.

Figure 8c: Run 6: from top to bottom: Average Lorentz factor per particle, for number density of the bunch n = 0.1, 0.2, 0.3, $0.4\:n_{{\rm GJ}\ast}$ (labelling as in Fig. 8a). Clearly, a small beam number density results in negligible losses, and therefore, no coherent emission.

Figure 9a: Run 7: spectrum for initial Lorentz factor of the beam particles $\gamma = 1000$ (solid), 500 (dotted0), 250 (short-dashed), 100 (long-dashed). The central frequency shift as $\gamma ^2$.

Figure 9b: Pulse power for run 7 plotted.

Figure 9c: Run 7: the average Lorentz factor versus the average distance z of the particles in simulations with changing Lorentz factor.

6 Discussion

We have investigated a specific form of a single-pass free-electron laser process as a possible mechanism to produce high brightness radio emission of pulsars. We have investigated the operation of a FEL in the presence of a wiggler which consists of a transverse electromagnetic disturbance as it is being overtaken by a relativistic electron beam.

We have shown that, in principle, a FEL can operate in the pulsar outer magnetosphere in the presence of a transverse wiggler. The deciding factors for particles to bunch and emission to be coherent within $3.3 \times 10^{-5}\:\mbox{s}$, are the following: a large beam particle density $n \geq 0.1\:n_{{\rm GJ}}(r_\ast)$, a Lorentz factor of the beam particles $\gamma > 100$ for K = 10, and a small background magnetic field $B_0 \leq 10^{-2}\:\mbox{T}$. The brightness temperature of the pulse depends sensitively on these parameters (Table 1).

The required beam particle density together with the wiggler's parameters mainly determine the timescale on which particles start to bunch. Bunching occurs after about 10 times the transit time through the wiggler (i.e. 10  $t_{{\rm eff}}$, Eq. (6)). This seem to set the level of inchorent radiation, which is then large enough for the ponderomotive force to act on the beam particles. Then the timescale over which a bunch stays together is about one  $t_{{\rm eff}}$, i.e. coherent radiation only occurs during this period. Due to an increase in axial velocity spread, debunching occurs. Although particles rebunch, the associated pulse is weaker due to the velocity spread in the bunch. In the pulse profile, the first bunching results in one pulse. Then for larger times, particles bunch more often and the pulse acquires more peaks, which are less powerful.

The overall fractional energy losses of the beam are $\sim$$5\%$. Together with the travel time, the calculations demonstrate that when applied to the pulsar magnetosphere for $B \leq 10^{-2}\:\mbox{T}$, the FEL interaction can produce coherent radiation. For a dipole model, this implies $R \gg 2 \times 10^3 R_\ast$, where we used $B_\ast = 10^8 \:\mbox{T}$, i.e. in the outer magnetosphere.

While the particles are bunched, most radiation is at $\nu _{{\rm res}}$, which in our runs varies from $0.7\:\mbox{GHz}$ to $90\:\mbox{GHz}$ as in the pulsar radio emission. The lower frequency is obtained with $\gamma = 250$.

The particles lose their energy most efficiently for K=1.

The measured FWHM bandwidth for the radiation $\Delta \nu$ ranges from $0.53\:\mbox{GHz}$ to $14.5\:\mbox{GHz}$, and is broadband $\Delta \nu/\nu \simeq 0.07 {-} 0.5$. The width of the spectrum highly depends on the number of bunching events occurring during the simulation. After the first bunching, the bandwidth is of the order of $1.5\:\mbox{GHz}$; After the second -less effective- bunching, $\Delta \nu \approx 2.9\:\mbox{GHz}$ and so on. Of course, the central radiation frequency agrees with radio pulsar emission by construction. However, it is interesting that the bandwidth is found to agree with pulsar radio observations where average pulses and microstructures are observed from about $100\:\mbox{MHz}$ to more than $10\:\mbox{GHz}$ (average pulses: e.g. Lyne et al. (1998); D'Amico et al. (1998); Manchester et al. (1996); Lyne & Manchester (1988), microstructures: e.g. Bartel & Sieber (1978); Rickett et al. (1975)).

As for the brightness temperature $T_{{\rm b}}$, which is a big obstacle for most radiation processes, we found a brightness temperature $T_{{\rm b}}$ at pulse maximum of $10^{30}\:\mbox{K}$ for $\beta_{{\rm w}} = 0.9$and $T_{{\rm b}} = 10^{29}\:\mbox{K}$ for a wiggler phase velocity $\beta_{{\rm w}} \la 0.8$, similar to the observationally derived pulsar brightness temperatures.

Apparently, the FEL is able to produce the required high brightness pulse. Further, the characteristic opening angle, given by $\tan \phi \approx \lambda/(2 \pi w_0)$, is $\phi \la 1^\circ$. Again, this agrees with the observed values in pulsar radio emission, as estimated from the microstructure duration relative to the pulsar rotation period: $\sim$ $100~ \mu \mbox{s}/P \times 360^\circ \la 1^\circ$ for $P = 0.1 \:\mbox{s}$.

The radiation pulse has a duration of $2\:\mbox{ns}$. The shortest elements of radio emission measured from pulsars are microstructures. These are quasi-periodic structures of $\sim$ $10^2~ \mu\mbox{s}$, which are broadband and highly (linearly) polarized (Lange et al. 1998; Rickett et al. 1975; Popov et al. 2002; Cordes & Hankins 1979).

Since the starting time for the particles to bunch is $10\: t_{{\rm eff}} = 10\lambda_{{\rm w}}/(\beta_z - \beta_{{\rm w}})$, the faster the wiggler wave is, the longer it takes for the particles to get bunched. Also, the characteristic pulse duration  $t_{{\rm eff}}$increases. This is clearly shown in run 4, where we varied the phase velocity of the electromagnetic disturbance between 0.2c and 0.9c. The computation time limitation forced us to consider only such relatively low values for electromagnetic wigglers in the relativistic outflow from a pulsar magnetosphere, but one can see from the particles' behaviour in Fig. 6c, that even more relativistic electromagnetic disturbances would lead to longer timescales, e.g. for $\gamma_{{\rm w}} = (1 - \beta_{{\rm w}}^2)^{-1/2} = 100$, the timescale would go up by a factor 5000, and become comparable to the observed micropulse durations. (The same argument holds for a lower beam number density n, Fig. 8c; i.e. when $n < 0.1 n_{{\rm GJ}\ast}$, bunching occur at $t > 3.3 \times 10^{-5}\:\mbox{s}$).

Finally, we note that the coherent emission from electrons and positrons in a transverse wiggler add constructively. Therefore, the total number of electrons and positrons in a bunch, and not the charge excess, determines the emission.

In summary, the operation of a single-pass high-gain FEL with a transverse electromagnetic wiggler within the pulsar magnetosphere in the radio regime, requires a mono-energetic beam of electrons/positrons at moderate Lorentz factors in a sufficiently small background magnetic field.

We expect a FEL process to be also possible much nearer to the pulsar, where the magnetic field strength is large and the dynamics of the particles are one-dimensional, when longitudinal instead of transverse wigglers are used. This would be the domain of single-pass, high-gain Cerenkov FELs, and a next logical step to study. Strong Langmuir turbulence has been studied in this context by (Schopper et al. 2002). By using a Particle-in-Cell method, these authors show that electron scattering on Langmuir turbulence, which is excited in a self-consistent way, results in high power output of radiation.

Acknowledgements

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). We want to thank the anonymous referee for his/her comments. J.K. likes to thank Dr. Marnix van der Wiel and Dr. Theo Schep for stimulating interest in the project. P.K.F. likes to thank Bas van de Geer, Marieke de Loos for their support in using the code GPT and Cyrille Thomas for the description of the Gaussian modes. Also, we want to thank Kees van der Geer for setting up the simulations.

References

 
• Asseo, E., & Melikidze, G. I. 1998, MNRAS, 301, 59
   
• Bartel, N., & Sieber, W. 1978, A&A, 70, 307
   
• Cordes, J. M., & Hankins, T. H. 1979, ApJ, 233, 981
   
• D'Amico, N., Stappers, B. W., Bailes, M., et al. 1998, MNRAS, 297, 28
   
• Erikson, W. L., & Singh, S. 1994, Phys. Rev. E, 49, 5778
   
• Lange, C., Kramer, M., Wielebinski, R., & Jessner, A. 1998, A&A, 332, 111
   
• Luo, Q., & Melrose, D. B. 1995, MNRAS, 276, 372
   
• Lyne, A. G., & Manchester, R. N. 1988, MNRAS, 234, 477
   
• Lyne, A. G., Manchester, R. N., Lorimer, D. R., et al. 1998, MNRAS, 295, 743
   
• Lyutikov, M. Y. 1998, Ph.D. Thesis, California Institute of Technology, Pasadena, California, USA
   
• Manchester, R. N., Lyne, A. G., D'Amico, N., et al. 1996, MNRAS, 279, 1235
   
• Melikidze, G. I., Gil, J. A., & Pataraya, A. D. 2000, ApJ, 544, 1081
   
• Melrose, D. B. 1978, ApJ, 225, 557
   
• Melrose, D. B., & Gedalin, M. E. 1999, ApJ, 521, 351
   
• Milonni, P. W., & Eberly, J. H. 1988, Lasers, 1st ed. (New York: John Wiley and Sons), 480
   
• Popov, M. V., Bartel, N., Cannon, W. H., et al. 2002, A&A, 396, 171
   
• Press, W. H., Flannery, B. P., Teukolsky, S. A., & Vetterling, W. T. 1992, Numerical Recipes in C: The Art of Scientific Computing, 2nd ed. (Cambridge University Press, 2nd edition)
   
• Rickett, B. J., Hankins, T. H., & Cordes, J. M. 1975, ApJ, 201, 425
   
• Schopper, R., Ruhl, H., Kunzl, T. A., & Lesch, H. 2002, in Neutron Stars, Pulsars, and Supernova Remnants, ed. W. Becker, H. Lesch, & J. Trümper, 193
   
• Usov, V. V. 1987, ApJ, 320, 333
  

 

  
7 Online Material

Appendix A: Derivation of Gaussian modes

The full derivation of the Gaussian modes can be found in Milonni & Eberly (1988) and Erikson & Singh (1994). As the radiation is produced by a narrow beam of relativistic particles the radiation propagates mainly in one direction (z-axis). As to the dependence of the radiation in the (x, y)-plane, we assume that, similar to a laser beam, the intensity has a Gaussian form inside the FEL; i.e. $I_\perp \propto {\rm e}^{-(r_\perp/w_0)^2}$ where $r_\perp = \sqrt{x^2 + y^2}$ and w₀ denotes the transverse distance where the intensity drops to 1/e of the peak value at z=0 (w₀ is called the waist). The electromagnetic field which has these properties and satisfies Maxwell's equations is the Gaussian mode. The derivation of Gaussian modes starts from Maxwell's equations in vacuum, which result in the (vector) wave equation:

$$\nabla^2 \vec{E}(\vec{r}, t) - \frac{1}{c^2}\frac{\partial^2 \vec{E}(\vec{r}, t)}{\partial t^2} = 0 .$$ (A.1)

For vacuum waves ( $\omega^2 = k^2c^2$) propagating mainly along the z-direction, the paraxial approximations can be applied:

$$\left\vert \frac{\partial^2 \vec{E}(\vec{r})}{\partial z^2}\right\vert \ll 2 k \left\vert \frac{\partial \vec{E}(\vec{r})}{\partial z}\right\vert$$ (A.2)

$$\left\vert \frac{\partial^2 \vec{E}(\vec{r})}{\partial z^2}\right\vert \ll 2k \left\vert \vec{E} (\vec{r})\right\vert .$$ (A.3)

With these, the vector wave equation becomes the vector paraxial wave equation:

$$\nabla^2_{\textrm{T}} \vec{E}(\vec{r}) + 2{\rm i}k \frac{\partial \vec{E}(\vec{r})}{\partial z} = 0$$ (A.4)

where $\nabla^2_{\textrm{T}} = \partial^2/\partial x^2 + \partial^2/\partial y^2$. This equation is valid for each component of $\vec{E}(\vec{r}) = E_x(\vec{r})\hat{x} + E_y(\vec{r})\hat{y} + E_z(\vec{r})\hat{z}$. The same holds for $\vec{B}(\vec{r})$.

   

Table A.1: The range of opening angles $\phi$ for different z₀ and wavelengths between 3 cm and 30 m according to Eq. (A.20).

z0	$\phi$
Polar cap radius Rpc = 100 m	24$^\prime$-3$^\circ$57$^\prime$
Stellar radius $R_{\ast}$ = 104 m	2$^\prime$23 $^{\prime\prime}$-24$^\prime$
Light cilinder Rlc = 107 m	4 $^{\prime\prime}$.5-45 $^{\prime\prime}$

The general solution for each component of this equation, which is axisymmetric, is the following:

$$u(r_\perp, z) = a \exp(-P(z)) \exp\left(-\frac{{\rm i}k r_\perp^2}{2 Q(z)}\right)$$ (A.5)

where a is an amplitude and P(z) and Q(z) are complex functions which specify the longitudinal and transverse mode behaviour. These can be retrieved by putting Eqs. (A.5) into (A.4), which results in two differential equations for P(z) and Q(z):

$$\frac{\mbox{d}Q}{\mbox{d}z}$$ (A.6)

$$\frac{\mbox{d}P}{\mbox{d}z} \frac{{\rm i}}{Q}\cdot$$ (A.7)

The solutions are simply given by:

Q(z) = z + q₀ (A.8)

$$-{\rm i} \ln (z + q_0)$$ (A.9)

where q₀ is an integration constant. Recall that function Q(z) gives the transverse behaviour of the mode and that the solution is paraxial. Therefore, Q(z) can be written in terms of a radius of curvature R(z) and a width w(z):

$$\frac{1}{Q(z)} = \frac{1}{z + q_0} = \frac{1}{R(z)} + \frac{-2{\rm i}}{k w^2(z)}\cdot$$ (A.10)

Furthermore, we assume that at a reference point z=0, the mode wavefront curvature is $R(0) \equiv \infty$, so that:

$$q_0 = \frac{{\rm i}kw^2(0)}{2}\cdot$$ (A.11)

Together with:

$$P(z) = -{\rm i}\ln(z + {\rm i} z_0) = -{\rm i}\ln\left[\left(z^2 + z_0^2\right){\rm e}^{{\rm i}\phi}\right]$$ (A.12)

$$\qquad = -{\rm i}\ln\left[z^2 + z_0^2\right] + \arctan(z_0/z)$$ (A.13)

$$z_0 = \frac{k w_0^2}{2}$$ (A.14)

$$\tan\phi = \frac{z_0}{z}$$ (A.15)

the general solution becomes:

$$u(r_\perp, z) = A \frac{w_0}{w(z)} Re\left[\exp\left(-r_\perp^2/w... ... kz + {\rm i} \omega t- {\rm i} k r_\perp^2/2 R(z) - {\rm i} \phi\right)\right]$$ (A.16)

where

$$w_0 \sqrt{1 + \left(\frac{\lambda z}{\pi w_0^2}\right)^2}$$                     w(z) = (A.17)

$$z + \frac{1}{z}\left( \frac{\pi w_0^2}{\lambda}\right)^2$$ R(z) = (A.18)

$$\tan\phi \frac{\lambda z}{\pi w_0^2}$$ = (A.19)

where $r_\perp$ and w₀ are as defined above, w(z) is the waist at z and $A, k, \omega$ are the amplitude, the z-component of the wave vector and the frequency of the electromagnetic wave respectively. The curvature radius of the wavefront is given by R(z) = z + z₀²/z, where $z_0 = (\pi w_0^2)/\lambda$ is (roughly) the separation between the near and the far field. The wavefronts of this mode change from planar in the near field to spherical in the far field (see Fig. 2).

The opening angle $\phi$ of the intensity of the mode is given by:

$$\tan{\phi} = \frac{1}{2}\sqrt{2}\frac{w(z)}{z} = \sqrt{\frac{\lambda}{2 \pi z_0}}\cdot$$ (A.20)

Within distance of length z₀, wavefronts of the Gaussian modes are considered planar. Different wavelengths and different z₀result in different opening angles. In table A.1 is listed the range of opening angles when we consider radio waves with wavelengths between 3 cm to 30 m and some typical distances in the pulsar magnetosphere.