Physics Letters A 363 (2007) 225–231
www.elsevier.com/locate/pla

Self-guiding electromagnetic beams in relativistic electron–positron plasmas
T. Tatsuno a,b , M. Ohhashi a , V.I. Berezhiani a,c,∗ , S.V. Mikeladze ca Graduate School of Frontier Sciences, The University of Tokyo, Kashiwa, Chiba 277-8561, Japanb Center for Scientific Computation and Mathematical Modeling, The University of Maryland, College Park, MD 20742-3289, USAc Institute of Physics, 6 Tamarashvili, Tbilisi 0177, Georgia

Received 15 September 2006; received in revised form 27 October 2006; accepted 30 October 2006
Available online 7 November 2006
Communicated by V.M. Agranovich

Abstract
Nonlinear interaction of an intense electromagnetic (EM) beam with relativistically hot electron–positron plasma is investigated by invoking thevariational principle and numerical simulation, resting on the model of generalized nonlinear Schrödinger equation with saturating nonlinearity.
The present analysis shows the dynamical properties including the possibilities of trapping and wave-breaking of EM beams. These properties of
EM beams may give a significant clue for the gamma-ray burst.
© 2006 Elsevier B.V. All rights reserved.
PACS: 52.27.Ep; 52.27.Ny; 52.35.Mw; 98.70.Rz
Keywords: Electron–positron plasma; Special relativity; Gamma-ray burst; Nonlinear Schrödinger equation; Saturating nonlinearity; Variational principle

1. Introduction
The problem of electromagnetic (EM) wave propagation andrelated phenomena in relativistic plasma have attracted considerable attention in the recent past. Relativistic electron–positron(e–p) dominated plasmas are created in a variety of astrophysical situations. Electron–positron pair production cascades are
believed to occur in pulsar magnetospheres [1]. The e–p plasmas are also likely to be found in the bipolar outflows (Jets) inactive galactic nuclei (AGN) [2], and at the center of our own
Galaxy [3]. In AGNs, the observations of superluminal motions are commonly attributed to the expansion of relativistic
e–p beams in a pervading subrelativistic medium. This model
implies copious pair production via γ –γ interactions creating
an e–p atmosphere around the source. The actual production of
e–p pairs due to photon–photon interactions occurs in the coronas of AGN accretion disks, which upscatter the soft photons
emitted by the accretion disks by inverse Compton scattering.

* Corresponding author.

E-mail address: vazhab@yahoo.com (V.I. Berezhiani).
0375-9601/$ – see front matter © 2006 Elsevier B.V. All rights reserved.
doi:10.1016/j.physleta.2006.10.096

The presence of e–p plasma is also argued in the MeV epochof the early Universe [4]. On the other hand the contemporary progress in the development of super-strong laser sources
with intensities I = 1021–23 W/cm2 has also made it possible
to create relativistic e–p plasmas in the laboratory by a variety
of experimental techniques [5]. Elucidation of the electromagnetic wave dynamics in a relativistic e–p plasmas will, perhaps,
be an essential determinant of the radiation properties of astrophysical objects as well as of the medium exposed to the fieldof super-strong laser radiation.
Wave self-modulation and soliton-formation is, perhaps, oneof the most interesting and significant features of the overall
plasma dynamics. The existence of stable localized envelope
solitons of EM radiation has been suggested as a potential
mechanism for the production of micro-pulses in AGN and pulsars [6–8]. In the early Universe localized solitons are strong
candidates to explain the observed inhomogeneities of the visible Universe [9,10], and vortex soliton can be considered toplay a significant role for the formation of the observed structure of the Universe [11]. The gamma-ray bursts (GRBs) andtheir afterglows are likely to be the result of energy dissipation from a relativistically expanding outflow [12]. A pointing-226

T. Tatsuno et al. / Physics Letters A 363 (2007) 225–231

flux driven outflow from a magnetized rotator is a promisingparadigm for GRB engines and there have been various implementations of this concept [13–17]. Recently Lyutikov andBlackman suggested that gamma-rays are emitted at the point
where MHD breaks down due to the overturn instability of large
amplitude electromagnetic wave [18].
In the present Letter we explore the mechanism for the localization of multi-dimensional intense EM radiation in pure e–p
plasmas. Assuming the plasma to be transparent to the beam,
and applying a fully relativistic hydrodynamical model, we
demonstrate the possibility of beam self-trapping leading to theformation of stable 2D solitonic structures. The high-frequencypressure force of the EM field (tending to completely expel thepairs radially from the region of localization) is overwhelmedby the thermal pressure force which opposes the radial expansion of the plasma creating conditions for the formation of thestationary self-guiding regime of beam propagation.

