SUPERLUMINAL WAVES IN PULSAR WINDS

The circularly polarized wave is described in the H-frame by three (phase-independent) parameters: the (constant) proper number density of each species n = n₀, the component of the four-velocity in the direction of propagation u_x = u_x0, and the magnitude of the four-velocity component perpendicular to the propagation direction $\sqrt{\vphantom{A^A}\smash{\hbox{${u_y^2+u_z^2}$}}}$, which, for this wave mode, equals the dimensionless strength parameter denoted by a and defined in terms of the electric field amplitude E₀ as

$$\begin{eqnarray} a&=&e|E_0|/(mc\omega). \end{eqnarray} \tag{ 1 }$$

Consequently, the Lorentz factor of the fluids is also phase independent: $\gamma =\gamma _0=\sqrt{\vphantom{A^A}\smash{\hbox{${1+u_{x0}^2+a^2}$}}}$. The frequency ω coincides with the (proper) plasma frequency, ω_p:

$$\begin{eqnarray} \omega &=&\omega _{\rm p}\,=\, \sqrt{\frac{8\pi n_0 e^2}{m}}. \end{eqnarray} \tag{ 2 }$$

In terms of these quantities, the particle flux density J and the (0, 0), (0, x), and (x, x) components of the stress energy tensor are given by (see Appendix B, Equations (A37)–(A39) and (A40)–(A42))

$$\begin{eqnarray} J &=& n_0\gamma _0 c \left(\frac{2 u_{x0}}{\gamma _0}\right) \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} T^{00} &=& n_0\gamma _0 mc^2 \left(2\gamma _0 + \frac{a^2}{\gamma _0}\right) \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} T^{01} &=& n_0\gamma _0 mc^2 \left(2u_{x0}\right) \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} T^{11} &=& n_0\gamma _0 mc^2 \left(\frac{2u_{x0}^2}{\gamma _0} + \frac{a^2}{\gamma _0}\right) \end{eqnarray} \tag{ 6 }$$

and are also phase independent.

The linearly polarized mode requires an additional parameter to fix the amplitude of the phase-averaged magnetic field. In the H-frame, the four parameters are the proper number density n₀ and the component u_x0 of the four-velocity in the direction of propagation, both taken at phase zero (the point in the wave where the magnitude of the electric field is maximal and u_z = u_y = 0), the nonlinearity parameter q, used by Kennel & Pellat (1976), and the ratio λ of the (phase-independent) magnetic field to the electric field E₀ at phase zero. Physically q is double the ratio at phase zero of the energy density in the fluids compared with that in the electric field, as measured in the H-frame:

$$\begin{eqnarray} q&=&\frac{32\pi n_0\gamma _0^2mc^2}{E_0^2}. \end{eqnarray} \tag{ 7 }$$

Defining the normalized field variable y = E/E₀, one finds

$$\begin{eqnarray} u_x&=&u_{x0}+4\gamma _0(1-y)\lambda /q \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} \gamma &=&\gamma _0+2\gamma _0(1-y^2)/q \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} &&| u_y |=\sqrt{\gamma ^2-u_x^2-1}=2\gamma _0\sqrt{N(y)}/q \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} u_z &=& 0, \end{eqnarray} \tag{ 11 }$$

where N(y) is a fourth-order polynomial, as follows from Equations (8) and (9), and $\gamma _0=\sqrt{\vphantom{A^A}\smash{\hbox{${1+u_{x0}^2}$}}}$. Periodic wave solutions exist provided N(y) has four distinct real roots y_1...4, in which case y oscillates between the values y₂ and y₃ (assuming the ordering y₁ < y₂ < y₃ < y₄). Because of the normalization of y, one has y₃ = 1 and |y₂| ⩽ 1. The strength parameter of the wave, as defined in Equation (1), is

$$\begin{eqnarray} a&=&\frac{2}{\pi } \int _{y_2}^1\,dy\,\gamma /\sqrt{N(y)} \end{eqnarray} \tag{ 12 }$$

and the wave frequency is given by

$$\begin{eqnarray} \omega &=&\frac{2\gamma _0}{a\sqrt{q}}\omega _{\rm p}. \end{eqnarray} \tag{ 13 }$$

The phase-averaged fluxes are

