The Astrophysical Journal, 697:1681–1694, 2009 June 1
C 2009. The American Astronomical Society. All rights reserved. Printed in the U.S.A.

doi:10.1088/0004-637X/697/2/1681

STABILITY OF RELATIVISTIC FORCE-FREE JETS
Ramesh Narayan1, Jason Li2, and Alexander Tchekhovskoy1
1 Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA; rnarayan@cfa.harvard.edu 2 Department of Astrophysical Sciences, Princeton University, Peyton Hall, Ivy Lane, Princeton, NJ 08544, USA Received 2009 January 28; accepted 2009 March 24; published 2009 May 13
ABSTRACT
We consider a two-parameter family of cylindrical force-free equilibria, modeled to match numerical simulations of relativistic force-free jets. We study the linear stability of these equilibria, assuming a rigid impenetrable wall at the outer cylindrical radius Rj. Equilibria in which the Lorentz factor γ (R) increases monotonically with increasing radius R are found to be stable. On the other hand, equilibria in which γ (R) reaches a maximum value at an intermediate radius and then declines to a smaller value γj at Rj are unstable. A feature of these unstable equilibria is that poloidal field line curvature plays a prominent role in maintaining transverse force balance. The most rapidly growing mode is an m = 1 kink instability which has a growth rate ∼ (0.4/γj )(c/Rj ). The e-folding length of the equivalent convected instability is ∼2.5γj Rj . For a typical jet with an opening angle θj ∼ few/γj , the mode amplitude grows only weakly with increasing distance from the base of the jet. The growth is much slower than one might expect from a naive application of the Kruskal–Shafranov stability criterion.
Key words: galaxies: jets – instabilities – MHD

1. INTRODUCTION
Relativistic jets in astrophysical sources have been known for many decades. Although their enormous power, large Lorentz factor, and strong collimation have been widely studied, these phenomena still lack an accepted explanation. An even greater mystery is the remarkable coherence and apparent stability of jets over very large length scales. This is the topic of the present paper.
The most promising models of relativistic jets involve acceleration and collimation by magnetic fields with footpoints attached to a spinning black hole/neutron star or an accretion disk. The forced rotation of the field lines induces a strong toroidal component of the magnetic field, which is responsible for accelerating the jet (e.g., Narayan et al. 2007; Tchekhovskoy et al. 2008, hereafter TMN08; and references therein). In this picture the toroidal component of the field dominates over other field components.
According to the well known Kruskal–Shafranov (KS) criterion (e.g., Bateman 1978), cylindrical magnetohydrodynamic (MHD) configurations in which the toroidal field dominates are violently unstable to the m = 1 kink instability (also called the screw instability). The KS criterion for instability is

Bφ > 2π Rj , Bp z

(1)

where Bφ and Bp are the toroidal and poloidal magnetic field strengths, Rj is the cylindrical radius, and z the length of the system. Typical jet models, including those described in this paper (see Section 2.2), have Bφ ∼ γj Bp, where γj 1 is the Lorentz factor of the jet, and they have jet angles θj ∼ Rj /z ∼ few/γj . Substituting these scalings in Equation (1), it is obvious that the KS instability criterion is
easily satisfied in relativistic jets. We therefore expect astro-
physical jets to be violently unstable, as argued for example by
Begelman (1998) and Li (2000). Yet, jets in nature are appar-
ently quite stable. How is this possible?
Many authors have investigated this question. They have
used jet models with a variety of velocity profiles, geometrical

shapes, composition, and boundary conditions (Kadomtsev 1966; Bateman 1978; Ferrari et al. 1978; Benford 1981; Payne & Cohn 1985), and applied both analytical and numerical methods (Istomin & Pariev 1996; Begelman 1998; Lyubarskii 1999; Li 2000; Lery et al. 2000; Appl et al. 2000; Tomimatsu et al. 2001; Mizuno et al. 2007; Moll et al. 2008; McKinney & Blandford 2009). As a result of this large body of work, several kinds of unstable modes have been identified: reflection modes, Kelvin–Helmholtz modes, current-carrying modes, etc. Unfortunately, it is difficult to synthesize the results to extract universal principles.
A fruitful approach in this field is to reduce relativistic jet models to their barest minimum. One such approach is to consider force-free jet models in which one ignores the inertia and pressure of the plasma and considers only charges, currents, and fields. The force-free approximation is valid whenever the energy density in fields dominates over matter energy density, as in pulsar magnetospheres (Goldreich & Julian 1969; Ruderman & Sutherland 1975). The force-free approximation is valid also in relativistic MHD jets, at least inside the fast magnetosonic surface (e.g., Tchekhovskoy et al. 2009; Lyubarsky 2009).
Theoretical studies of force-free jets have led to the identification of two distinct stability criteria. In a detailed analysis, Istomin & Pariev (1996) showed that cylindrical force-free jets in which Bz is independent of R are stable. Lyubarskii (1999) considered models with nonconstant Bz and showed that forcefree jets are unstable if

dBz < 0, dR

(2)

i.e., if the poloidal field decreases with increasing distance from

the axis. We refer to Equation (2) as the IPL criterion for

instability. On the other hand, Tomimatsu et al. (2001; TMT)

showed via an approximate analysis3 that force-free jets are

unstable if

Bφ

|Ω|R >,

Bp c

(3)

3 They effectively restricted their analysis to the region of the jet inside the light cylinder. Therefore, the flow speeds they considered were at best only quasi-relativistic.

1681

1682

NARAYAN, LI, & TCHEKHOVSKOY

Vol. 697

where Ω is the angular velocity of the field line. This criterion— the TMT criterion—differs from the KS criterion (1) in that it explicitly accounts for rotation. It is also apparently very different from the IPL criterion.
We describe, in this paper, a class of force-free cylindrical jet equilibria that closely match the numerical force-free jet simulations reported in TMN08, including the important effect of poloidal field curvature. Within the context of force-free jets from rigidly rotating stars, we believe that this two-parameter family of equilibria is generic and fairly complete. We study the stability properties of these equilibria and attempt to relate our results to the KS, IPL, and TMT criteria (Equations (1)–(3)).
In Section 2, we summarize the numerical simulation results of TMN08 (Section 2.1) and we describe an analytical forcefree jet model that matches the simulation data very closely (Section 2.2). In Section 3, we carry out a linear stability analysis of these equilibria and show that the linear modes of the system are obtained by solving an eigenvalue problem involving two coupled differential equations, with appropriate boundary conditions. In Section 4, we numerically solve the equations and identify the unstable modes in the system. We also derive an approximate estimate for the growth rate of the instability. We conclude in Section 5 with a summary and discussion.
We use (r, θ, φ) for spherical coordinates and (R, φ, z) for cylindrical coordinates. The idealized jet equilibria we consider are axisymmetric and cylindrical, with rigid wall boundary conditions at radius R = Rj . Each equilibrium is defined in terms of two parameters, Rm and γm, as described in Section 2.2. When Rm/Rj < 1, the Lorentz factor γ reaches a maximum value equal to γm at R = Rm and decreases to a smaller value γj at R = Rj . On the other hand, when Rm/Rj > 1, γ increases monotonically out to R = Rj . These two types of models are described in Sections 2.2.1 and 2.2.2, and Figures 1 and 2, and their stability properties are studied in Section 4.

2. FORCE-FREE JET EQUILIBRIUM

2.1. Structure of Force-Free Jets

TMN08 considered a rigidly rotating star of unit radius (r = 1) surrounded by a differentially rotating infinitely thin disk extending from R = 1 outward. The star was threaded by a
uniform radial magnetic field Br, and the disk was threaded by a power-law distribution of vertical field,

Bz ∝ Rν−2.

(4)

Using a relativistic force-free code (Gammie et al. 2003;

McKinney & Gammie 2004; McKinney 2006; Mignone &

McKinney 2007; Tchekhovskoy et al. 2007), TMN08 numeri-

cally evolved the system and obtained the equilibrium configu-

ration of the magnetic field.

Following TMN08, we will call the field lines that emerge

from the star as the “jet” and the field lines from the disk as the

“wind.” The critical field line that emerges from the star–disk

boundary defines the boundary between the jet and the wind.

This boundary starts off at θ = π/2 at the surface of the star

(r = R = 1) but decreases to smaller values of θ with increasing

r (or z).

We are primarily interested in the “jet”—the bundle of field

lines attached to the central star. Since all of these field lines

rotate at the angular velocity Ω of the star, the Alfve´n surface for

these lines takes the form of a cylinder—the “light cylinder”—

with radius

RA = c/Ω.

(5)