In this section we apply our general formulation to the problem of self-trapping of EM beams in pure e–p plasmas withrelativistic temperatures. We assume that the equilibrium stateof the plasma is characterized by an overall charge neutrality+−+n−∞ = n∞ ≡ n∞ , where n∞ and n∞ are the unperturbed number densities of the electrons and positrons in the far region of
the EM beam localization. In most mechanisms for creating e–p
plasmas, the pairs appear simultaneously and due to the symme− =T+ ≡T ,try of the problem it is natural to assume that T∞∞∞−+
where T∞ and T∞ are the respective equilibrium temperatures.
We shall assume that for the radiation field of interest,
the plasma is underdense and transparent, i.e.,  = ωe /ω  1,
where ω is the mean frequency of EM radiation and ωe =
(4πe2 n∞ /m0e )1/2 is the plasma frequency. Following Ref. [19]we introduce a temperature dependent momentum ± =G± p± and relativistic factor Γ ± = G± γ ± , where G = K3 (z)/K2 (z) is the temperature dependent factor with Kn the nth ordermodified Besselfunction of the second kind (z = m0e c2 /T ),p± and γ ± = 1 + (p± /me c)2 are respectively momentum
and relativistic factor of e–p particles. Introducing the dimensionless quantities t˜ = ωt, r̃ = (ω/c)r, T̃ ± = T ± /m0e c2 ,˜ ± = ± /(m0e c), andÃ = eA/(m0e c2 ), φ̃ = eφ/m0e c2 , ñ± = n± /n∞ , we arrive at the dimensionless equations [19],



∂
∂ 2A
− A + ∇φ −  2 J+ − J− = 0,
2
∂t
∂t
where
  K2 (1/T ± )T ±  f T± =
exp G± T ± /T ± ,G± (T ± )Q = Q0 (ξ, x1 , y1 , z2 ) + Q1 (ξ, x1 , y1 , z2 ),(1)±0⊥ = ∓A0⊥ .

(5)

(9)

To the next order (in ), the transverse component of Eq. (1)reads:
−a

∂±
∂A1
1
+ ∇⊥ Γ0± = ±a
∓ ∇⊥ φ0 .
∂ξ
∂ξ

(10)

Averaging Eq. (10) over the fast variable ξ we obtain ∇⊥ Γ0± =∓∇⊥ φ0 . Using these relations, we obtain ∇⊥ φ0 = 0 andΓ0 ≡ Γ0± = const.

(11)

Note that from the lowest order of Eqs. (2) and (4), we conclude−+−
that n+0 = n0 ≡ n0 and T0 = T0 ≡ T0 . The relation between
EM field and temperature can be found by invoking Eq. (11).
Using Eqs. (8) and (9) and determining integration constantsuch that A → 0 and T0 → T∞ , we obtain
G2 (T0 ) = G2 (T∞ ) − |A|2 .

(12)

The equation for slowly varying envelope A of EM beam canbe obtained from Eq. (5). To the lowest order we find
a1

(4)

(8)

where A is the slowly varying envelope of the EM beam, x̂ andŷ denote unit vectors, and c.c. is the complex conjugate.
We now give a short summary of the steps in the standard
multiple-scale methodology (Ref. [20]). In the lowest orderin , the transverse (to the direction of EM wave propagation z)
component of Eq. (1) reduces to

(2)
(3)

(7)

where (x1 , y1 , z2 ) = (x, y,  2 z) denote the directions of slow
change, and a1 = (a 2 − 1)/ 2 ∼ 1. We further assume that thehigh-frequency EM field is circularly polarized,
1
A0⊥ = (x̂ + i ŷ)A exp(iξ/a) + c.c.,2