$$\begin{eqnarray} \langle J\rangle &=& n_0\gamma _0 c \left<\frac{2 u_x}{\gamma }\right> \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} \langle T^{00}\rangle &=& n_0\gamma _0 mc^2 \left<2\gamma + \frac{4\gamma _0(y^2+\lambda ^2) }{q}\right> \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} \langle T^{01}\rangle &=& n_0\gamma _0 mc^2 \left<2u_x + \frac{8\gamma _0y \lambda }{ q}\right> \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} \langle T^{11}\rangle &=& n_0\gamma _0 mc^2 \left<\frac{2u_x^2}{\gamma } + \frac{4\gamma _0(y^2+\lambda ^2)}{q}\right> \end{eqnarray} \tag{ 17 }$$

and the phase-averaged electric field is

$$\begin{eqnarray} \frac{\langle E\rangle^2}{4\pi }&=&n_0\gamma _0mc^2 \left(\frac{4\gamma _0 \langle y\rangle^2}{q}\right). \end{eqnarray} \tag{ 18 }$$

Ampere's law enables these averages to be expressed as integrals over the normalized electric field y:

$$\begin{eqnarray} \langle A(y)\rangle &=& \frac{\int _{y_2}^1\,dy\,A(y)\gamma /\sqrt{N(y)}}{ \int _{y_2}^1\,dy\,\gamma /\sqrt{N(y)}} \end{eqnarray} \tag{ 19 }$$

(see Appendix A).

If the phase-averaged magnetic field vanishes (λ = 0) one finds y₂ = −1, and the integrals in Equations (13)–(18) can be expressed in closed form in terms of elliptic integrals. For λ ≠ 0, closed forms can be found to lowest order in an expansion in the small parameter q (see Kennel & Pellat 1976). However, these are not adequate to describe the pulsar case, as we show in the following section. In the general case, one must resort to numerical integration, noting that the integrands have integrable singularities at each of the limits.

In order to identify those parts of the wind of a given pulsar in which superluminal waves can propagate, one has to associate the wave properties with wind quantities inferred from the observational data. A pulsar wind transports particles, energy, and magnetic flux and exerts a ram-pressure on its surroundings. Across a stationary shock front the quantities J', T'⁰¹, and T'¹¹ are conserved. Also, from Faraday's law, the transverse electric field component E' is conserved. Across a thin transition front where one wave-mode converts into another, and which is stationary on the oscillation timescale, the phase-averaged values of these quantities are conserved. "Thin" in this connection means that the extension of the front in the propagation (radial or the x-direction) is small compared with its transverse (y and z) dimensions, but not necessarily compared to a wavelength. This is guaranteed in the plane-wave approximation. On the other hand, "stationary" means that the location of the front does not change on timescales that are long compared with the oscillation timescale. The transverse electric field is proportional to the rate at which magnetic flux enters and leaves the front. This is easily seen in the case of the subluminal wave, where E' = β_<B'. In the linearly polarized superluminal wave, it follows from the phase-averaged equation of motion for the transverse momentum (A6) written in the laboratory frame: 〈E'〉 = 〈u'_xB'/γ'〉.

However, it is the relative magnitudes of the fluxes, rather than their absolute values, that determine the physics of wave propagation. We therefore introduce three dimensionless quantities:

$$\begin{eqnarray} &&\mu = \frac{\langle T^{\prime 01} \rangle}{mc \langle J^{\prime } \rangle} \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} &&\nu = \frac{\langle T^{\prime 11} \rangle}{mc \langle J^{\prime } \rangle } \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} &&\eta = \frac{\langle E^{\prime } \rangle ^2/4\pi }{mc \langle J^{\prime } \rangle }, \end{eqnarray} \tag{ 22 }$$

which are conserved across a transition layer and are independent of the density parameter n₀. μ is the mass-loading parameter introduced by Michel (1969), which corresponds to the Lorentz factor each particle would need if they were to carry the entire energy flux.

For the superluminal wave modes discussed above, these quantities can be computed by Lorentz transforming their values in the H-frame into the pulsar or laboratory frame. A boost of (dimensionless) speed β_> = c/v_ϕ in the negative x-direction is required, for which explicit expressions are given in Appendix B. (We use the notation β_> and Γ_> for the speed and associated Lorentz factor of the H-frame of the superluminal wave, and β_< and Γ_< for the phase (or bulk) speed and associated Lorentz factor of the subluminal wave, both as seen from the laboratory (pulsar) frame.) Applying these boosts, one obtains, for linearly polarized modes, expressions for μ, ν, and η, as functions of the four input parameters u_x0, q, λ, and β_> that must be evaluated by numerical integration. The case of circular polarization, where η = 0, is much simpler, since μ and ν are given in closed form as functions of u_x0, a, and β_>.