Field lines become strongly toroidal once they are outside the Alfve´n surface, which is where most of the collimation and acceleration occurs (TMN08).
TMN08 showed that the structure of the jet is strongly affected by the radial pressure profile of the region surrounding the jet. Specifically, if we write the radial variation of the confining pressure as p ∝ r−α, the jet properties are determined by the value of α. In the numerical experiments, the pressure was caused by a force-free disk wind, and α was determined by the index ν (Equation (4)) according to

α = 2(2 − ν).

(6)

At distance z along the axis, the cylindrical radius Rj of the jet is approximately given by

Rj ∼ zα/4.

(7)

For all α < 4, the jet collimates as it moves away from the star (Lynden-Bell 2006). As a result, at a sufficiently large distance from the central star (r 1), the jet is nearly cylindrical in shape.
In the asymptotic nearly cylindrical region of the jet, force balance in the R direction is described by the following equation (TMN08; see Equation (23) below for the cylindrical version of this equation):

d B2 − E2 + Bφ2 − E2 + Bp2 − E2 = 0, (8)

dR 8π

4π R

4π Rc

where B is the total magnetic field strength, Bφ is the toroidal field strength, and Bp is the poloidal field strength:

B = Bp2 + Bφ2,

Bp = BR2 + Bz2.

(9)

The electric field E is given by

E = −(ΩR/c) φˆ × B,

(10)

where Ω is the angular velocity of the field line. The electric field has only a poloidal component: Ep = ΩRBp/c.
Each of the three terms on the left-hand side of Equation (8) represents a force in the −R direction. The quantity (B2 − E2)/8π is the pressure of the force-free fluid in the comoving
frame; therefore, the first term describes the inward force due
to the gradient of pressure. The second term arises from the
toroidal curvature of the field line. The toroidal magnetic field
Bφ contributes an inward force due to “hoop stress,” while the poloidal electric field E contributes an outward force.4 The third
term in (8) gives analogous contributions from the poloidal
curvature of the field line, where Rc is the poloidal radius of curvature; once again there is an inward force due to the poloidal
magnetic hoop stress and an outward force due to the electric
field. Note that the contributions involving E are important
only for relativistic flows. In standard nonrelativistic MHD, one
neglects these terms and keeps only the terms involving B.
TMN08 showed that models with α > 2, i.e., ν < 1, are good
analogs of relativistic jets found in nature, especially gamma-
ray burst jets. The main features of the numerical jet solutions

4 As Equation (10) shows, E is directly proportional to Ω, so the latter term results from rotation and may loosely be viewed as a sort of “centrifugal force” (V. Beskin 2007, private communication).

No. 2, 2009

STABILITY OF RELATIVISTIC FORCE-FREE JETS

1683

are as follows (see Figure 1 below). (1) Bp is essentially constant inside the jet, showing hardly any variation with R. As we show
in Appendix A, this is required for force-free jet solutions that smoothly connect to the central compact object.5 (2) Since Ω is constant and E ∝ ΩRBp (Equation (10)), we have E2 ∝ R2. (3) Bφ is almost equal to E and so Bφ also varies primarily as R2. (4) Bφ is slightly larger than E with Bφ2 − E2 ∝ R4. Because of this property, the first two terms in Equation (8)
both represent inward forces. Therefore, (5) the third term in
Equation (8), which involves the poloidal curvature of the field
line, is important for force balance. This is the only outward
force in the balance equation—it is outward because E is of
order Bφ and is much greater than Bp outside the light cylinder. This force has to balance the other two terms.
The velocity of a force-free flow is usually identified with the
drift velocity,

v c

=

E×B B2 ,

(11)

2.2. Analytical Jet Model
The axisymmetric numerical jet models described in Section 2.1 have magnetic and electric field components that are functions of both R and z. This is not convenient for linear perturbation analysis. Since the numerical models are nearly cylindrical at a large distance (i.e., dR/dz 1), we now consider an idealized jet equilibrium model which is perfectly cylindrical and in which all quantities are functions only of R. We choose the following specific functional forms:

B0R = 0,

(16)

B0φ = − 2 γm2 − 1 (R/Rm)2 + (R/Rm)4 1/2 ≡ −f (R), (17)

B0z = exp −3R2/4 γm2 − 1 Rm2 ≡ g(R),

(18)

and the Lorentz factor is defined correspondingly. Since E ·B = 0, it is easily shown that

E0R = − 2(γm2 − 1) 1/2 (R/Rm) ≡ −h(R), (19)

γ2

=

B2 .
B2 − E2

TMN08 derived two approximate relations for γ ,

(12)

E0φ = 0, E0z = 0,

(20) (21)

γ1 = [1 + (ΩR/c)2]1/2,

(13)

γ2 = (3Rc/R)1/2,

(14)

and they showed that the net γ of the fluid is given to good accuracy by the following simple formula (this result is derived later; see Equation (29)):

1 11 γ 2 = γ12 + γ22 .

(15)

Along each field line, γ is initially determined mainly by rotation, and so γ ≈ γ1. This is a region of efficient acceleration which TMN08 called the first acceleration regime. However,
beyond a certain distance from the star, the effect of poloidal
field line curvature becomes important, and γ switches to the
less efficient γ2, the second acceleration regime. Relatively near the star, all field lines in the jet are
in the first acceleration regime and γ (R) behaves like γ1 (Equation (13)). Specifically, γ increases more or less linearly
with R and reaches its maximum value at the edge of the jet at R = Rj . However, when we consider the jet at a larger distance from the star, some of the field lines have already switched to the second acceleration regime, where γ ∼ γ2 ∝ 1/R (see Equation (14), coupled with Equation (22) below). In this
case, the maximum Lorentz factor occurs at some radius Rm inside the jet, not at the boundary; we have γ ∼ γ1 ∝ R for R Rm and γ ∼ γ2 ∝ R−1 for Rm R < Rj . TMN08 discuss in detail the physics of the two acceleration
regimes.

5 Asymptotic force-free jet configurations with nonconstant profiles of Bp are certainly possible (Lyubarskii 1999), but there exists no solution that would
smoothly connect them to the compact object.

Rc = 2 γm2 − 1 Rm2 /3R.

(22)

The zeros in the subscripts are meant to indicate that all these
quantities refer to the unperturbed model. The model has two
parameters, γm and Rm, whose meanings are explained below. For simplicity, we have chosen units such that B0z = 1 at the jet axis (R = 0). Also, we have assumed that B0z and Ω are positive, so both B0φ and E0R are negative, i.e., magnetic field lines are swept backward with respect to the rotation and the
electric field is pointed radially inward. With this choice of
signs, the three functions f (R), g(R), and h(R) are positive.
Note that, in all cases of interest, g(R) is practically equal to
unity. The particular exponential form given in Equation (18) is
designed to handle small higher-order terms in the force balance
Equation (23), but the deviations of g(R) from unity are tiny and unimportant.6
By direct substitution it is easily verified that the above model
satisfies the radial force balance Equation (8). Under cylindrical
symmetry, this equation takes the form

d B02φ + B02z − E02R + B02φ − E02R + B02z − E02R = 0.

dR 8π

4π R

4π Rc

(23)

The first two terms are positive, i.e., both represent inward

forces, with the first term providing twice as much force as

the second. The third term is negative and its magnitude is equal

to the sum of the other two terms.

An important feature of the above model is that the outward

force from the third term involves the poloidal curvature radius

Rc. Technically, a perfectly cylindrical model has Rc → ∞. To

get around this problem, we treat Rc as an externally imposed

6 As a test, in the stability analysis described later we have done the
calculations both with the full expression for g(R) given in Equation (18) and with the simpler choice g(R) = 1. The results are practically the same.

1684

NARAYAN, LI, & TCHEKHOVSKOY

Vol. 697

property of the solution which is adjusted so as to reproduce the poloidal curvature force present in the numerical jet model. In other words, even though we have straightened out field lines in the z-direction by enforcing cylindrical geometry, we still retain the effect of poloidal curvature by means of an artificial external force. This procedure is analogous to the widely used shearing sheet approximation in accretion disk studies (Goldreich & Tremaine 1978; Narayan et al. 1987) in which fluid streamlines are straightened out in the azimuthal direction, but the effect of azimuthal curvature is still retained via a Coriolis force. Note that, apart from the extra term due to poloidal field curvature, Equation (8) is identical to the standard balance condition derived in other papers in the literature, e.g., Equation (6) in Istomin & Pariev (1996) or Equation (16) in Lyubarskii (1999).
To get a better idea of the nature of the above analytical model, we now make a couple of simplifications. As already mentioned, B0z is practically independent of R inside the jet. Also, for highly relativistic jets, we invariably have Bφ2 −E2 Bφ2. We therefore replace Equations (17), (18), and (19) by the following simpler formulae:

B0z ≈ 1,

(24)

B0φ ≈ E0R = − 2 γm2 − 1 1/2 (R/Rm),

(25)