2. Basic equations∂±∂A+ ∇Γ ± = ∓
∓ ∇φ,
∂t
∂t
n±
= const,±Γ f (T ± )∂n±+ ∇ · J± = 0,∂tφ =  2 n− − n+ ,with J± = n± ± /Γ ± and Γ ± = (G± )2 + (± )2 . Here A
and φ are the vector and scalar potentials (∇ · A = 0). The tilde
is suppressed for convenience.
Of various techniques that could be invoked to investigate
Eqs. (1)–(6) to study the self-trapping of high-frequency EM
radiation propagating along the z-axis, we choose the method
presented in the excellent paper by Sun et al. [20]. The method
is based on the multiple scale expansion of the equations in thesmall parameter . Assuming that all variations are slow compared to the variation in ξ = z − at, we expand all quantities
Q = (A, φ, ± , n± , . . .) as∂ 2 A0⊥∂ 2 A0⊥n0 (T0 )2− ∇⊥A0⊥ − 2+2
A0⊥ = 0.
2
∂ξ ∂z2
G(T∞ )∂ξIn deriving this equation, we have used the relationΓ0 = G2 (T0 ) + |A|2 = G(T∞ ).

(13)

(14)

While from Eq. (2), we have(6)n0 (T0 ) =f (T0 ).
f (T∞ )(15)

T. Tatsuno et al. / Physics Letters A 363 (2007) 225–231

227

Substituting Eq. (8) into Eq. (13) we find

∂A
2 
2
(16)
A+
+ ∇⊥
1 − n(T0 ) A = 0,∂z
G(T∞ )
where subscripts for variables (x1 , y1 , z2 ) are dropped for simplicity. We also assumed without loss of generality that (a 2 −1)/ 2 a 2 = 2/G(T∞ ), which in dimensional units (provided
that a = ω/kc) coincides with the linear dispersion relation of
the EM wave in an e–p plasma:
2i

2ωe2
.
(17)
G(T∞ )
Thus, the dynamics of EM beams in hot relativistic e–p plasma
has become accessible within the context of a generalized nonlinear Schrödinger equation (16).
We first seek the localized 2D soliton solutions of Eq. (16),and analyze the stability of such solutions. Making the selfevident
renormalization of variables z → zG(T∞ ), r⊥ →√r⊥ G(T∞ )/2, Eq. (16) can be written asω2 = k 2 c2 +

∂A
2
(18)
A + Ψ A = 0,
+ ∇⊥∂zwhere Ψ = 1 − n0 (T0 ) represents the generalized nonlinearity.
The companion equation (15) can be viewed as a transcendental
algebraic relation between T0 and |A|2 , i.e. Ψ is an implicit
function of |A|2 [Ψ = Ψ (|A|2 )]. We note that Eq. (18) can be
written in the Hamiltonian form i∂z A = δH /δA∗ , where the
Hamiltonian isH = dr⊥ |∇⊥ A|2 − F |A|2 ,
(19)
i

t

and F (t) = 0 Ψ (t ) dt . This implies that Eq. (18) conserves
the Hamiltonian H in addition to the power (photon number)
N = dr⊥ |A|2 .
To find a stationary, nondiffracting axisymmetric solution, we use the representation of vector potential A/Ac =U (r) exp(iλz), where Ac = [G(T∞ )2 − 1]1/2 , r = (x 2 + y 2 )1/2
and λ is the nonlinear wave-vector shift. The radially dependentenvelope U (r) obeys a nonlinear ordinary differential equation:
 
1 dU
d 2U
+
− λU + Ψ U 2 U = 0.
(20)
2
r dr
dr
The profiles of the field U (r), the plasma density n0 (r), andthe temperature T0 (r) of the stationary solution can be foundin Fig. 3 of Ref. [19] for λ = 0.1; The plasma temperature anddensity is reduced in the region of field localization. Similarplots could be obtained for all allowed values of λ. When λ →λc , where λc is the upper bound of the propagation constant,the plasma cavitation takes place, i.e. the plasma density andtemperature tends down to zero at r = 0. Appearance of zerotemperature is not surprising since the corresponding region is
the “plasma vacuum”; all particles are gone away.
The stability of obtained soliton solutions can be investigated using the stability criterion of Vakhitov and Kolokolov[21]. According to this criterion the soliton is stable against
small arbitrary perturbations if
dN
> 0,
dλ

(21)

Fig. 1. The beam power N versus λ (T∞ = 1).