For the cold subluminal wave modes, such as the striped wind (in the absence of dissipation), μ and ν are independent of radius in the inner parts of the wind where a single-fluid MHD description is adequate (Kirk & Mochol 2011a). Their values can be estimated from standard pulsar models (e.g., Lyubarsky & Kirk 2001). However, these two quantities are almost equal in highly relativistic winds, and it proves more convenient to use, instead of ν, the magnetization parameter σ, defined as the ratio of magnetic to particle energy flux according to

$$\begin{eqnarray} \mu &=& \Gamma _< (1+\sigma) \end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray} \nu &=& \frac{\Gamma _<^2 (1+\sigma) - (1+\sigma /2)}{\sqrt{\Gamma _<^2-1}}. \end{eqnarray} \tag{ 24 }$$

In the striped wind, the third parameter, η, depends on the phase-averaged magnetic field value

$$\begin{eqnarray} b&=&\langle B^{\prime }\rangle/\langle B^{\prime 2}\rangle^{1/2}, \end{eqnarray} \tag{ 25 }$$

via the expression:

$$\begin{eqnarray} \eta &=& b^2\sigma \beta _<\Gamma _<. \end{eqnarray} \tag{ 26 }$$

Without loss of generality, we may choose 〈B'〉 ⩾ 0, so that 0 ⩽ b ⩽ 1. The above expression has been calculated taking into account that E' = β_<B' in the striped wind. The corresponding quantity for the superluminal wave follows from the expressions given in Appendix B.

Both b and η are functions of latitude in the striped pulsar wind, since the spacing in phase of the current sheets varies. At latitudes greater than the inclination angle between the magnetic and rotation axes, conventionally denoted by α, the sheets vanish (zero spacing), whereas on the equator, they are equally spaced (separated in phase by π).

Thus, for given μ, σ, and b, the jump conditions can be solved for three of the four input parameters u_x0, q, b, and β_>, provided the fourth is specified. The wave frequency, normalized to ω_p, can then be constructed from Equation (13).

For circular polarization, the jump conditions can be solved analytically (Kirk 2010). An example is plotted in Figure 1, which shows the refractive index ck'/ω' = β_>, the strength parameter a, and the component of the particle speed in the propagation direction, u_x0, as a function of laboratory frame frequency ω' normalized to the (proper) plasma frequency ω_p. In this example, we have chosen moderate values of the energy and momentum fluxes in the laboratory frame: μ = 12, σ = 3 (and, therefore, Γ_< = 3), in order to display the structure of the solutions near the frequency cutoff. Since the wavenumber k vanishes in the H-frame, the refractive index is trivially obtained from the Lorentz transformation of the wave number and frequency: ω' = Γ_>ω and k' = β_>Γ_>ω/c, (this also applies, of course, to the case of linear polarization) and the dispersion relation ω = ω_p implies ω'/ω_p = Γ_>. However, the corresponding strength parameter a is found only after solving the full set of equations. As can be seen from this figure, there are no physical solutions for ω' < ω₀, two solutions for ω₀ < ω' < ω₁, and only one physical solution for ω' > ω₁. (The full expressions for ω_{0, 1} are cumbersome, but for μ ≳ σ^3/2 ≫ 1 (corresponding to a mildly supermagnetosonic outflow), one finds ω₀ ≈ ω_p[1 + (4μ² − σ³)σ/(32μ⁴) + O(μ^−5/3)] and ω₁ ≈ ω_p(μ/8)^1/4[1 + O(μ^−1/3)], (cf. Kirk 2010).) Thus, the imposition of a finite energy flux per particle prohibits wave propagation at frequencies close to the cutoff ω' = ω_p found from linear theory. At intermediate frequencies, waves of two different amplitudes, and different values of u_x0, are available to transport the required fluxes. One of these disappears (a → 0) at finite ω'/ω_p, leaving only a single solution of the jump conditions at high frequency.

Zoom In Zoom Out Reset image size

In the case of a pulsar wind, the wave frequency is dictated by the rotation of the neutron star, but the plasma density decreases with increasing distance from the star. Thus, rather than solving the jump conditions as a function of frequency, more insight can be gained by solving as a function of radius. The connection between flux density and radius is found in terms of the luminosity per unit solid angle of the (radial) wind, dL/dΩ_s:

$$\begin{eqnarray} \langle F^{\prime }\rangle &=& \frac{1}{r^2}\frac{dL}{d\Omega _{\rm s}}, \end{eqnarray} \tag{ 27 }$$

which can be expressed in dimensionless form in terms of the quantity

$$\begin{eqnarray} a_{\rm L}&=&\sqrt{\frac{4\pi e^2}{m^2c^5}\frac{dL}{d\Omega _{\rm s}}}. \end{eqnarray} \tag{ 28 }$$

Physically, a_L is the strength parameter of the circularly polarized vacuum wave that would be needed to carry the entire pulsar luminosity at the surface ρ = 1. For a spherically symmetric wind a_L = 3.4 × 10¹⁰L^1/2₃₈, where L₃₈ = L/(10³⁸ erg s⁻¹). This quantity is related to the multiplicity parameter κ, defined as the Goldreich–Julian charge density at the light cylinder divided by e (see, for example, Lyubarsky & Kirk 2001; Kirk & Mochol 2011a) by

$$\begin{eqnarray} \kappa &=& a_{\rm L}/(4\mu). \end{eqnarray} \tag{ 29 }$$

Writing the energy flux density as μ〈J'〉 leads to the expression

$$\begin{eqnarray} \langle J^{\prime }\rangle &=& \frac{2 n_0c\Gamma _>^2\omega ^2 a_{\rm L}^2}{\rho ^2\omega _{\rm p}^2\mu }. \end{eqnarray} \tag{ 30 }$$

Applying this to the case of circular polarization, where 〈J'〉 = 2n₀u'_x0 and ω = ω_p, one finds u'_x0 = Γ² _>a²_L/(μρ²). Then, the inequalities u'_x0 < γ' ⩽ μ and Γ_> > 1 imply that circularly polarized waves propagate only when ρ > a_L/μ. We therefore introduce the dimensionless scaled radius

$$\begin{eqnarray} &&R=\rho \mu /a_{\rm L} = \rho /(4\kappa), \end{eqnarray} \tag{ 31 }$$

which can be constructed after solving the jump conditions from the expression

$$\begin{eqnarray} R&=& \frac{\Gamma _>\omega }{\omega _{\rm p}} \sqrt{\frac{2\mu }{\gamma _0(\langle J^{\prime }\rangle /n_0\gamma _0)}} \end{eqnarray} \tag{ 32 }$$

for both circularly and linearly polarized modes.

Figure 2 shows the properties of circularly polarized waves as functions of R for parameters appropriate to pulsar winds: σ = 100 and μ = 10100 (and, therefore, Γ_< = 100). Near the cutoff radius, which lies close to R = 1, the mode properties are similar to those seen in Figure 1, except that the refractive index, now plotted as a function of radius, is no longer degenerate. The detailed properties of these solutions have been discussed in Kirk (2010). However, note that these curves do not describe the radial evolution of a wave packet, as erroneously suggested in that paper and assumed by Asseo et al. (1984, who considered linear polarization), but simply specify the wave modes into which a wind of given μ and σ can convert at a given radius. The difference arises because the radial momentum flux density, ν, is not a conserved quantity in the radial evolution of the superluminal modes (see also Kirk & Mochol 2011b).

Zoom In Zoom Out Reset image size

When the phase-averaged magnetic field vanishes, the properties of linearly polarized modes are similar to those of circularly polarized modes. This is illustrated in Figure 3, where one sees that the refractive indices are practically identical.

Zoom In Zoom Out Reset image size

This figure also shows the nonlinearity parameter q. For linear polarization, q is defined in terms of quantities measured at phase zero. This point in the wave is, however, special. It corresponds to the turning points of the "saw-tooth" waveform (Max & Perkins 1971), where the fluid velocity either vanishes or lies precisely in the propagation direction. These points do not play a large role in determining the average properties of the wave. As a result, q varies substantially, especially close to points where u_x0 changes sign (log (R) ≈ 0.4 in the figure). One consequence is that an expansion of the fluxes using q as a small parameter is inadequate over most of the range relevant for pulsars, despite the fact that the waves are highly nonlinear (see the discussion in Section 4). In the case of circular polarization, the quantities that enter into the definition of q (Equation (7)) are phase independent. The value of q is a useful characterization of the mode in this case, and Figure 3 clearly illustrates that one of the two solutions for given R is a relatively weak wave, with q ∼ few, whereas the other is stronger, with q ∼ 10⁻².