B02φ − E2 = (R/Rm)4.

(26)

By Equation (10), the angular velocity of rotation of the field lines is given by

Ω = −cE0R/B0zR ≈ 2 γm2 − 1 (c/Rm),

(27)

and is the same for all field lines, as required for a rigidly rotating star at the base of the jet.7 Using these simpler expressions, we
obtain the following result for the Lorentz factor:

1 γ 2(R)

=

B02 − E02 B02

=

B02φ + B02z − E02R B02φ + B02z

≈ 1 + (R/Rm)4 1 + 2 γm2 − 1 (R/Rm)2

(28)

≈

1

+

1 (ΩR/c)2

+

R 3Rc

,

(29)

where we have made use of Equations (22) and (27). We thus
reproduce the result given earlier in Equation (15).
The simple analytical model described here can reproduce
all the features seen in numerical simulations of force-free jets.
Figure 1 shows numerical results obtained by TMN08 for their “fiducial” model (ν = 0.75 or α = 2.5). Panels (a) and (b) correspond to a relatively near region of the jet at z = 102, where the jet is entirely in the first acceleration regime. Panels (c) and (d) correspond to a more distant region at z = 107, where part of the jet has made the transition to the second acceleration
regime (the region where γ decreases with increasing R). In each
panel, the abscissa gives the cylindrical radius R normalized by

7 Since we have designed our analytical model to match the numerical models of TMN08, all of our models have constant Ω(R). It would be straightforward to generalize the model to nonconstant Ω(R), using additional parameters.

8 (a)
6
4
2
0
-2
-4 -2 -1.5 -1 -0.5

0 0.5

3 (b)
2
1
0 -2 -1.5 -1 -0.5 0 0.5

10 (c)
8
6
4
2
0
-2 -2 -1.5 -1 -0.5

0 0.5

5 (d)
4
3
2
1 -2 -1.5 -1 -0.5 0 0.5

Figure 1. Numerical results from a force-free jet simulation with α = 2.5 (“fiducial model” corresponding to ν = 0.75 in TMN08). Panel (a): field components as functions of normalized cylindrical radius R/Rj at z = 102. The solid horizontal line shows Bp2 (R)/Bp2 (0), the short-dashed line shows [Bφ2(R) − E2(R)]/Bp2 (0) and the long-dashed line shows E2(R)/Bp2 (0). The dotted lines are the corresponding results for the analytical model with γm = 6, Rm = 1.2 (Model A, Section 2.2.1). The vertical solid line shows the boundary between the jet and the external confining medium. Panel (b): Lorentz factor
γ (R) (solid line), and the two approximations, γ1(R) (short-dashed line) and γ2(R) (long-dashed line), at the same z. The dotted lines are from the analytical model. Panels (c) and (d): similar to (a) and (b), but at z = 107. The dotted lines in these panels correspond to the analytical model with γm = 2600, Rm = 0.18 (Model C, Section 2.2.2).

the local “jet radius” Rj, which is the cylindrical radius of the field line that separates the jet from the surrounding disk wind.
In the analytical model, it is easily shown that γ reaches a maximum at R = Rm and that its value at this radius is equal to γm. Thus, the two model parameters Rm and γm allow us to control the basic features of the equilibrium. To obtain a jet that is entirely in the first acceleration regime (similar to Figures 1(a) and (b)), we must choose Rm > Rj , while for a jet that is partly in the second acceleration regime (Figures 1(c) and (d)), we require Rm < Rj . We consider these two cases in the following subsections.
2.2.1. Rm > Rj : Maximum Lorentz Factor Located at the Jet Boundary
A jet in which all field lines are in the first acceleration regime has its maximum Lorentz factor at the boundary of the jet, R = Rj . This corresponds to choosing Rm > Rj in the analytical model, so that the term (R/Rm)4 in the numerator of Equation (28) can be neglected. In this case, the profile of γ has two segments: the region of the jet inside the light cylinder (R < RA = c/Ω) which is not accelerated very much, and the region outside the light cylinder which has γ increasing linearly with radius,

γ (R) ≈ 1, R/RA ≈ (R/Rj )γj ,

R < RA = Rm/ 2 γm2 − 1 , RA < R < Rj ,
(30)

No. 2, 2009

STABILITY OF RELATIVISTIC FORCE-FREE JETS

1685

3. LINEAR PERTURBATION ANALYSIS

We now consider linear perturbations of the cylindrical equilibrium described in Section 2.2. The unperturbed state has magnetic and electric fields

B0 = B0RRˆ + B0φφˆ + B0zzˆ,

(37)

E0 = E0RRˆ + E0φφˆ + E0zzˆ,

(38)

where the various components are given by the expressions
in Equations (16)–(21). As mentioned previously, we choose B0z and Ω to be positive, so B0φ and E0R are negative. The unperturbed current and electric charge are

J0

=

c 4π

∇

×

B0

=

c 4π

−

dB0z + B0z dR Rc

φˆ

+

1 R

d dR

(RB0φ )zˆ

,

(39)

Figure 2. Profiles of γ vs. R/Rj for the analytic Models A, B, and C.

where γj is the Lorentz factor at the jet boundary,
γj ≈ γ1(Rj ) = 1 + 2 γm2 − 1 (Rj /Rm)2 1/2 . (31)
The dotted lines in panels (a) and (b) in Figure 1 show the dependences of various quantities as a function of R/Rj for a model with γm = 6 and Rm = 1.2Rj . The agreement with the numerical simulation results at z = 102 is striking. We call the analytical model with this particular choice of γm and Rm Model A:

Model A :

γm = 6, Rm = 1.2Rj .

(32)

The solid line in Figure 2 shows the variation of γ as a function of R for this model.

2.2.2. Rm < Rj : Maximum Lorentz Factor Located Inside the Jet
As we described in Section 2.1, a jet in which some field lines have switched to the second acceleration regime has its maximum Lorentz factor inside the jet. This means Rm < Rj . In this case, the Lorentz factor γj at the jet boundary is roughly equal to

γj ≈ γ2(Rj ) = 2 γm2 − 1 1/2 Rm/Rj . Now the profile of γ has three segments:

(33)

γ (R) ≈

1, R < RA,

R/RA ≈ (R/Rm)γm,

RA < R < Rm, (34)

(Rm/R)γm ≈ (Rj /R)γj , Rm < R < Rj .

The dotted lines in panels (c) and (d) of Figure 1 show model results corresponding to γm = 2600 and Rm = 0.18Rj . We find excellent agreement with the numerical results at z = 107.
We call the analytical model with these values of γm and Rm Model C. For completeness we also consider a less extreme model called Model B in which γm = 6 and Rm = 0.3Rj :

Model B:

γm = 6,

Rm = 0.3Rj , (35)

Model C:

γm = 2600, Rm = 0.18Rj . (36)

The dashed and dotted lines in Figure 2 show the variations of γ as a function of R for these two models.

ρ0

=

1 4π

∇

·

E0

=

1 4π R

d dR (RE0R)

+

1 4π

E0R Rc

.

(40)

Note that we have included terms involving Rc in the unperturbed current and charge density. These terms describe the
contributions of poloidal field curvature to the quantities ∇ × B0 and ∇ · E0, respectively. By including these terms, we retain the forces associated with poloidal field curvature without actually having curved field lines in the model.

3.1. The Eigenvalue Problem
We now consider small perturbations. Let us write the perturbed electric field as E = E0 + E1, where E1 is a small perturbation of the form

E1 = [E1R(R)Rˆ + E1φ(R)φˆ + E1z(R)zˆ] exp(−iωt + imφ + ikz). (41)
Let us similarly write B = B0 + B1, J = J0 + J1, ρ = ρ0 + ρ1. Each of these small perturbations can be expressed in terms of the perturbed electric field via Maxwell’s equations. From

1 ∂B = −∇ × E, c ∂t

(42)

we obtain From we obtain

ic B1 = − ω ∇ × E1. 1 ∂E = ∇ × B − 4π J , c ∂t c

(43) (44)

J1

=

iω 4π

E1

−

ic2 ∇ 4π ω

×

(∇

×

E1).

(45)

Finally, from

∇ · E = 4πρ,

(46)

we obtain

ρ1

=

1 4π

∇

·

E1.

(47)

Since the perturbed system is force-free, it must satisfy

E · B = 0. The zeroth-order terms satisfy this trivially (as they

1686

NARAYAN, LI, & TCHEKHOVSKOY

Vol. 697

should). The first-order terms give the condition E1 · B0 + E0 · B1 = 0. Substituting for the various quantities, this condition allows
us to solve for E1φ in terms of E1z:

E1φ = C1E1z,

(48)

where the function C1 is given by

C1

=

ωRg − cmh .
ωRf − ckRh

(49)

Successive differentiations give

