K. Asano - F. Takahara
Department of Earth and Space Science, Osaka University, Toyonaka 560-0043, Japan
Received 17 September 2003 / Accepted 28 November 2003
Motivated by the particle acceleration problem in pulsars, we numerically investigate electrostatic instability of electron-positron pairs injected in an external electric field. The electric field is expected to be so strong that we cannot neglect effects of spatial variation in the 0th order distribution functions on the scale of the plasma oscillation. We assume that pairs are injected mono-energetically with 4-velocity u0>0 in a constant external electric field by which electrons (positrons) are accelerated (decelerated). By solving linear perturbations of the field and distribution functions of pairs, we find a new type of electrostatic instability. The properties of the instability are characterized by u0 and the ratio R of the braking time-scale (determined by the external electric field) to the time-scale of the plasma oscillation. The growth rate is as large as a few times the plasma frequency. We discuss the possibility that the excited waves prevent positrons from returning to the stellar surface.
Key words: instabilities - plasmas - stars: pulsars: general
Waves in homogeneous plasmas are well described by the linear perturbations that have the Fourier-harmonic dependence of the form . The properties of various wave modes have been extensively studied for various physical situations. Plasma instabilities, such as the two-stream instability, the Weibel instability (Weibel 1959), and many others have been recognized as important processes in astrophysics as well as in laboratory situations. For inhomogeneous plasmas, the Fourier-harmonic dependence is not assured in a strict sense. If the wavelength is short enough compared with the scale of the inhomogeneities, we may neglect effects of the spatial gradient of the distribution functions of particles, and carry out the Fourier-harmonic expansion. We call this treatment the local approximation hereafter.
When a static electric field exists in a plasma, it accelerates particles and leads to inhomogeneous velocity distributions of particles. Wave properties in an electric field were studied in geophysical researches, adopting the local approximation (e.g. Misra et al. 1979; Das & Singh 1982; Misra & Singh 1977). In many high-energy astronomical phenomena, electric field is expected to be so strong that the local approximation is not adequately applied. Such a situation typically appears in the magnetosphere of pulsars.
A spinning magnetized neutron star provides huge electric potential differences between different parts of its surface as a result of unipolar induction (Goldreich & Julian 1969). A part of the potential difference will be expended as an electric field along the magnetic field somewhere in the magnetosphere. Although a fully consistent model for the pulsar magnetosphere has yet to be constructed, several promising models have been considered. Among them, the polar cap model (Ruderman & Sutherland 1975; Sturrock 1971) assumes that an electric field parallel to the magnetic field lines exists just above the magnetic poles. The electric field accelerates charged particles up to TeV energies, and resultant curvature radiation from these particles produces copious electron-positron pairs through magnetic pair production. These pairs may provide gamma-ray emission by curvature radiation or synchrotron radiation as well as coherent radio emission and a source for the pulsar wind.
The localized potential drop is maintained by a pair of anode and cathode regions. In the cathode region the space charge density deviates from the Goldreich-Julian (GJ) density negatively, where is the angular velocity of the star and Bz is the magnetic field strength along the rotation axis. On the other hand, deviates positively for the anode. Outside the accelerator the electric field is screened out. In the polar cap model, especially for space charge limited flow model (Arons & Scharlemann 1979; Scharlemann et al. 1978; Fawley et al. 1977), where electrons can freely escape from the stellar surface, i.e., on the stellar surface, the formation mechanism of a static pair of anode and cathode, which can sustain enough potential drop for pair production, is a long-standing issue. Current flows steadily along the magnetic field line so that the charge density is determined by the magnitude of the current and field geometry with suitable boundary conditions. Good examples for space charge limited flow are in Shibata (1997). When and the electron density ( , where B is the magnetic field strength) is larger than the GJ number density ( , where -e is the electronic charge) on the stellar surface, a cathode is provided on the stellar surface. The cathode accelerates electrons. When the field lines curve away from the rotation axis, n deviates nagatively so much more for "away'' curvature, which enhances the cathode. Hence electrons continue to be accelerated, and potential drop becomes large enough to produce pairs. The mechanism of the electric field screening, i.e., a way to provide an anode, has been considered to be provided by pair polarization. Although most papers take it for granted that copious pair production can instantly screen the field, recently Shibata et al. (1998,2002) casted doubt on this issue; the electric field screening is not an easy task as considered usually.
Shibata et al. (1998,2002) investigated the screening of electric fields in the pair production region. They found that the thickness of the screening layer is restricted to be as small as the braking distance for which decelerating particles become non-relativistic, where is the electron mass. If the above condition does not hold, too many positrons are reflected back and destroy the negative charge region (cathode). In order to screen the electric field consistently, huge number of pairs should be injected within the small thickness . The required pair multiplication factor per one primary particle is enormously large and cannot be realized in the conventional pair creation models. Thus, some other ingredients are required for the electric field screening.
In the previous studies of the screening, pairs were assumed to accelerate or decelerate along the 0th order trajectories determined by . However, if an electrostatic (longitudinal) instability occurs, the excited waves may produce effective friction. Friction on particles change the charge polarization process. Thus, instability in the presence of an external electric field may have a relevance to this problem, which motivates us to make an exploratory study in this paper. Various instability mechanisms outside the accelerator have been studied for pair plasmas along with a primary beam in relation to coherent radio emission mechansims (Cheng & Ruderman 1977; Asséo et al. 1983; Hinata 1976; Lyubarskii 1992; Asséo et al. 1980). However, plasma instability inside the accelerator has not been studied. The Lorentz factor of the primary beam ( -107) is much larger than that of electron-positron pairs ( -103). In such a case, it is difficult to induce the two-stream instability. However, it is not clear whether pairs stably flow in the electric field or not. For typical pulsar parameters, the braking distance is cm, while the length-scale of plasma oscillation is cm (Shibata et al. 1998,2002), where . Particles are accelerated or decelerated in a period that is shorter than the typical time scale of plasma oscillation. Therefore, the distribution function is not uniform on the scale we consider. The local approximation is not adequate to deal with plasma oscillation.
Properties of a pair plasma in such a strong electric field have not been studied. Such studies may bring us a new key to understanding astrophysical phenomena. One of our purposes is to examine if electrostatic instability makes the screening easier. As the first step toward this purpose, in this paper we investigate electrostatic instability of pairs injected in an external electric field. Investigations of electrostatic waves, when we cannot adopt the local approximation, may be important not only for pulsars but also for other high-energy astronomical phenomena. Since an analytical treatment is difficult in this case, we simulate electrostatic waves numerically in idealized situations. In Sect. 2 we mention the two-stream instability without an electric field for comparison. In Sect. 3 we describe the situation we consider. The most simplified physical condition is adopted; the electric field, injection rate, injection energy are constant. In Sect. 4 we explain our numerical method. Our method can treat only linear waves. Numerical results are summarized in Sect. 5. Our results show a new type of plasma instability due to electric field. Section 6 is devoted to summary and discussion.
First of all, for reference, we consider one-dimensional (1-D)
homogeneous flows of electrons and positrons in the absence of
electric field. The pair-distributions are functions of the
For simplicity, we assume that the distribution functions
are expressed by the step function
|Figure 1: Distribution function of the toy model in Sect. 2.|
|Open with DEXTER|
When u2=u3, electrons and positrons distribute
continuously, though the average velocity is different.
In this case, the dispersion relation is reduced to a quadratic
equation, and we obtain
In a strong magnetic field as in the pulsar polar caps, transverse momenta of relativistic electrons and positrons are lost during a very short time via synchrotron radiation. These particles move along the magnetic field lines and their distribution functions are spatially 1-D. In this section we consider 1-D distribution functions of electron-positron pairs injected in an external electric field that is parallel to magnetic field lines. As the conventional theories for the two-stream instability implicitly assume, we neglect the toroidal magnetic field due to the global current of plasma. Only within this treatment the 1-D approximation is adequate. In order to simplify the situation, we assume the external electric field E0 is constant. In pulsar models, there exists a primary beam which produces current flow. The external electric field is determined by the complicated combination of the beam current, injected pair plasma, and GJ density. The constant E0 requires that the charge density of the beam, pairs, and GJ density cancels out, which may be an artificial situation. As will be discussed in Sect. 6, the approximation of constant E0is justified for a smaller rate of pair injection, which may be realized for actual pulsar parameters.
Anyway, we depart from actual pulsar physics, and
deal with plasma physics in an idealized situation hereafter.
We adopt the 1-D approximation and assume the existence of
the background charge which leads to constant E0.
In our treatment we totally neglects effects of the existence
of the background
on development of waves, and consider the behaviour of pair plasma only.
Pairs are assumed to be injected between
at a constant rate .
In our calculation the pair injection is monoenergetic with
Let us start from assuming the steady state of flows of
The distribution function f0(z,u)
satisfies the Boltzmann-Vlasov equation
Then, injected positrons (q=e>0) will
be decelerated as
Electrons (q=-e<0) will continue to be accelerated as
|Figure 2: Distributions of electrons and positrons in the phase space. Electrons (positrons) are in the region ( ). The electric field (schematically shown as the central corrugation in the figure) generates the disturbances F at the injection of pairs (u=u0), , and . The arrows schematically show the transports of F (see text).|
|Open with DEXTER|
In the region the pair distribution has separate two streams. Since this region is peculiar in our idealized model, we do not consider the waves in this region hereafter.
We consider linear perturbations of the distribution function
and electric field as
We directly solve the perturbations f1 and E1from the linearly perturbed Boltzmann-Vlasov equation
We set up grids along the 0th order trajectories of pairs
in the 2-D phase space
Following the Lagrangian method
we follow time evolution of f1 in these grids from Eq. (15).
Eq. (15) is rewritten as
On the other hand, the evolution of electric field is calculated by
the Eulerian method;
We have ascertained that results obtained from our numerical code satisfy the Gauss law. In addition we have checked our code by reproducing two stream instability in the absence of electric field, using the distribution functions in Sect. 2.
We have simulated electrostatic waves from various initial conditions and parameter values. We are interested in a parameter region R <1. In this region the typical wavelength of plasma oscillation is longer than the braking distance . We give an initial disturbance in a spatially limited region. As will be shown below, when , we find an absolute instability in which disturbance grows in amplitude but always embraces the original region, where the initial disturbances of F and are given. The condition means that the distance injected positrons move forward before they turn back, , is larger than the length-scale of the plasma oscillation . On the other hand, for , we find a convective instability in which disturbance grows while propagating away from the original region. The waves excited from the convective instability propagate backward. Empirically, the results do not largely depend on the spatial size of the pair injection region , which determines the minimum 4-velocity . In this section we show some examples of the instabilities found in our simulations. The parameters and initial conditions are summarized in Table 1.
Table 1: Parameters and initial conditions in computation.
The initial conditions are taken to satisfy the Gauss law.
Given the parameters R, u0, and ,
we set the initial values of the disturbances F and
The initial disturbances are confined within a small spatial region of one wavelength ( ). In the other region there is no disturbance. The initial perturbation of has a single sign with a form of a cosine curve. On the other hand, F- and F+ have the form of a sine curve, satisfying the Gauss law. The parameter induces asymmetry of the charge density of electrons and positrons. The total charge density ( ) does not depend on , and also has a form of a sine curve with . We have tried various values of in our simulations and find that the ratio of F- to F+ is settled as time passes irrespective of . When instabilities occur, growing wave modes end up dominating other modes. Therefore, results do not largely depend on the initial conditions.
First we describe the results for R=0.1 (the calculations RUN1-RUN3) and see the behaviour of the linear perturbations for various values of the parameter u0. In Figs. 3 and 4 we plot electric field for RUN1 for which R u0=1. In this calculation, the initial disturbance exists from to . Since positrons will turn the direction of motion after , we must follow the disturbance much longer than that. As is illustrated in Fig. 3, at the disturbance remains in the originally disturbed region. As time passes, the amplitude around the original region of the disturbance grows, and the wave packet spreads backward little by little. In the forward region , we do not observe any growing wave. Although particles move almost at the light velocity, the disturbances remain around the original region and the wave packet does not spread at the light velocity. In order to show the growth of the amplitude, in Fig. 5 we plot the time evolution of that is the electric field for the maximum amplitude. The maximum electic field oscillates over positive and negative regions. Initially changes complicatedly because of the initial conditions we artificially set. As time passes, smoothly grows while oscillating. The period of the oscillation of is . The growing time ti, where , is .
|Figure 3: Electrostatic waves for RUN1 for , 165, and 195.|
|Open with DEXTER|
Let us look into the the behaviour on a shorter time scale for RUN1. At a fixed position s, the local electric field grows while oscillating. However, when we see spatio-temporal behaviour, we notice that the spatial pattern propagates backward while changing their amplitude. As is shown in Fig. 4, waves propagate from backward with a growing amplitude. The amplitude becomes maximum around , and then the amplitude declines while propagating backward. This decline leads to the confinement of the wave packet. Even though waves pass the disturbed region many times, waves exist only in a spatially limited region. In the wave packet of there are multiple peaks and bottoms, and the most prominent one of them corresponds to . If we define the "phase velocity'' as the velocity of peaks (or bottoms), the phase velocity (-2.8 c) turns out to be faster than the velocity of light. As peaks propagate backward, the peak or bottom associated with alternates one after another, so that the position of hangs around the original region. When we define the group velocity by averaging the velocity of the position of for a longer time scale than the oscillation period, the group velocity turns out to be almost zero. The wavelength is about which is almost the same as the initial wavelength of the disturbance. Even if we start from another , the growing mode dominates others and the final wavelength is the same as this result, . The final wavelength of growing waves is unchanged for different initial conditions. Figure 6 shows charge density distributions at . The number density of electrons (positrons) has opposite (same) sign of the charge density in Fig. 6. The phases of number densities of electrons and positrons are the same. The amplitude of the positron density is always larger than that of the electron density. The difference in the number densities of positrons and electrons is proportional to the total charge density which satisfies the Gauss law.
|Figure 4: Electrostatic waves for RUN1 around .|
|Open with DEXTER|
|Figure 5: Time evolution of for RUN1.|
|Open with DEXTER|
|Figure 6: Charge-density distributions for in RUN1. The dashed line is the electric field at that time.|
|Open with DEXTER|
Next we discuss the results of RUN2 ( R u0=30>1). The initial distubance ranges from to . In Fig. 7 we show electric fields at several epochs. The wave profiles are not so simple compared to the case of RUN1 (R u0=1). There are waves propagating both forward and backward. This may be because the distance positrons move forward before they returns ( from their injection point) is longer than the typical wavelength . Though the properties of the waves are complex, we can see that there is an instability in this case, too. The waves around the original region grow while diffusing both forward and backward. We note that a separate component of the disturbance appears around (s<300 in this case). This disturbance may be due to two stream instability, which grows faster than in the other region.
|Figure 7: Electrostatic waves for RUN2. The thick lines are for and 90. The thin lines are for and 340.|
|Open with DEXTER|
Although we do not show here, for a much larger value of R u0, absolute instabilities are found to occur in our simulation. However, we do not further pursue this issue because we need to calculate over a much wider region of s and resultant memory in computation becomes large for a large value of u0.
In cases for R u0=0.3<1 (the calculations RUN3), absolute instability does not occur, but convective instability propagating backward occurs (see Fig. 8). In RUN3 the initial disturbance is given from s=420 to with a wavelength . The disturbance of the electric field is seen to propagate backward, and before long characteristic waves grow. The wavelength of the growing wave ( ) is slightly longer than the initial length and almost constant. As the disturbance propagates, the wave packet spreads and the number of waves in the packet increases. The growing time of is about (see Fig. 9). The phase velocity, , is faster than the velocity of light. The amplitude of each peak in the packet initially grows. As the peak approaches the head of the wave packet, the growth of the amplitude turns over to damping. This behaviour is similar to that in the absolute instability for RUN1. We obtain the group velocity . Figure 10 shows that the charge density is dominated by positrons. We have simulated for R u0=0.3 with various initial conditions. However, in any case there is no sign of wave instability propagating forward.
|Figure 8: Backward waves for RUN3.|
|Open with DEXTER|
|Figure 9: Time evolution of for RUN3.|
|Open with DEXTER|
|Figure 10: Charge-density distributions for in RUN3. The dashed line is the electric field at that time.|
|Open with DEXTER|
We have investigated for R=0.01 also (the calculations RUN4-RUN6), and confirmed that the qualitative results are determined by the value of R u0 (see Table 2). The absolute instability and convective instability are induced for R u0=1 and R u0=0.3, respectively. For R u0=0.03 (RUN6), however, we do not find any instability. A smaller value of R means lower density of pair plasma for a given value of the external electric field. In such a low density plasma, particles tend to be less affected by forces from other particles compared to the external electric field. Therefore, too small value of makes the plasma stable. From the above results, we may conclude that electrostatic instability is mainly determined by the parameter R u0. The wavelength and the growing time are a few or ten times the typical scales of the plasma oscillation, and , respectively.
Table 2: Rough values of the growing time, wavelength, and phase velocity. The characters "A'' and "C'' represent the absolute instability and convective instability, respectively.
The value R has been assumed to be smaller than 1 heretofore. We have also simulated for and found instabilities. The wave properties are as complicated as the example in RUN2, so that we do not report the details of the results in this paper. When , E0 is so small that . In the limit of E0=0 ( ), the instability does not occur as was shown in Sect. 2. However, we have not tried the cases , because of poor computational capacity so far.
Judging from the backward spread of wave packets and the negative phase velocity, it is seen that returning positrons play a decisive role in the instabilities. The dominance of positrons in the charge density in comparison with electrons also suggests that the instabilities are due to returning positrons. It is remarkable that positrons pass the same region twice, forward and backward. The typical wavelength of plasma oscillation ( ) may resonate with the distance positrons move forward, . As we have mentioned in Sect. 4, the excited electric field generates the disturbances of pairs, F(u0), at their injection. The value of F(u) is transported along the trajectories of pairs, conserving their value. Since F+ and F- acquire the same value at their injection, the contribution to the charge density is almost canceled out as long as electrons and positrons move forward together at . When positrons turn around, the charge is polarized. As for positrons, F+(u,s) at conserves the information on the electric field in the past ( ), where . The displacement between s and s' conforms to the trajectories of positrons. Therefore, the charge density of positrons (except for ) is proportional to a superposition of displaced at different times. This superposition will increase the amplitude of the charge density in response to evolution of .
Since is constant along the characteristics, the value F+ remains to be finite even after positrons enter the region where , while changes as long as exists. Our numerical results imply that cancels out the charge density due to in the region where . This process prevents the waves from spreading backward at the speed of light. The propagation of the disturbance by electrons is relatively simple. The contribution due to prevents the waves from spreading forward, as does. In RUN3 and RUN5 (R u0 <1), positrons turn back and move backward quickly before cancels out the charge density, so that growing waves propagate backward.
In this paper we have found a new type of instability in electron-positron pair flows injected in an external electric field, which is assumed to be spatially constant. The properties of the instability are characterized by the ratio R (the braking time-scale to the typical oscillation time-scale of the plasma ) and 4-velocity u0 at injection. For absolute instability is induced, while convective instability propagating backward is excited for . The growing time in amplitude is as short as a few times the time-scale . The wavelength is also several times . The instabilities are caused by returning positrons. For , the pair plasma turns out to be stable. A small value of R implies that the plasma density is so low compared to the electric field E0. For the collective interaction of the pair plasma is not important, so that each particle moves along the trajectory determined by E0 independently of other particles.
Growing electrostatic waves may work as frictional forces.
In this paper we have treated waves as linear perturbations,
following the propagation of disturbances in the distribution function
and electric field.
Our method does not allow us to follow processes
of gaining or losing kinetic energy of each particle from the waves.
The quasi-linear theory is not applied
to deal with the reaction of particles as it is,
because the disturbances do not have the Fourier-harmonic dependence.
Thus, we consider the qualitative character of the effective reaction force
from a numerical treatment as follows.
The spatial averages of F and
oscillate with time
in our simulations.
Therefore, the expectation values of F and can be considered to be zero.
On the other hand, when waves grow or are attenuated,
the spatial average of the cross term,
may have a finite value.
As is the case with the quasi-linear theory,
the 0th distribution function may change,
following the 2-nd order order approximations of the Boltzmann-Vlasov equation:
We plot G(u) in Fig. 11 for in RUN1. The modulation pattern of G(u) does not change, but the amplitude grows with time. Apparently, the modulation pattern is asymmetric for particles of u>u0=10 (electrons) and u< u0 (positrons). These qualitative behaviour is common for the other RUNs. As Fig. 2 and Eqs. (17) and (18) show, perturbations are generated at u=u0and (or ). Figure 11 shows G(u) for a region around u=u0 only, and outside of this region G(u) has also significant value due to the disturbances generated at and . However, the modulation pattern of G(u) for such regions oscillates with time.
|Figure 11: The "average force'' G(u) for in RUN1.|
|Open with DEXTER|
In the usual two-stream instability, the excited waves accelerate background fluid, and decelerate beam fluid. As is shown in Fig. 11, the direction of the reaction force depends on u even in the same species of particles. The amplitude of G(u) takes always the maximum value at u=u0. If the effective reaction force grows enough, positrons (electrons) just injected ( ) feel positive (negative) force as is shown in Fig. 11. Thus, the reaction force may make particles tend to stay around the regions of u=u0. The integral of G(u) around u=u0 (roughly from u=-100 to 10 for Fig. 11) for positrons is also positive. Such positrons are accelerated by the reaction force on average. However, the integral becomes negative all the time, if we include the contribution due to the disturbances generated at , though G(u) for a large |u| oscillates with times. The returning positrons injected at feel a negative reaction force on average. On the other hand, the absolute value of the integral for electrons is much smaller than those for positrons. Therefore, the reaction force does not work as usual "frictional'' force between electrons and positrons. The particles just injected are most affected by the reaction force, and lose (or gain) their energy owing to waves. If the reaction force is strong enough, positrons just injected may have difficulty to turn back. If positrons suffer from such frictional force, the distribution function f0 should be largely altered. This may help to solve the problem of the electric field screening in the pulsar polar caps.
If excited waves grow enough to change trajectories of particles, the waves cannot be treated as the linear perturbations. In a strong electric field, |f1| can be as large as f0 before achieves. We may simplify the energy-loss process of particles due to perturbed electric field as follows; particles from u=u0 to lose their energy owing to the waves at the same rate, where is the equivalent width of particles interacting the waves. Then, the growth rate of the field energy is roughly considered as the average energy loss rate of pairs, which means , where is the temporal change of due to the reaction force. Here, we assume that the energy density of the waves attributed to induced particle motions is comparable or negligible to the energy density of the electric field E1. When , the frictional force is sufficient to alter trajectories of pairs. The above condition is rewritten as . Assuming , , and in esu (these values may be typical for the pulsar polar cap), |E1| is needed to be larger than in esu to change the distribution of pairs. Even if |E1| < |E0|, the excited disturbances may change the 0th distribution of pairs. However, we need numerical simulations, which can deal with non-linear process, in order to check how the reaction force modifies the distribution function f0.
As a first step to deal with behaviour of pairs in an electric field, in this paper we have assumed that the background charge distribution cancels out the modification of E0 due to injected electron-positron pairs. Of course, this simplification may not be appropriate for pulsars, while it makes computation easier. Inhomogeneous electric field might play an important role in plasma instability. Let us check whether our model can be used when the background charge density is constant for s>s0. The charge density changes for s>s0 owing to the pair injection and electric field should be modified. The charge density decreases with distance as n0 s for s>s0. In this case the variation of over the moving distance of injected positrons before they turn back is . Therefore, in the case of stable plasma ( ) the constant electric field is a good approximation even for the constant background. Although our simulations show a new possibility of plasma instability around pulsars, we need to simulate with an inhomogeneous electric field for in order to conclude whether instability occurs in actual situations on pulsars. In any case, the condition is an necessary (but not sufficient so far) condition to induce the instability.
Let us consider implications for the pulsar polar cap.
We suppose that the primary electron-beam is accelerated from z=0,
and its Lorentz factor becomes
at the pair production front (PPF) (z=L).
In this case the average electric field is
The braking time is expressed as
By curvature radiation an electron of
photons per unit time,
where r is the radius of curvature of the field line.
We express the number density of the primary beam as
Assuming curvature gamma-rays immediately turn into pairs,
the pair-injection rate is approximated as
Adopting the average electric field, the ratio R becomes
At present it is not clear that the electrostatic instability we have considered in this paper is an important process for the screening of electric field above the pulsar polar cap. However, there may be extreme environments (magnetars etc.) where R u0 is as large as one. We expect that studies on the plasma instability in electric fields lead to opening a new approach to high-energy astrophysics.
This work is supported in part by a Grant-in-Aid for Scientific Research from Ministry of Education and Science (No.13440061, F.T.). One of the authors (K.A.) is supported by the Japan Society for the Promotion of Science.