The properties of linearly polarized waves that carry a non-zero phase-averaged magnetic field in the toroidal direction are shown in Figure 4. The refractive index is plotted for superluminal waves that correspond to μ = 10, 100, σ = 100, and four different values of the parameter b, that describes the relative amount of phase-averaged magnetic flux carried in the subluminal wave (see Equation (25)). It can be seen that the cutoff moves to greater R as b rises. In the striped wind model, this corresponds to moving away from the equator, where b = 0.

Zoom In Zoom Out Reset image size

Although the refractive indices for non-zero b have similar shape to those for b = 0, there is one potentially important difference: for b > 0.65 the lower branch cuts across the β_> = 0 axis, and continues toward solutions with negative phase speed. Zero refractive index indicates that the H-frame coincides with the laboratory frame. These solutions are standing waves in the pulsar wind, with constant magnetic field and wavelength ≫r_LC. Negative refractive index indicates propagation toward the pulsar, in the sense that the velocity of the H-frame is inward propagating. Nevertheless these modes still carry the particle, energy, momentum, and magnetic fluxes outward toward the nebula.

Conversion of a subluminal wave into a superluminal wave is accompanied by a substantial transfer of energy from the fields to the particles. This is quantified in Figure 5, where the phase-averaged Lorentz factor 〈γ'〉 seen in the laboratory frame is plotted as a function of R, for the same values of μ, σ, and b as in Figure 4. For the smaller b values, b = 0 and b = 0.2, the upper branch of the solution has 〈γ'〉 ≈ μ, which implies that almost all of the field energy is converted to particles in the equatorial regions of the wind. At higher latitudes, only a part of the Poynting flux can be converted, since the phase-averaged magnetic flux cannot be dissipated. However, even at b = 0.95, the upper branch of the dispersion relation carries particles with a mean Lorentz factor one order of magnitude larger than those of the incoming, striped wind solution.

Zoom In Zoom Out Reset image size

According to the above figures, superluminal waves cannot propagate in the inner parts of a pulsar wind. The radius at which propagation becomes possible is close to R = 1 for vanishing phase-averaged magnetic field and rises with this quantity. However, in the striped wind model, the luminosity of the wind, and, hence, the parameter a_L is a function of latitude. In particular, the energy flux carried by the fields is proportional to sin ²θ, where θ is the colatitude (measured from the rotation axis). The flux carried by the particles, on the other hand, is expected to be much smaller, but its angular dependence is uncertain. Assuming, for simplicity, that the ratio σ of these fluxes is independent of colatitude, one has

$$\begin{eqnarray} a_{\rm L}&=&a_{\rm L, eq}\sin \theta, \end{eqnarray} \tag{ 33 }$$

where a_{L, eq} = 4.2 × 10¹⁰L^1/2₃₈ is the value of a_L in the equatorial plane.

The dependence of the average magnetic field carried by the striped wind on latitude depends on the inclination angle α between the magnetic and rotation axes. Two current sheets are present in the striped wind, located at phases

$$\begin{eqnarray} &&\phi =\pm \phi _{\rm sheet}+2n\pi=\pm \hbox{arccos}\left(\hbox{cot}\,\alpha \,\hbox{cot}\,\theta \right) +2n\pi, \end{eqnarray} \tag{ 34 }$$

where n is an integer (Kirk et al. 2002). At θ = atan(cotα), the sheets merge and vanish at higher latitudes (smaller colatitudes θ). At these sheets, the magnetic field reverses direction, so that

$$\begin{eqnarray} b(\theta,\alpha)&=&\left(-\int _{0}^{\phi _{\rm sheet}}\,d\phi +\int _{\phi _{\rm sheet}}^\pi \,d\phi \right)/\pi \nonumber \\ &=&1-2\hbox{ arccos } (\hbox{cot}\,\alpha \,\hbox{cot}\,\theta )/\pi. \end{eqnarray} \tag{ 35 }$$

At those latitudes where there are no sheets, the striped wind model predicts an azimuthally symmetric magnetic field, i.e., no wave.

From curves similar to those plotted in Figure 4 it is possible, for each value of b, to identify the critical radius inside which the mode cannot propagate. This can then be interpolated to give a function R_cr(b). Together with the definition R_cr = ρ_crμ/a_L, this leads to