appropriate boundary conditions, we obtain ω for given values of m and k.
The singular points of the equations are located at the radii where D1(R) and D4(R) vanish. Anticipating later discussion, we write down here the expression for the quantity D1D4:
g D1D4 = − ωR
(ωRf −ckRh)2 +(ωRg−mch)2 −(ckRg−mcf )2 ×.
(ωRf −ckRh)
(58)

E1φ = C1E1z + C1E1z,

(50) Also, from Equation (11), the perturbed velocity is

E1φ = C1E1z + 2C1E1z + C1E1z,

(51)

where primes refer to derivatives with respect to R.
We now consider the force balance condition: ρE+(1/c)J ×B = 0. The zeroth-order terms give

1 ρ0E0 + c J0 × B0 = 0,

(52)

which is simply the equilibrium force balance condition (23).
Notice that the poloidal curvature terms in J0 and ρ0 are necessary to satisfy equilibrium in the unperturbed solution. From the first-order terms in the force balance equation we obtain

ρ1E0

+

ρ0E1

+

1 c J1

×

B0

+

1 c J0

×

B1

=

0.

(53)

The φˆ component of this equation gives a relation between E1φ, E1φ, E1z, and E1z. Eliminating E1φ using Equation (50), we obtain a first-order differential equation for E1z(R):

D1E1z + D2E1z + D3E1R = 0,

(54)

where D1, D2, and D3 are functions of R. The expressions are relatively long and we give them in Appendix B. The zˆ component of Equation (53) has no new information; it just gives back Equation (48). The Rˆ component, however, gives
a new relation between the various components of E1 and their derivatives. Eliminating E1φ, E1φ, E1z, and E1z using Equations (50), (51), (54), and the derivative of (54), we obtain
a differential equation for E1R(R):

D4E1R + D5E1z + D6E1R = 0,

(55)

where D4, D5, and D6 are again functions of R and are given in Appendix B.
We have thus reduced the linear mode analysis problem to
a pair of first-order differential equations, Equations (54) and
(55), for E1z(R) and E1R(R). For convenience, we write down the two equations again:

E1z

=

−

D2 D1

E1z

−

D3 D1

E1R

,

(56)

E1R

=

−

D5 D4

E1z

−

D6 D4

E1R

.

(57)

These equations constitute an eigenvalue problem, where ω is the eigenvalue. By numerically solving the equations with

v1 c

=

E1

× B0 + E0 B2

× B1 .

(59)

In component form this gives

v1R c

=

gE1φ f2

+ +

f E1z g2

,

(60)

v1φ c

=

−

gE1R f2

− hB1z + g2

,

(61)

v1z c

=

−

f

E1R f2

+ +

hB1φ g2

.

(62)

We now consider boundary conditions. A physically valid perturbation will be well behaved on the axis (R = 0) and will satisfy suitable boundary conditions at the jet boundary (R = Rj ). The condition on the axis is different for axisymmetric (m = 0) and nonaxisymmetric (|m| 1) perturbations, so we consider each of these cases in turn. At the jet boundary, we
assume that the jet is constrained by a “rigid wall” and we write
down the corresponding boundary condition. In the following,
we employ the specific forms of f (R), g(R), and h(R) given in
Equations (17), (18), and (19).

3.2. Boundary Condition on the Axis: m = 0

Setting m = 0 and substituting the expressions for f (R), g(R), and h(R) in D1 − D6, we find that the leading terms of the differential equations, Equations (56) and (57), at small R
are given by

E1z

=

azz R

E1z

+

azR E1R ,

(63)

E1R

=

aRz R2

E1z

+

aRR R

E1R

,

(64)

where

2ω

azz =

−

, ck

(65)

azR

=

i(c2k2 − c2k

ω2) ,

(66)

4iω

aRz

=

k(ck

−

, ω)

(67)

(ck + 2ω)

aRR =

. ck

(68)

No. 2, 2009

STABILITY OF RELATIVISTIC FORCE-FREE JETS

Requiring the perturbation to be analytic as R → 0 immediately gives the following solution near the axis:

E1z = KR2,

(69)

2ic E1R = − K (ck − ω) R, where K is an arbitrary normalization constant.

(70)

3.3. Boundary Condition on the Axis: |m| > 0 When m = 0, we obtain

azz = 0,

(71)

1687

azR

=

imA(ck − ω) ,
(Acm − ωRm)

aRz

=

−

im(Acm − ωRm) , A(ck − ω)

(72) (73)

aRR = − 1,

(74)

where the constant A is defined to be A = 2 γm2 − 1 1/2 .

(75)

The physically relevant solution close to the axis is then

E1z = KR|m|,

(76)

E1R

=

− K sgn(m) i(Acm − ωRm) A(ck − ω)

R|m|−1.

(77)

3.4. Boundary Condition at the Jet Boundary: Rigid Wall
We assume that our cylindrical jet is terminated at R = Rj by a rigid impenetrable wall. By impenetrable we mean that no energy flows across this boundary, i.e., the Poynting flux lies in the φ − z plane. Equivalently, the velocity vector has no radial component.
The equilibrium Poynting flux of course lies in the φ−z plane.
The perturbed Poynting flux is proportional to E1×B0 + E0×B1. Since E0 is parallel to Rˆ, the term E0 × B1 is automatically in the φ − z plane. The term E1 × B0 will also be in this
plane if E1 is precisely radial, i.e., both E1φ and E1z vanish. By Equation (48), E1φ is proportional to E1z. We thus obtain the following boundary condition at the outer wall:

E1z = 0,

R = Rj .

(78)

4. NUMERICAL RESULTS
We have computed frequencies of modes by numerically solving the differential equations, Equations (56) and (57), along with the boundary conditions described in Sections 3.2–3.4. For each choice of k and m, a countable infinity of solutions exists which may be ordered by the number of zeros of Re[E1z(R)], not counting zeros at the boundaries.8 The lowest-order solution (the “fundamental mode”) is such that Re[E1z(R)] has no zeros between R = 0 and R = Rj , the next solution has one zero

8 Re() stands for the real part of a complex quantity.

Figure 3. Dispersion relation for axisymmetric modes (m = 0) in Models A (solid lines), B (dashed lines), and C (dotted lines). The curves corresponding to radial mode numbers n = 0, 1, and 2 are labeled. All the modes are stable.
inside the jet, and so on. In the following we identify each mode by its radial mode number n which we define to be the number of zeros. As one might expect, the mode frequency increases with increasing n.
We solve for the frequencies via a shooting method. We start with a guess value of ω, make use of the expressions given in Section 3.2 or Section 3.3 (depending on the value of m) to set up the initial solution at small R, and integrate Equations (56) and (57) to R = Rj . We then adjust ω in the complex plane until the outer boundary condition given in Section 3.4 is satisfied. The only subtle point is that the quantities D1 and D4 appear in the denominators of various coefficients in Equations (56) and (57), and so their zeros correspond to singularities in the solution. To avoid these singularities, we treat R as a complex variable and integrate the equations over a “safe” trajectory in complex-R space. Since the solution is analytic, the exact track that we follow is unimportant so long as it lies above all singularities in the R-plane. Istomin & Pariev (1996) give a detailed discussion of this point in connection with currentdriven instabilities in force-free jets. The reader is also referred to standard discussions of the topic in plasma physics texts in the context of Landau damping, or Goldreich et al. (1986) for a discussion in the context of accretion disk instabilities.
4.1. Axisymmetric Modes: m = 0
Figure 3 shows the dispersion relation—the variation of ω with k—of axisymmetric modes (m = 0). Results are shown for Models A, B, and C (Equations (32), (35), and (36)) for three radial mode numbers: n = 0, 1, and 2. For all k, we find that the mode frequency is real, which means that all these modes are stable.
For large values of kRj, the mode frequency asymptotes to ω = ±ck, so the modes behave like electromagnetic waves moving parallel or antiparallel to the z-axis. At small k, however, the frequency asymptotes to a constant value. There is thus a minimum frequency for propagating modes inside the jet. The minimum frequency is of order the inverse of the lightcrossing time across a radial wavelength of the mode (e.g., ωmin ∼ 2π c/Rj for the mode with n = 0).

1688

NARAYAN, LI, & TCHEKHOVSKOY

Vol. 697

Figure 4. Imaginary part of ω for modes in Model B with m = 1 and n = 0. Growing modes have Im(ω) > 0, while decaying modes have Im(ω) < 0.

We find that the dispersion relations of modes with positive and negative k are not quite the same. The difference arises because the background has a nonzero velocity in the z-direction, which breaks the symmetry between waves propagating toward +z and −z. The effect is, however, quite weak.

4.2. Nonaxisymmetric Modes: m = ±1

The most interesting modes are those with m = ±1. These

modes are stable in Model A, but unstable in Models B and C.