where N is the photon number or more precisely the power of
the trapped mode. In our case, the dependence of N on λ isshown in Fig. 1. One can see that dN/dλ > 0 everywhere andconsequently the corresponding solution is stable for 0 < λ <λc .
3. Nonlinear dynamics based on variational approach
The complex dynamics of a beam governed by Eq. (18) canbe analyzed by the variational approach [22]. This approach
determines the relations between the characteristic parametersof the localized solution approximated by a trial function. The
variational method gives qualitatively good results, provided
the beam does not undergo structural changes during its evolution. The first standard step is to construct the Lagrangian
L=

∂A
1
A∗
− c.c. − H,
2
∂z

(22)

where H is the Hamiltonian density [H = dr⊥ H, see
Eq. (19)]. In the optimization procedure, the first variation of
the variational function must vanish on a suitably chosen trial
function. As a trial function, we will use the Gaussian-shapedbeam,
A = Λ(z) exp −r2+ ir 2 b(z) + iφ(z) ,2a(z)2(23)

with the amplitude Λ, the beam radius a, the wave front curvature b and the phase φ as the unknown functions of the propagation coordinate z respectively, which will be furthermore usedin order to make the variational functional an extremum. Substituting expression (23) into Eq. (22) and demanding that thevariation of the spatially averaged Lagrangian with respect toeach of these parameters is zero, we obtain the corresponding
set of Euler–Lagrange equations,d  2 2
Λ a = 0,dz

(24)

 2  K(Λ2 )
4
2
d 2a
Λ −
=
−
K
,
dz2
a3 a
Λ2

(25)

b=

1 da
,
4a dz

(26)

228

T. Tatsuno et al. / Physics Letters A 363 (2007) 225–231

to be solved for the three functions Λ, a, and b, where the function K(u) is defined as∞K(u) = 42
dp pF ue−p .

(27)

0

Eq. (24) is nothing but a statement of the fact that during the
EM beam evolution its power is conserved,
N = πΛ2 a 2 = πΛ20 a02 ,(28)

where Λ0 and a0 are respectively the initial amplitude and theinitial “radius” of the EM beam at z = 0.
Using Eq. (24), the integration of Eq. (25) leads to
1 da
2 dz

2

+ V (a) = H = V (a0 ),(29)

Λ20 a022a 2
2−
Ka 2 a02 Λ20
a2

(30)

Fig. 2. Potential profile of V (a) versus a in the case of N = 8.79 and N = 9.74.

where
V (a) =

plays the role of an effective potential for the evolution of the radius a. We have assumed that the initial beam have a plane front(or zero curvature) [da/dz|z=0 = 0 = b(0)]. Using the analogy
with a particle in a potential well, we can acquire a deeper physical understanding of light beam dynamics. Choosing the initial
beam radius a0 to be equal to the equilibrium radius ae , a stationary solution of Eq. (29) is obtained if ∂V /∂a|a=ae = 0. Notethat −∂V /∂a is equal to the right-hand side of Eq. (25). The
equilibrium radius of the beam is readily found to be  K(Λ20 ) −1.
ae2 = 2 K Λ20 −
Λ20

(31)