$$\begin{eqnarray} \frac{\mu \rho _{\rm cr}}{a_{\rm L, eq}}&=&\sin \theta R_{\rm cr} \left[b(\theta,\alpha)\right]. \end{eqnarray} \tag{ 36 }$$

This surface is plotted in Figure 6, for the case of the perpendicular rotator (magnetic inclination angle α = π/2) and, in Figure 7, for various values of α.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The location of the termination shock is determined by the balance between the momentum flux carried by the wind and the external pressure. If we make the reasonable assumption that the latter is independent of latitude, the approximate location ρ_ts(θ) of the shock is given for σ ≫ 1 by (Lyubarsky 2002, Equation (9))

$$\begin{eqnarray} \left(\frac{d \rho _{\rm ts}}{d\theta }\right)^2+\rho _{\rm ts}^2&=&\rho _0^2\sin ^2\theta, \end{eqnarray} \tag{ 37 }$$

where ρ₀ is a constant. The normalization of the radius of the termination shock surface depends on the external pressure and has to be calculated for each individual object. In Figure 6 three examples are shown, labeled by the ratio of the equatorial radius of the shock to that of the critical surface, ρ_{eq, ts}/ρ_{eq, cr}. It can be seen that if this ratio is significantly larger than unity, the cutoff surface falls almost entirely within the termination shock, except for a very small region close to the pole, where both surfaces have a cusp. A qualitatively similar conclusion applies also for different magnetic inclination angles, as shown in Figure 7: for ρ_{eq, ts}/ρ_{eq, cr} ≳ a few, almost the whole critical surface falls within the termination shock radius, apart from a small region at latitudes close to the magnetic inclination angle α. The plots have been calculated using μ = 10, 100 and σ = 100, and show one quadrant of the poloidal plane. The three-dimensional surface is generated by reflection in the ρ axis and rotation about the z-axis.

When the cutoff surface falls inside the termination shock, superluminal waves can propagate in the wind. It is therefore possible, if the ratio ρ_{eq, ts}/ρ_{eq, cr} > 1, that different modes coexist in the same wind but at different latitudes. The MHD description is more appropriate for the higher latitudes (close to α), and the superluminal wave solutions exist closer to the equatorial plane of the wind. If r_{eq, ts}/r_{eq, cr} < 1, no superluminal modes can be supported by the wind at any latitude. In the extreme case where r_{eq, ts}/r_{eq, cr} ≫ 1, practically the whole wind can support these waves.

The above calculations have been conducted for μ ≈ 10⁴. In this case the critical surface in the equatorial plane is located at ρ_cr ∼ 10⁷ for the Crab pulsar, whereas the termination shock is located at roughly ρ_ts ∼ 10⁹. For the Vela pulsar, the critical radius lies at ρ_cr ∼ 10⁶, which is also well within the termination shock, located for this pulsar at ρ_ts ∼ 10⁸ (Kargaltsev & Pavlov 2008). This value of μ corresponds to pair multiplicities of κ = 2.5 × 10⁶ and κ = 2.5 × 10⁵ for the Crab and Vela, respectively, much higher than conventional estimates (de Jager 2007), which would place the critical surface even closer to the pulsar. Thus, at least for young, energetic pulsars, ρ_{eq, ts}/ρ_{eq, cr} ≫ 1 and superluminal waves can propagate in essentially the entire equatorial wind.

However, the situation is different for pulsars which are members of binary systems, since the wind of the companion star can provide an obstacle capable of sustaining the termination shock relatively close to the pulsar. One example is the system PSR B1259-63, containing a 48 ms pulsar in an eccentric orbit about a B2e star. At periastron, the separation of the stars is ρ ∼ 4.3 × 10⁴, while at apastron it becomes ρ ∼ 6 × 10⁵. Very high energy (TeV) gamma-rays were predicted to emerge from the termination shock of this pulsar (Kirk et al. 1999) and subsequently were detected (Aharonian et al. 2005, 2009; Abdo et al. 2011) at binary phases close to periastron. Assuming μ ≈ 10⁴ in this young pulsar, the critical radius at the equator is ρ_cr ∼ 2.4 × 10⁵, which is greater than periastron but smaller than apastron separation. Therefore, if the termination shock created by the interaction of the stellar wind with the pulsar wind lies close to the B2e star, superluminal modes may play a role in determining its properties near apastron, but not near periastron. It would, therefore, be interesting to search for an observable diagnostic of this transition, such as an orbital modulation of the radiation from the shock front.