Figure 4 shows Im(ω)9 as a function of kRj for a sequence of modes in Model B; the modes correspond to m = 1 and n = 0. In

this sequence, modes with kRj < 0.65 and those with kRj > 28 are stable and have ω real. However, for 0.65 < kRj < 28, we

find a pair of modes with complex values of ω. The branch with

Im(ω) > 0 corresponds to growing modes, and the branch with

Im(ω) < 0 to decaying modes.10

Figure 5 shows eigenfunctions corresponding to a few of

the growing modes. Plotted is Re(E1z) as a function of the

scaled radius R/Rj . The mode corresponding to kRj = 5 is

representative of all modes with kRj

5. These modes

have eigenfunctions with no zero crossings between R = 0

and R = Rj . By our definition, the modes correspond to n = 0.

Each of the remaining eigenmodes in Figure 5 has a pronounced

dip in Re(E1z) which causes a zero crossing. These dips result

from a singularity in the equations, as we discuss below. If

we discount the singularity-induced zero crossings, then these

eigenfunctions may also be identified as n = 0 modes.

Figures 6 and 7 show similar results for Model C. Growing

modes (and their decaying counterparts) are present for m = 1

and all kRj < 2.1 × 104. Modes with m = −1 are also unstable (see Figure 6).11 Figure 7 shows a few eigenfunctions. All modes

with kRj 1.3 × 103 have eigenfunctions with the standard n = 0 shape (see the mode with kRj = 1.25 × 103). For larger

values of k, the eigenfunctions develop negative spikes due to

the presence of a singularity. However, we still view them as

n = 0 modes.

9 Im() refers to the imaginary part of a complex quantity. 10 We discuss at the end of Section 4.2 the reason for the sudden jump in the
value of Im(ω) in the decaying branch. 11 In the case of Model B, modes with m = −1 appear to be stable, and only the m = +1 modes show an instability.

Figure 5. Eigenfunctions corresponding to growing modes in Model B with m = 1, n = 0 and kRj = 5, 10, 15, 20, and 25. The real part of E1z is plotted.

Figure 6. Imaginary part of ω for growing modes in Model C. The solid line corresponds to modes with m = 1 and n = 0 and the dotted line corresponds to modes with m = −1 and n = 0.

We have determined numerically that the singularities which

cause the dips in the eigenfunctions are due to zeros in the
function D4(R) defined in Section 3.1. This function appears in the denominator of the differential equation, Equation (57), and hence its zeros behave like poles.12
Equation (58) gives the analytic form of the quantity D1D4. Since the modes of interest to us have Re(ω) very nearly equal to ck, let us substitute ω = ck in this equation. Then, setting

D1D4 equal to zero gives the following relation between the wavenumber k of the mode and the radius Rsing of the singularity:

kRsing

=

g(Rsing)

m[f (Rsing) + h(Rsing)] + sgn(m) g2(Rsing) + [f 2(Rsing) −

. h2(Rsing)]

(79)

Figure 8 shows the position of the singularity Rsing as a function of kRj for modes with m = ±1 in Models B and C, as
calculated with this equation. For comparison, the dots show

12 In contrast, although the function D1(R) appears in the denominator of Equation (56), its zeros correspond to removable singularities in the
differential equation.

No. 2, 2009

STABILITY OF RELATIVISTIC FORCE-FREE JETS

1689

Figure 7. Eigenfunctions corresponding to growing modes in Model C with m = 1, n = 0, and kRj = 1.25 × 103, 2.5 × 103, 5 × 103, 104, and 2 × 104.
The real part of E1z is plotted.

outward and the growth rate of the mode increases (Figure 4). At kRj ∼ 5, when the singularity reaches the wall, the growth rate is close to its maximum value. At yet smaller values of k,
the singularity moves outside the outer wall, but its presence
is still felt and there is continued instability. The growth rate,
however, decreases with decreasing k.
A similar pattern is seen in Model C. Unstable modes are present only for kRj 2 × 104. With decreasing k the growth rate increases and reaches its maximum approximately when the singularity reaches the jet boundary (Rsing = Rj ), which happens at kRj ∼ 1.3 × 103. In contrast to Model B, however, the growth rate remains large even for smaller values of k, and the instability survives as k → 0.
We finally discuss the peculiar behavior of Im(ω) in the
decaying branch of modes in Figure 4. As we mentioned earlier,
in numerically solving for the eigenvalue we must integrate the
differential equations, Equations (56) and (57), along a path in
the complex-R plane that lies above the poles in the solution. For
growing modes, the pole is located below the real R-axis. We can
therefore integrate along the real R-axis without any difficulty.
For decaying modes, however, the pole is above the real
R-axis and we must choose the integration path with care. If the
singularity has Re(Rsing) > Rj , i.e., the singularity is outside the jet, there is no problem and we can simply integrate along
the real axis. However, when 0 < Re(Rsing) < Rj , we have to deform the integration path. In our calculations, we integrated from R = 0 along a path with Im(R) = Re(R) until the point Im(R) = Re(R) = Rj and we then integrated down to R = Rj . The jump in Im(ω) in Figure 4 is the result of the singularity
moving into the jet. To the left of the break, the singularity
is located at R > Rj . Here the eigenvalues of the growing and decaying modes are complex conjugates of each other. However,
to the right of the break, the singularity has moved inside the jet
(R < Rj ) and now the complex conjugate symmetry is broken. We note that eigenfunctions and eigenvalues of decaying
modes are not very meaningful. This can be shown by analyzing
the initial value problem along the lines of Landau’s treatment
of plasma damping. The reader is referred to Istomin & Pariev
(1996) for a detailed discussion of this topic.

Figure 8. Location of the singularity Rsing as a function of kRj in Models B and C, calculated using Equation (79). The solid curves correspond to modes with m = 1 and n = 0 and dotted lines to modes with m = −1 and n = 0. The solid dots show the locations of minima in the eigenfunctions plotted in Figures 5 and 7.
the radii at which the functions Re[E1z(R)] reach their minima in the eigenfunctions plotted in Figures 5 and 7. The agreement between the analytical curve and the dots is excellent, showing that Equation (79) captures the physics of the singularity.
From Figure 8 we see that the singularity lies inside the jet (Rsing < Rj ) only for a finite range of k above a certain minimum value. For values of k smaller than this minimum, the singularity is outside the jet (for very small k it is well outside the jet). In the case of Model B, the singularity enters the jet from outside when kRj ∼ 5 and it disappears at R = 0 when kRj ∼ 28 (for m = +1). This is the primary range of k over which an unstable mode is present. At kRj ∼ 28, the singularity is barely present near the center of the jet and we have a very weakly growing mode. With decreasing k, the singularity moves

4.3. Why is m = ±1 Special?
We have not exhaustively explored modes with |m| > 1. However, in spot tests with various choices of m, n, and kRj in Models A, B, and C, we found all modes to be stable. We believe that, if at all, there are only weakly unstable modes for |m| > 1; there is no sign of the kind of vigorous instability described in Section 4.2 for modes with m = ±1. So why is m = ±1 special? The answer to this question is well known in the magnetic confinement literature (e.g., Bateman 1978). We discuss it briefly here for completeness.
Consider the radial component of the perturbed velocity v1R near the axis of the jet. Equation (60) gives the expression for v1R in terms of the perturbed electric field components E1z and E1φ, and Equation (48) shows the relation between these two field components. For small values of R near the axis, we have

g(R) ≈ 1, f (R) ≈ h(R) = A R , Rm
where the quantity A is defined in Equation (75), and

(80)

C1

≈

(ωRm − Acm) AR(ω − ck)

.

(81)

1690

NARAYAN, LI, & TCHEKHOVSKOY

Vol. 697

Consider first modes with m = 0. Equation (69) shows that E1z ≈ KR2 near the axis. Substituting this in Equation (60) and
using the other approximations given above, we find

v1R

≈

ωRm KR A(ω − ck)

+

O(R3).

(82)

By symmetry, the velocity goes to zero on the axis, and the flow consists of a simple radial divergence.
Consider next modes with m = 0. Using Equation (76) for E1z, we find

v1R

≈

(ωRm − Acm) KR|m|−1 A(ω − ck)

cos mφ

+

O(R|m|+1),

(83)