In the subsequent analysis we will consider the small temperature case T∞  1, and the case of moderately high temperature T∞ ∼ 1. We do not consider here the ultrarelativistic
1) because applied model equations
temperature case (T∞fail to adequately describe the plasma dynamics due to the neglect of heavy particle production.
In the small temperature case the nonlinear term in Eq. (18)
reduces to a simple analytic expression


3/2

,
Ψ |A|2 = 1 − 1 − |A|2
(32)
which gives the function F (u) as2
2
F (u) = u + (1 − u)5/2 − .
(33)
55
We will apply the established general formalism to the saturating nonlinearity given by Eq. (32). The function K becomes[see Eq. (27)]
46 22u 2u24(1 − u)1/2
−
+
5
15
15
5
8
4
4 46
− arctan(1 − u)1/2 − log u −− log 4 .
5
5
5 15
(34)

K(u) = 2u +Fig. 3. The parameter regime where solitary wave can be trapped for T∞  1.

Note that normalized strength of the field A (= A/Ac ) is
restricted from above |A|  1. Above this value the wavebreaking of the field takes place. Substituting Eq. (34) into
Eq. (25), we can investigate the nonlinear dynamics of the beamfor different initial conditions.
Using Eq. (34) we can find the effective potential V (a) [seeEq. (30)]. The shape of effective potentials in the case of N =8.79 and N = 9.74 are illustrated in Fig. 2. From the shape of
the potential we conclude that the “effective particle” (i.e. thebeam) can be trapped in the potential well. The endpoints ofthese lines (for small a) correspond to the occurrence of zerodensity. As a decreases the beam amplitude Λ increases andapproaches the unity at endpoints. Note that the potential at theendpoint is negative for N = 9.74. Thus, the cavitation can takeplace when the initial radius a0 is large enough. It is obviousthat a(z) is an oscillating function provided that a0 is in thetrapping region.
The beam never diffracts if H < 0 which coincides withthe general criterion established by Zakharov et al. [23]. Fig. 3
shows the parameter regime where the solitary wave could betrapped. The central line denotes the equilibrium curve, andthe shaded region between two dashed curves denotes a region
where EM beam is trapped. The lower dashed line corresponds

T. Tatsuno et al. / Physics Letters A 363 (2007) 225–231

to the zero Hamiltonian. Thus, below this line, the electromagnetic beam with any parameters will diffract. Above the upperdashed line, the electromagnetic beam amplitude will grow up
to wave-breaking limit. Therefore, the beam will be trapped inoscillatory regime provided that its parameters are situated inshaded region in Fig. 3.
4. Numerical simulation
In the preceding analysis, we applied a variational approach
involving a Gaussian trial function. The main limitation of thisapproach is that it is valid only in the aberration-less approximation, i.e. the approach is unable to account for structuralchanges in the beam shape. Such aspects of the beam dynamics are better delineated by numerical simulations. The detailed
dynamics of arbitrary field distribution must be studied by direct simulations of Eq. (18). The guidelines for simulation are
still provided by approximative analytical approaches.
The initial profile of the beam is taken to be Gaussian A =Λ0 exp(−r 2 /2a02 ). In the case of small temperature T∞ = 0.01
the initial parameters are taken in the trapping region shown in
Fig. 3. The initial focusing or defocussing before reaching theequilibrium, can be seen in Fig. 4, where the z-propagations of
the beam fields |A(r = 0, z)| with same initial power N = 9.5
but with different amplitudes Λ0 = 0.37 and Λ0 = 0.95 are simultaneously drawn. Corresponding initial states of the beams
are respectively situated on the left and right side of equilibrium curve in Fig. 3. Fig. 4 shows that the beams parameters are oscillating around the equilibrium with initial focusing(Λ0 = 0.37) or defocussing (Λ0 = 0.95) in agreement with theprediction of variational approach.
Fig. 5 shows the field distribution versus radius r and the
propagating coordinate z for Λ0 = 0.37. Note, however, that
due to the appearance of the radiation spectrum the amplitudesof the field oscillations are monotonically decreasing with increasing z. For larger z the formation of a ground solitonic state
may take place due to the damp-out of the oscillations. If the initial profile of the beam is close to the equilibrium one, then the
beam quickly reaches the profile of ground-state equilibrium,

Fig. 4. The beams parameters are oscillating around the equilibrium with initial
focusing (Λ0 = 0.37) or defocussing (Λ0 = 0.95).

229

and propagates for a long distance without much distortion ofits shape. Since variational approach is unable to account thestructural changes in the beam shape and corresponding formation of radiation spectra we can expect that the trapping regionin (N, Λ0 ) plane obtained by numerical simulations will bedifferent from the one shown in Fig. 3. The result of these simulations is presented in Fig. 6. The curve ‘e’ corresponds tothe equilibrium state, the curve ‘h’ the zero Hamiltonian, andthe curve ‘t’ the trapping boundary, respectively. The Gaussian
beam with initial parameters (i.e. N and Λ0 ) in the region below line h has a positive Hamiltonian and will be diffracted. The
trapping region of the beam is the area between lines h and t.
The beam with the parameters in this region will be trapped inself-guiding regime of propagation and will either focus or defocus to the ground state, exhibiting damped oscillations around
it. The initial focusing (defocussing) takes place if the parameters are in the region between lines t and e (between e and h).
Note that the area of the trapping region shown in Fig. 6 is larger
than it follows from variational approach (see Fig. 3). This enlargement of the trapping area is related to the radiation losses
during the beam convergence to the equilibrium state.
The beam with parameters in region above trapping line t
will focus until wave-breaking and plasma cavitation takes

Fig. 5. The field distribution versus r and z for Λ0 = 0.37.

Fig. 6. The trapping region in (N, Λ0 ) plane obtained by numerical simulations
for T∞  1. The curve ‘e’ corresponds to the equilibrium state, the curve ‘h’
the zero Hamiltonian, and the curve ‘t’ the trapping boundary, respectively.

230

T. Tatsuno et al. / Physics Letters A 363 (2007) 225–231

Fig. 7. The field |A(r = 0)| versus z for the cases when the beams parameters
are in below (the curve 1) and above (the curve 2) the line t.

Fig. 9. The plasma temperature T (r) versus z during the trapped beam evolution
toward the stable equilibrium.

5. SummaryFig. 8. The trapping region in (N, Λ0 ) plane for T∞ = 0.3 (the notations aresame as in Fig. 6).

place. Fig. 7 shows the plot |A(r = 0)| versus z for the cases
when the beams parameters are in the regions below (thecurve 1) and above (the curve 2) the line t. One can see thatin the first case the beam trapping takes place while, in the second case the beam amplitude increases up to the wave-breaking
limit (A = 1).
For relativistic high temperature case the EM beam dynamics is similar to what we observed in the low temperature case.
Fig. 8 shows the trapping region for T∞ = 0.3 (the notations aresame as in Fig. 6). Quantitative difference stems from the factthat in the high temperature case the mass of the e–p pairs are
modified by the temperature dependent G-factor. It is interesting to mention that in this case in the region of field localization
the plasma temperature can be decreased considerably. In Fig. 9
we demonstrate T (r) versus z during the trapped beam evolution toward the stable equilibrium.
It is thus confirmed by numerical simulations that the beam
trapping is robust especially for low background temperature
regime. The beam can achieve the equilibrium state from a wide
range of parameters, even far from the equilibrium, as it is predicted by our analytic approach. However, the area of trapping
region reduces as temperature increases.

We have investigated the nonlinear propagation of strong 2D
EM radiation in a relativistic, unmagnetized electron–positronplasma. The treatment is fully relativistic—in the thermal motion as well as in the coherent motion of the plasma particles.
The fact that relativistically hot e–p plasmas are capable of sustaining high amplitude localized structures of high amplitude
electromagnetic fields should be important to understand thecomplex radiative properties of different astrophysical objects
where such plasmas are considered to exist.
By applying a variational technique with a Gaussian trial
function, we have demonstrated the possibility of different
regimes of EM beam propagation. In particular, the trapped
beam exhibits an oscillatory behavior around the stationary solution. However, since the system contains an upper bound of
the field amplitude, the parameter regime is restricted for theseoscillations. If the initial beam intensity is large enough, theponderomotive force acts strongly leading to the beam overfocusing, and eventually the EM field gives rise to the plasma
cavitation and wave-breaking.
The region in parameter space where the beam can be regularly trapped and oscillate is found. We have also performed thenumerical simulations and demonstrated the parameter regime
where beam trapping takes place. Comparing the results of variational approach and the results of numerical simulations wehave confirmed that analytical estimate gives a qualitatively
good prediction for nonlinear wave dynamics. Due to the radiation and self-deformation of the beam, however, it is shownthat the beam may be more robust than it follows from variational method.
In our consideration effects related to the temporal reshapingof the radiation have been ignored. However, we can easily generalize our results by allowing a temporal variation of the field.

T. Tatsuno et al. / Physics Letters A 363 (2007) 225–231

If we assume that A weakly depends on τ = t − z/a, Eq. (18)acquires an additional term related to the wave group velocity
dispersion [∼ (ωe2 /ω2 )∂ 2 A/∂τ 2 ]. In transparent plasma casethis term can affect the long time dynamics of self-guiding
channel. In particular due to weak modulation instability [24]the self-trapped beam eventually will break into a train of spatiotemporal solitons, i.e. the light bullets. Since general dynamical properties of the light bullet formation should not besensitive to the type of saturating nonlinearity [25,26], we expect similar behavior of the EM beams in our case as well.
We may apply the present results to verify variety of GRB
models. The obtained results could be useful to understand theradiative properties of astrophysical objects like the AGN, thepulsars and extragalactic electron–positron Jets.
Acknowledgements
The work of V.I.B. and S.V.M. was partially supported by
ISTC and Georgian NSF.
References[1] P.A. Sturrock, Astrophys. J. 164 (1971) 529;
M.A. Ruderman, P.G. Sutherland, Astrophys. J. 196 (1995) 51;
F.C. Michel, Theory of Neutron Star Magnetospheres, University of
Chicago Press, Chicago, 1991.
[2] M.C. Begelman, R.D. Blandford, M.J. Rees, Rev. Mod. Phys. 56 (1984)
255.
[3] M.L. Burns, in: M.L. Burns, A.K. Harding, R. Ramaty (Eds.), Positron–
Electron Pairs in Astrophysics, American Institute of Physics, New York,1983.
[4] S. Weinberg, Gravitation and Cosmology, Wiley, New York, 1972;

231

P.J.E. Peebles, Principles of Physical Cosmology, Princeton Univ. Press,
Princeton, 1993.
[5] M. Perry, G. Mourou, Science 264 (1994) 917;
S.P. Hatchett, et al., Phys. Plasmas 7 (2000) 2076.
[6] A.C.L. Chian, C.F. Kennel, Astrophys. Space Sci. 97 (1983) 9.
[7] R.E. Kates, D.J. Kaup, J. Plasma Phys. 41 (1989) 507.
[8] R.T. Gangadhara, V. Krishan, P.K. Shukla, Mon. Not. R. Astron. Soc. 262(1993) 151.
[9] T. Tajima, T. Taniuti, Phys. Rev. A 42 (1990) 3587;
K. Holcomb, T. Tajima, Phys. Rev. D 40 (1989) 3909.
[10] V.I. Berezhiani, S.M. Mahajan, Phys. Rev. Lett. 73 (1994) 1110;
V.I. Berezhiani, S.M. Mahajan, Phys. Rev. E 52 (1995) 1968.
[11] T. Tatsuno, V.I. Berezhiani, S.M. Mahajan, Phys. Rev. E 63 (2001)
046403;
T. Tatsuno, V.I. Berezhiani, M. Pekker, S.M. Mahajan, Phys. Rev. E 68
(2003) 016409.
[12] M.J. Rees, astro-ph/9701162.
[13] V.V. Usov, Nature 357 (1992) 472.
[14] C. Thompson, Mon. Not. R. Astron. Soc. 270 (1994) 480.
[15] E.G. Blackman, L. Yi, G.B. Field, Astrophys. J. 473 (1996) L79.
[16] P. Mészáros, M.J. Rees, Astrophys. J. 482 (1997) L29.
[17] W. Kluzniak, M.A. Ruderman, Astrophys. J. 505 (1998) 113.
[18] M. Lyutikov, E.G. Blackman, Mon. Not. R. Astron. Soc. 321 (2001) 177.
[19] V.I. Berezhiani, S.M. Mahajan, Z. Yoshida, M. Ohhashi, Phys. Rev. E 65
(2002) 047402.
[20] G.Z. Sun, E. Ott, Y.C. Lee, P. Guzdar, Phys. Fluids 30 (1987) 526.
[21] N.G. Vakhitov, A.A. Kolokolov, Izv. Vyssh. Uchebn. Zaved. Radiofiz. 16(1973) 1020, Sov. Radiophys. 9 (1973) 262.
[22] V. Skarka, V.I. Berezhiani, R. Miklaszewski, Phys. Rev. E 56 (1997) 1080.
[23] V.E. Zakharov, V.V. Sobolev, V.C. Synakh, Zh. Eksp. Teor. Fiz. 60 (1971)
136, Sov. Phys. JETP 33 (1971) 77.
[24] V.E. Zakharov, A.M. Rubenchik, Zh. Eksp. Teor. Fiz. 65 (1973) 997, Sov.
Phys. JETP 38 (1974) 494.
[25] S.M. Mahajan, V.I. Berezhiani, R. Miklaszewski, Phys. Plasmas 5 (1998)
3264.
[26] P.K. Shukla, M. Marklund, B. Eliasson, Phys. Lett. A 324 (2004) 193.