The simplest solution of these equations is the case of a circularly polarized wave. Then B = 0, and we have

$$\begin{eqnarray} &&u_{\perp } = u_{\perp 0} {\rm e}^{i\omega _{\rm p}t} \end{eqnarray} \tag{ A11 }$$

$$\begin{eqnarray} &&E = E_{0} {\rm e}^{i\omega _{\rm p}t} \end{eqnarray} \tag{ A12 }$$

$$\begin{eqnarray} &&u_{\parallel } = u_0 \end{eqnarray} \tag{ A13 }$$

$$\begin{eqnarray} &&n=n_0. \end{eqnarray} \tag{ A14 }$$

The frequency coincides with the (proper) plasma frequency:

$$\begin{equation} \omega _{\rm p}= \sqrt{\frac{8\pi n_0 e^2}{m}} \end{equation} \tag{ A15 }$$

and the magnitude of the perpendicular component of the four-momentum equals the strength or nonlinearity parameter a, defined using the electric field amplitude as

$$\begin{eqnarray} &&u_{\perp 0}=a\equiv \frac{e\left|E_0\right|}{mc\omega _{\rm p}}. \end{eqnarray} \tag{ A16 }$$

The parallel component of the four-velocity is a constant.

In the laboratory frame, which moves at speed −c/v_ϕ as seen from the H-frame, the frequency of the wave is

$$\begin{equation} \omega ^{\prime } = \Gamma _>\omega _{p0}, \end{equation} \tag{ A17 }$$

where Γ_> is the Lorentz factor of the transformation between the H-frame and the laboratory frame.

Near the equatorial plane of a pulsar wind, the wave is linearly polarized and the relevant field components are the polar electric field and the toroidal magnetic field. In the H-frame then, this requirement corresponds to E_z = 0 and B_y = 0, i.e., E = E_y becomes real and B = iB_z purely imaginary (and constant). It follows that the z-component of the four-velocity is constant Im(u_⊥). In a pulsar, a non-zero u_z corresponds to an azimuthal current loop and, therefore, is proportional to the integrated magnetic flux crossing the enclosed surface. In the split-monopole geometry, this quantity falls off inversely with radius, and we set it to zero, so that u_⊥ is a real variable.

Ampere's law becomes

$$\begin{eqnarray} &&\frac{d E}{d \phi } = -\frac{mc\omega }{e}\frac{u_\perp }{\alpha \gamma }, \end{eqnarray} \tag{ A18 }$$

where ϕ = ωt and

$$\begin{eqnarray} &&\alpha = \frac{\omega ^2}{\gamma _0 \omega ^2_p}. \end{eqnarray} \tag{ A19 }$$

Multiplying the equations of motion (A5) and (A6) by n and integrating, the four-velocity components can be expressed in terms of the variable y = E/E₀, the electric field normalized to the maximum value of its modulus (i.e., −1 ⩽ y ⩽ 1) as follows:

$$\begin{eqnarray} &&u_{\parallel } = u_0 + \alpha a^2 B_z (1-y)/E_0 \end{eqnarray} \tag{ A20 }$$

$$\begin{eqnarray} &&\gamma = \gamma _0 + \frac{\alpha a^2}{2}(1-y^2) \end{eqnarray} \tag{ A21 }$$

$$\begin{eqnarray} &&u^2_{\perp } = \gamma ^2 - u^2_{\parallel } - 1. \end{eqnarray} \tag{ A22 }$$

Having expressed all other unknowns in terms of the electric field, it remains to calculate the phase dependence of y, which is given by the equation

$$\begin{eqnarray} &&\alpha ^2 a^2 \left(\frac{dy}{d\phi } \right)^2 = \frac{u_{\perp }^2}{\gamma ^2} \end{eqnarray} \tag{ A23 }$$

or, explicitly

$$\begin{equation} \alpha ^2 a^2 \left(\frac{dy}{d\phi } \right)^2 = \frac{N(y)}{[q/2 +(1-y^2)]^2} \end{equation} \tag{ A24 }$$

(Kennel & Pellat 1976) with q, Q, and b given by the expressions

$$\begin{eqnarray} &&q = \frac{4 \gamma _0}{\alpha a^2} \end{eqnarray} \tag{ A25 }$$

$$\begin{eqnarray} &&Q = 2 \left(1 - \frac{u_0}{\gamma _0} \lambda \right) \end{eqnarray} \tag{ A26 }$$