where we have included cos mφ to show the angular dependence of the mode. The leading term goes like R|m|−1, which corresponds to R0 when m = ±1. This means that the mode has a finite radial velocity, and hence a finite radial displacement, on the axis when |m| = 1. The cos φ dependence of v1R, coupled with the fact that v1φ has the same amplitude but a sin φ dependence, ensures that the velocity vector is unique and analytic at R = 0. By writing the velocity vector in Cartesian coordinates, it is easily seen that the complex phase of v1R determines the orientation of the velocity vector in the x–y plane. If we consider values of |m| 2, the velocity vanishes on the axis, just as in the case of m = 0.
This then reveals what is special about |m| = 1 modes. These are the only modes in which fluid perturbations communicate
across the axis and cause the jet to shift bodily across the axis. In modes with |m| = 1 the center of mass of the jet itself shifts into a spiral shape, which is the characteristic feature of the kink or screw mode. For all other values of |m|, the center of mass remains on the axis and the perturbations are concentrated on
the outside. In helical MHD configurations in the laboratory, the |m| = 1
kink mode is known to be highly unstable and to be the
greatest threat to the stability of equilibria (Bateman 1978).
Not surprisingly we see the same feature in our force-free jet
equilibria.

Figure 9. Numerically calculated growth rates of modes with m = 1 and Rsing = 1.1Rj . These modes have among the largest growth rates. The solid lines show results for a series of models with fixed γm and varying Rm/Rj . The dotted lines are the growth rates predicted by Equation (85). Note the very good agreement except near the top of the plot, where the models are nonrelativistic. The dashed lines are the numerical growth rates for modes with k = 0.
kRj = 4.4. Since Model B has γj = 2.5, Equation (84) predicts kRj ≈ 4.0, which is close. Similarly, the mode with the maximum growth rate in Model C has Re(ω) ≈ kRj = 1000, whereas Equation (85) with γj = 660 predicts kRj ≈ 1060. We see that the approximate formula (84) is quite good.
Our numerical results indicate that the growth rate Im(ω) of the fastest growing mode is proportional to Re(ω)/γj2. We also know that unstable modes are present only when Rm < Rj ; for instance, Model A with Rm = 1.2Rj has no unstable modes, whereas Model B with Rm = 0.3Rj and Model C with Rm = 0.18Rj both have unstable modes. With these clues in mind, we obtain the following empirical estimate for the growth rate of the fastest growing mode:

4.4. Growth Rate of the Instability

The growth rate of the fastest growing mode is a matter of
practical interest since it limits the lifetime of an unstable sys-
tem. As discussed in Section 4.2, for the models we have con-
sidered here the most unstable mode generally has a singularity close to the outer wall: Rsing ∼ Rj . Knowing this, we estimate here the fastest growth rate by assuming that the pole is located at Rsing = 1.1Rj . (We locate the singularity slightly outside the jet, since this speeds up the numerical integrations.)
Given an assumed value of Rsing, we can substitute this value in Equation (79) and make use of the expressions for f (R),
g(R), and h(R) given in Section 2.2. Recalling that Models B
and C are in the regime described in Section 2.2.2, we note that f 2(Rj ) − h2(Rj ) g2(Rj ). In addition, f and h are nearly equal to each other and γj is given by Equation (33). We then find that kRj ≈ 1.6γj . Also, the real part of the frequency is nearly equal to ck. Thus, we estimate

Mode with Rsing = 1.1Rj :

k ≈ 1.6γj , Rj

Re(ω) ≈ 1.6γj c . Rj (84)

These estimates should apply to the fastest growing mode. The

mode with the maximum growth rate in Model B has Re(ω) ≈

Mode with Rsing = 1.1Rj :

Im(ω) ≈ 0.4 1 − 2Rm γj Rj

c .
Rj (85)

The coefficients 0.4 and 2 are very approximate (to emphasize

this, we give their values to only one significant digit). Nev-

ertheless, as we show in Figures 9 and 10, this approximate

formula does quite a good job of fitting the numerical results

for a wide range of models. The only region of parameter space

where the formula fails is when the underlying equilibrium be-

comes “nonrelativistic” and γj approaches unity. Such models,

which are located near the upper end of Figures 9 and 10, have

extremely large growth rates.

Although modes with Rsing ∼ Rj have the largest growth rates, these modes have relatively short wavelengths Rj along

the z-axis (see Equation (84)). With such short wavelengths it

is not clear if the instability can grow to a large amplitude.

It is therefore interesting to consider modes with k → 0.

Figure 4 shows that Model B is stable as k → 0, whereas

Figure 6 indicates that the k = 0 mode in Model C is nearly as

unstable as the fastest growing mode.

The dashed lines in Figures 9 and 10 show numerical

results for the growth rates of modes with k = 0 for various

combinations of the model parameters γm and Rm. For small

No. 2, 2009

STABILITY OF RELATIVISTIC FORCE-FREE JETS

1691

Figure 10. Numerically calculated growth rates of modes with m = 1 and Rsing = 1.1Rj . The solid lines show the results for a series of models with fixed Rm/Rj and varying γm. The dotted lines are the growth rates predicted by Equation (85). The dashed lines are the numerical growth rates for modes with k = 0.

values of Rm 0.1Rj , the results are nearly identical to those we described above for the fastest growing mode (Rsing =
1.1Rj , solid lines). This is to be expected based on the results shown in Figure 6 for Model C, which has Rm = 0.18Rj . With increasing Rm, however, the k = 0 modes become less unstable than the modes with Rsing ∼ Rj . By Rm ∼ 0.3Rj , the k = 0
modes are fully stable, thus explaining the stability of Model B
as k → 0 (Figure 4).

4.5. Spatial Growth of Unstable Modes

The discussion so far was limited to modes with real k and complex ω. An equally interesting problem is to consider modes with real ω and complex k. From Equation (41), we see that the eigenfunctions take the form

E1 = [E1R(R)Rˆ + E1φ(R)φˆ + E1z(R)zˆ] × exp[−iωt + imφ + iRe(k)z − Im(k)z]
∝ exp[iRe(k)z] exp(z/Z),

(86)

where ω is real and Z = −1/Im(k) is the scale length on which

the mode e-folds in the z-direction. Such spatially growing

“convected” modes are particularly relevant for sources with

long-lived steady jets.

As discussed in Payne & Cohn (1985) and Appl & Camenzind

(1992), there is a strong symmetry between modes with real

k and complex ω, and those with complex k and real ω. In

particular, the growth rates of the two kinds of modes are related

by

Im(k) = −Im(ω)/vg,

(87)

where vg = ∂Re(ω)/∂Re(k) is the group velocity of the mode. Our unstable m = 1 modes have vg very nearly equal to c. Therefore, we immediately obtain from Equation (85) the
following estimate for the spatial e-folding scale Z of the fastest
growing convected mode:

Mode with Rsing = 1.1Rj :

Z Rj

≈

2.5γj . (1 − 2Rm/Rj )

(88)

Figure 11. Numerically calculated e-folding scales Z for modes with m = 1, ω real and Rsing = 1.1Rj . These modes have among the largest growth rates. The solid lines show results for a series of models with fixed γm and varying Rm/Rj . The dotted lines are the growths predicted by Equation (88). Note the very good agreement except near the bottom of the plot, where the models are
nonrelativistic. The dashed lines are the numerical values of Z for modes with ω = 0.

Zero-frequency modes (ω = 0) should have almost the same
Z for small values of Rm/Rj , but the growth should cut off at a somewhat smaller value of Rm compared to the modes with Rsing = 1.1Rj .
Figure 11 shows numerical results. Modes with Rsing = 1.1Rj have growths consistent with Equation (88), and the modes with ω = 0 have similar growths except that the instability cuts off
at somewhat smaller values of Rm. The results are as expected and are very similar to those shown in Figure 9.
The above results correspond to an idealized cylindrical jet.
In the case of real jets we must allow for a finite opening angle θj ≡ dRj /dz. Many force-free jet models have opening angles that vary inversely as the Lorentz factor: θj ∼ few/γj (TMN08). Using this scaling we can estimate approximately the evolution
of the mode amplitude a with distance:

da = 1 da ≈ γj da ≈ γj a ≈ γj

a

dRj θj dz few dz few Z few 2.5γj Rj

≈1

a .

few × 2.5 Rj

(89)

Solving this differential equation, and using Rj ∝ zα/4 ∼ z0.5−1 (Equation (7)), we obtain

a(z) ∝ z ,

(90)

where 0.1. This estimate is very crude, but it does suggest that, in realistic jets, the unstable kink mode we have studied in this paper grows only weakly with increasing distance.

4.6. Toward an Improved Instability Criterion
In Section 1, we introduced three different instability criteria, of which the IPL and TMT criteria refer specifically to rotating force-free jets. Since Bz is practically constant in our equilibria, all of our models are close to the boundary between stability

1692

NARAYAN, LI, & TCHEKHOVSKOY

Vol. 697

and instability according to the IPL criterion (Equation (2)). Similarly, since Bφ2 ≈ Ω2R2Bp2/c2 in our equilibria, our models are marginally stable according to the TMT criterion (Equation (3)).13 It is thus not possible to understand from either
of these criteria why Model A is stable and Models B and C are
unstable.
It is, of course, not surprising that the above instability criteria
fail. Our jet equilibria include the effects of poloidal field
curvature, which were not considered by the previous authors.
For easier comparison with previous work, let us rewrite our
balance condition (23) as follows:

1d R2 dR

(B02φ − E02R)R2 8π

=− d dR

B02z 8π

+

E02R − B02z 4π Rc

.

(91)

If we leave out the last term, this is equivalent to Equation (6) in

Istomin & Pariev (1996) and Equation (16) in Lyubarskii (1999).

The IPL instability criterion states that the quantity dB02z/dR should be negative. We might wish to generalize this by saying

that the right-hand side of Equation (91), including the poloidal

curvature term, should be positive. Unfortunately, this simple

modification is not sufficient since the right-hand side is positive

for all of our models, whereas not all our models are unstable;

not only should the right-hand side be positive, its magnitude

should be larger than some amount. The same seems to be true

with the TMT criterion. This criterion indicates that all our

equlibria should be unstable, whereas only some of them are.

Qualitatively, what distinguishes the unstable Models B and

C from the stable Model A is that the former have made the

transition to the second acceleration regime, where poloidal field

curvature has a strong effect on transverse force balance. This is

reflected in the γ (R) profiles (Figure 2) which have dγ /dR < 0

at larger radii. Thus, one might guess that instability requires

the jet to be in the second acceleration regime, i.e., poloidal

field curvature to be important, or the jet to have a declining

γ (R). Once again, these related conditions by themselves are

not sufficient. To have an instability, γ (R) should decline over

a sufficiently broad range of radius, e.g., Rm should be less than ∼ 0.45Rj in our models.

An alternate approach which we have found useful is to focus

on the left-hand side of Equation (91). From Equations (17)–

(19), we see that for our equilibria we have B02φ − E02R = (R/Rm)4B02z. Furthermore, we have seen that modes with Rsing ∼ Rj (the fastest growing modes) are unstable so long as Rm 0.45Rj , while modes with k → 0 (long-wavelength
modes) are unstable for Rm 0.3Rj . From this, we obtain the
following approximate instability criteria:

Modes with Rsing ∼ Rj: B02φ − E02R 1/2 > 5|B0z|, (92)

Modes with k → 0: B02φ − E02R 1/2 > 12|B0z|, (93)
where we have set R = Rj to obtain the numerical coefficients on the right. These conditions are easier to interpret if we boost to the comoving frame of the fluid (V. Pariev 2008, private communication), where the electric field vanishes. In this frame, the criterion for instability becomes

Modes with Rsing ∼ Rj: |B0φ,comov| > 5|B0z,comov|, (94)

13 A strict application of the TMT criterion would indicate that our models are unstable, since Bφ2 > E2 = Ω2R2Bp2 /c2. However, Bφ2 − E2 Bφ2, so the models deviate only slightly from marginal stability.

Modes with k → 0: |B0φ,comov| > 12|B0z,comov|. (95)
That is, in the comoving frame, the toroidal field must dominate the poloidal field by more than a certain critical factor.14 Written in this form, the condition resembles the KS criterion (Equation (1)).
5. SUMMARY AND DISCUSSION
The relativistic jet model we have considered in this paper is particularly simple: it is cylindrical, it assumes force-free conditions, and it assumes rigid rotation. Within the limitations of these reasonable approximations, we have attempted to be as close to numerically simulated jets as possible. A unique feature of our model is that we include the effect of poloidal field curvature, which is known to play an important role in numerical force-free jets (Section 2.1). Also, we choose functional forms for the various field components in the equilibrium (Section 2.2) to match as closely as possible our previous force-free simulations (TMN08).
Our equilibrium model is described by two parameters: the maximum Lorentz factor, γm, and the radius at which this maximum is achieved, Rm. The ratio of the latter to the jet radius, Rj, determines the basic physics of the equilibrium. Models in which Rm/Rj > 1 have γ (R) increasing monotonically with radius R out to some maximum Lorentz factor γj < γm at the outer edge of the jet. Model A (Figure 2) is an example. In these models the entire jet is in the first acceleration regime (see Sections 2.1, 2.2.1 and TMN08 for details). We find that all these models are perfectly stable.
Models with Rm/Rj < 1 are more interesting. Here, γ (R) increases up to a maximum value γm at R = Rm and then decreases down to a Lorentz factor γj < γm at R = Rj . Models B and C (Figure 2) are examples. In these models, the jet fluid at R < Rm is in the first acceleration regime, while the fluid at Rm < R < Rj is in the second acceleration regime. We find that the subset of these models with Rm/Rj 0.45 are linearly unstable. For Rm/Rj just below 0.45, all the unstable modes have short wavelengths in the z-direction: λ = 2π/kz ∼ 2π Rj /γj . With decreasing Rm, a wider range of kz becomes unstable, and for Rm/Rj 0.3, waves with kz = 0, i.e., with arbitrarily long wavelengths, become unstable. The latter modes are perhaps of most interest since they are likely to grow to the largest amplitudes. The numerical results are summarized in Figures 3–11.
The unstable modes we find are all kink modes with azimuthal wavenumber m = ±1. These are nonaxisymmetric modes in which the jet is distorted helically. A key feature is that, at each z, the center of mass of the jet is shifted away from R = 0. It is well known that MHD configurations with toroidal fields are especially susceptible to the kink mode (Bateman 1978), and our models follow this trend. However, because our equilibria both rotate and move relativistically along z, the criterion for instability is different from the usual KS criterion (Equation (1)).
14 Since our equilibria assume a constant Ω for all field lines, the criteria (94) and (95) are technically valid only for such models. However, since the criteria have been written without any explicit reference to Ω, they may be valid more generally even when Ω varies with R. Also, in the case of a spinning black hole, the most important example of a central “star” with nonconstant Ω, the angular velocity of the horizon varies by only a factor of 2 between the pole and the equator (Blandford & Znajek 1977; McKinney & Narayan 2007a, 2007b). Such a modest variation of Ω across the jet is unlikely to modify our results substantially, though this remains to be demonstrated.

No. 2, 2009

STABILITY OF RELATIVISTIC FORCE-FREE JETS

1693

We find that the typical growth rate of the unstable kink mode in our jet models is given by Equation (85): the e-folding time is of order γj times the light-crossing time Rj /c across the jet. The extra factor of γj is easily understood—it arises from time dilation. In the comoving frame of the jet (at R = Rj ), we expect the growth timescale to be simply Rj /c. Lorentz transforming to the lab frame, this becomes γj Rj /c. For convected modes with a real frequency, the e-folding length scale is of order γj times the jet radius Rj, the extra γj arising in this case from Lorentz contraction. Since jets typically have opening angles ∼ few/γj , the net result is that the unstable modes grow only slowly with distance from the base of the jet (Equation (90)). Of course, relativistic jets in astrophysical sources propagate over many decades, so in principle even this slow growth might lead to a noticeable perturbation amplitude. Nevertheless, the fact that the growth is very slow reduces the seriousness of the kink instability.
Our jet equilibria turn out to be close to the boundary between stability and instability according to either the IPL or the TMT criterion (Equations (2) and (3)), so these criteria are not useful for interpreting the results. In addition, since our models include the effects of poloidal field curvature, they lie outside the range of validity of the IPL and TMT criteria. The most useful instability criterion we have come up with is that, in the comoving frame of the jet fluid, the tangential field should be an order of magnitude or more larger than the poloidal field (Equations (94) and (95)). Expressed thus, the criterion is similar to the KS criterion (1). However, we should not apply the KS criterion in the lab frame. Rather, we should apply it in the comoving frame of the jet fluid, which is reasonable, and we should set z ∼ Rj , which is again reasonable because of Lorentz contraction (z ∼ γj Rj in the lab frame implies z ∼ Rj in the comoving frame).
All the work described here assumes a rigid wall enclosing the jet at the boundary R = Rj . We have done some calculations with a constant pressure boundary and we find unstable modes with much larger growth rates compared to the rigid wall case. However, since we are dealing with a force-free jet, it is not clear that a constant pressure boundary is particularly meaningful. For instance, if the pressure is from a nonrelativistic gaseous envelope or cocoon, the gas would have substantial inertia and would probably behave to first approximation like a rigid wall.
Various authors have discussed mechanisms by which instabilities might be suppressed in astrophysical jets. Hardee et al. (2007, and references therein) have shown that an external wind or cocoon can stabilize the Kelvin–Helmholtz mode in MHD jets, though it is not clear if this is relevant for force-free jets. Moll et al. (2008) show that lateral expansion causes instabilities to grow more slowly. In a sense, we have already included this effect when we derived the growth-rate estimate given in Equation (90). In addition, we note that some of the growth suppression seen by Moll et al. is probably because expansion causes different parts of the jet to lose causal contact with one other. This is not an issue for force-free models, where signals propagate at the speed of light.
It would be interesting to simulate numerically the unstable modes described in this paper. Apart from verifying the linear theory, such calculations will reveal the nonlinear development of the mode. Does the kink mode saturate at a finite amplitude and lead to a more-or-less coherent helical pattern or does it destroy the initial equilibrium? This important question can be answered only with three-dimensional simulations. Since the kink mode involves lateral motion of the jet across the