$$\begin{eqnarray} &&\lambda = \frac{B_z}{E_0}\hspace{5.0pt} \end{eqnarray} \tag{ A27 }$$

and

$$\begin{eqnarray} &&N(y)=(1-y)^2[(y+1)^2-4\lambda ^2-q] +(1-y)qQ\nonumber\\ \end{eqnarray} \tag{ A28 }$$

is a quartic in y. Periodic waves exist if N has four real roots, which are (in order)

$$\begin{eqnarray} &&y_1=\zeta \cos (\xi -4\pi /3)-1/3 \end{eqnarray} \tag{ A29 }$$

$$\begin{eqnarray} &&y_2=\zeta \cos (\xi -2\pi /3)-1/3 \end{eqnarray} \tag{ A30 }$$

$$\begin{eqnarray} &&y_3=1 \end{eqnarray} \tag{ A31 }$$

$$\begin{eqnarray} &&y_4=\zeta \cos (\xi)-1/3, \end{eqnarray} \tag{ A32 }$$

where

$$\begin{eqnarray} &&\zeta =2(4+12\lambda ^2+3q)^{1/2}/3 \end{eqnarray} \tag{ A33 }$$

$$\begin{eqnarray} \xi &=&\arccos \left[\frac{3\sqrt{3}}{\sqrt{2}(\sqrt{3} u)^{3/2}}\right. \nonumber\\ &&\times\, \left.(16/27+2q/3-2 \lambda q u_0/\gamma _0-16\lambda ^2/3)\vphantom{\frac{3\sqrt{3}}{\sqrt{2}(\sqrt{3} u)^{3/2}}}\right] \end{eqnarray} \tag{ A34 }$$

provided ξ is real.

The nonlinear dispersion relation follows by demanding that the change in phase in a half-cycle of the field oscillation equals π:

$$\begin{equation} \pi = \int _{y_1}^1 \frac{dy}{dy/d\phi }. \end{equation} \tag{ A35 }$$

The phase-averaged value of a quantity A that depends on y is given by

$$\begin{equation} \langle A \rangle = \frac{1}{\pi } \int _{y_1}^1 \frac{A(y)}{dy/d\phi }dy. \end{equation} \tag{ A36 }$$

If the phase-averaged magnetic field vanishes (λ = 0), N(y) reduces to a quadratic in y², with y₂ = −1 and y₃ = 1. In this case, the integrals in (A35) and (A36) can be expressed in closed form in terms of elliptic integrals, at least for the functions of interest. For non-zero λ, on the other hand, y oscillates between y₂ > −1 and y₃ = 1. The dispersion relation and the phase-averages must then be evaluated by numerical integration.

The components of the stress energy tensor consists of a fluid and a field contribution. In the H-frame these are

$$\begin{eqnarray} T^{00}_{\rm part} &=& 2 m c^2 \langle n \gamma ^2 \rangle \, =\, 2 n_0\gamma _0 mc^2 \langle \gamma \rangle \end{eqnarray} \tag{ A37 }$$

$$\begin{eqnarray} T^{01}_{\rm part} &=& 2 m c^2 \langle n \gamma u_x \rangle \,=\, 2 n_0\gamma _0 mc^2 \langle u_x \rangle \end{eqnarray} \tag{ A38 }$$

$$\begin{eqnarray} T^{11}_{\rm part} &=& 2 m c^3 \big\langle n u_x^2 \big\rangle \,=\, 2 n_0 \gamma _0 mc^3 \left< \frac{u_x^2}{\gamma }\right> \end{eqnarray} \tag{ A39 }$$

for the fluids and

$$\begin{eqnarray} T^{00}_{\rm EM} &=& \frac{\langle E^2 \rangle + B ^2}{8 \pi } \,=\, \frac{\langle y^2\rangle + \lambda ^2}{8 \pi } E^2_0 \end{eqnarray} \tag{ A40 }$$

$$\begin{eqnarray} T^{01}_{\rm EM} &=& \frac{c\langle E \rangle B}{4\pi } \,=\, \frac{c\langle y \rangle \lambda E^2_0}{4\pi } \end{eqnarray} \tag{ A41 }$$

$$\begin{eqnarray} T^{11}_{\rm EM} &=& T^{00}_{\rm EM} \end{eqnarray} \tag{ A42 }$$

for the fields. For the circularly polarized modes, λ = 0, <y > =0, and <y² > =1. In addition, u_x = u_x0 and γ = γ₀ are phase independent.