axis R = 0, the numerical technique used must be flexible enough to allow such motions (e.g., as described by McKinney & Blandford 2009).
We conclude by reminding the reader that the work described here refers to a particularly simple model of relativistic jets which is based on the force-free approximation. In real jets, once the flow crosses the fast magnetosonic point, the inertia of the gas starts to play a role and the force-free approximation is no longer valid (e.g., Tchekhovskoy et al. 2009; Lyubarsky 2009). We must then consider the full MHD equations.

The authors thank Alison Farmer for assistance during the early stages of this work and Jonathan McKinney for numerous helpful discussions and comments on the paper. This work was supported in part by NASA grant NNX08AH32G.

APPENDIX A

CONSTANCY OF POLOIDAL MAGNETIC FIELD ACROSS FORCE-FREE JETS

Figure 1 shows that in numerical simulations of force-free jets

Bp hardly changes with R. Here we show that this is a common

feature of all jet solutions that smoothly connect to a spinning

compact object at the base.

Consider the force balance equation, Equation (8). Suffi-

ciently near the compact object, where the jet is in the first

acceleration regime (see Section 2.1), we can drop the terms proportional to Rc−1 and (Bφ2 − E2) since in the first acceleration regime γ12 γ22 leading to E2/Bp2 Rc/R and Bφ2 −E2 Bp2 (TMN08). Then the force balance equation, Equation (8), sim-

plifies to

d Bp2 ≈ 0. dR

(A1)

Therefore, sufficiently near the compact object the poloidal field

is nearly constant,

Bp(R) ≈ const.

(A2)

Each field line is labeled by the amount of poloidal magnetic

flux Φ that it encloses. Due to Equation (A2) this flux can be

written simply as

Φ ≈ π BpR2.

(A3)

These relations are valid throughout the first acceleration regime. We now show that they actually hold asymptotically in all parts of the jet.
Recall that in force-free magnetospheres the enclosed poloidal current I is preserved along each field line (Mestel 1961; Okamoto 1978; Thorne et al. 1986; Beskin 1997; Narayan et al. 2007; Tchekhovskoy et al. 2008),

I

=

I (Φ)

=

c 2

RBφ

≈

−Ω 2

Bp R 2 ,

(A4)

where the approximate equality is due to Bφ ≈ −E =

−ΩRBp/c for R

RA (see Equations (10) and (12)).

Comparing Equations (A3) and (A4) and recalling that Ω is

conserved along field lines, we obtain

Φ ≈ − 2π I (Φ). Ω

(A5)

Since both sides of this relation depend only on Φ, this relation is valid everywhere in the solution, even though we derived it only in the first acceleration regime. Using Equation (A4) to

1694

NARAYAN, LI, & TCHEKHOVSKOY

Vol. 697

substitute for I yields back Equation (A3). Thus Equations (A2) and (A3) are valid everywhere in the jet.
APPENDIX B COEFFICIENTS IN EQUATIONS (56) AND (57)
Here, we give explicit expressions for the coefficients D1(R)– D6(R) defined in Section 3.1. The functions are

ickm ick2 2ic

C6(R) =

f− ωR

ω g − ωR2 g

+

iω g

−

im h,

cR

C7(R)

=

2ic f
ωR

+

ic f
ω

2ic + C5C1 + ω gC1,

(B12) (B13)

D1(R) =

−

ck g−
ω

cm ωR gC1,

(B1)

D2(R) =

−

cm ωR2 f

−

cm f
ωR

+

C2C1

−

cm ωR gC1,

(B2)

D3(R) =

ick2 ω

+

icm2 ωR2

−

iω c

g,

(B3)

D4(R)

=

ck f
ω

+

cm g
ωR

−

h

−

C3

D3 D1

,

(B4)

D5(R)

=

C4

+

C3

D22 D12

+

C6

+

C5 R

C1

+

C5C1

−

C7

D2 D1

+ C3

D2D1 D12

−

D2 D1

ic + ω gC1 ,

(B5)

D6(R)

=

C8

+

C3

D2D3 D12

−

C7

D3 D1

+ C3

D3D1 D12

−

D3 D1

,

(B6)

where primes denote derivatives with respect to R, and the functions C1(R)–C8(R) are given by

ωRg − cmh

C1(R)

=

ωRf

−

. ckRh

(B7)

ck cm h h ck

C2(R)

=

f ωR

−

ωR2 g

−

R

−

Rc

+

f ω

−h,

(B8)

ic

C3(R) =

(f ω

+ gC1),

(B9)

icm2 iω ickm

C4(R) = − ωR2 f +

f+ c

ωR

g − ikh,

(B10)

ic ic ic

C5(R)

=

g ωR

+

g ωRc

+

g ω

,

(B11)

C8(R)

=

2ck f
ωR

−

cm g
ωR2

−

2h R

+

cm g
ωRcR

−

h Rc

+

ck f
ω

+

cm g

−h.

(B14)

ωR

REFERENCES

Appl, S., & Camenzind, M. 1992, A&A, 256, 354 Appl, S., Lery, T., & Baty, H. 2000, A&A, 355, 818 Bateman, G. 1978, MHD Instabilities (Cambridge, MA: MIT Press) Begelman, M. C. 1998, ApJ, 493, 291 Benford, G. 1981, ApJ, 247, 792 Beskin, V. S. 1997, Sov. Phys. Usp., 40, 659 Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 Ferrari, A., Trussoni, E., & Zaninetti, L. 1978, A&A, 64, 43 Gammie, C. F., McKinney, J. C., & To´th, G. 2003, ApJ, 589, 444 Goldreich, P., Goodman, J., & Narayan, R. 1986, MNRAS, 221, 339 Goldreich, P., & Julian, W. H. 1969, ApJ, 157, 869 Goldreich, P., & Tremaine, S. 1978, ApJ, 222, 850 Hardee, P., Mizuno, Y., & Nishikawa, K. I. 2007, Ap&SS, 311, 281 Istomin, Y. N., & Pariev, V. I. 1996, MNRAS, 281, 1 Kadomtsev, B. B. 1966, Rev. Plasma Phys., 2, 153 Lery, T., Baty, H., & Appl, S. 2000, A&A, 355, 1201 Li, L. X. 2000, ApJ, 531, L111 Lynden Bell, D. 2006, MNRAS, 369, 1167 Lyubarskii, Y. E. 1999, MNRAS, 308, 1006 Lyubarsky, Y. 2009, ApJ, submitted (arXiv:0902.3357) McKinney, J. C. 2006, MNRAS, 368, 1561 McKinney, J. C., & Blandford, R. D. 2009, MNRAS, 394, L126 McKinney, J. C., & Gammie, C. F. 2004, ApJ, 611, 977 McKinney, J. C., & Narayan, R. 2007a, MNRAS, 375, 513 McKinney, J. C., & Narayan, R. 2007b, MNRAS, 375, 531 Mestel, L. 1961, MNRAS, 122, 473 Mignone, A., & McKinney, J. C. 2007, MNRAS, 378, 1118 Mizuno, Y., Hardee, P., & Nishikawa, K. I. 2007, ApJ, 662, 835 Moll, R., Spruit, H. C., & Obergaulinger, M. 2008, A&A, 492, 621 Narayan, R., Goldreich, P., & Goodman, J. 1987, MNRAS, 228, 1 Narayan, R., McKinney, J. C., & Farmer, A. J. 2007, MNRAS, 375, 548 Okamoto, I. 1978, MNRAS, 185, 69 Payne, D. G., & Cohn, H. 1985, ApJ, 291, 655 Ruderman, M. A., & Sutherland, P. G. 1975, ApJ, 196, 51 Tchekhovskoy, A., McKinney, J. C., & Narayan, R. 2007, MNRAS, 379, 469 Tchekhovskoy, A., McKinney, J. C., & Narayan, R. 2008, MNRAS, 388, 551 Tchekhovskoy, A., Mc Kinney, J. C., & Narayan, R. 2009, ApJ, submitted
(arXiv:0901.4776) Thorne, K. S., Price, R. H., & Mac Donald, D. A. 1986, Black Holes: The
Membrane Paradigm (New Haven, CT: Yale Univ. Press) Tomimatsu, A., Matsuoka, T., & Takahashi, M. 2001, Phys. Rev. D, 64, 123003