EFFICIENCY OF MAGNETIC TO KINETIC ENERGY CONVERSION IN A MONOPOLE MAGNETOSPHERE

We idealize the central neutron star as a perfectly conducting sphere that we refer to as the "star." We assume that the star has a split-monopole magnetic field configuration, with unit field strength at the stellar surface (this choice sets the energy scale). Exterior to the star, we initialize the system with a low but finite rest-mass density atmosphere, which is done because the code cannot accurately evolve a large density contrast between the initial and injected density or a large σ value in the initial atmosphere. The density of the atmosphere is chosen so that it is dynamically unimportant (the kinetic energy of the piled-up atmosphere is much less than the kinetic energy of the wind). This atmosphere is easily swept away by the outflowing MHD wind and has no effect on the final results. This was confirmed by considering otherwise identical models but where the atmosphere density was 20 times lower. We find all our results are converged indicating negligible impact by the atmosphere on the injected wind.

The initial system has no rotation, so field lines are perfectly radial and both B_θ and B_φ vanish. Starting with this initial configuration, we impose a uniform rotation on the star and study the time evolution of the external magnetosphere. As the star spins up, the footpoints of magnetic field lines are forced to rotate, and this generates a set of outgoing waves traveling at nearly the speed of light. A short distance behind the outgoing wave front, the magnetosphere settles down to a steady state. We are interested in the properties of this steady MHD wind.

The computational domain in our simulations is the upper hemisphere, 0 < θ < π/2, with radius extending from the surface of the star, r = 1, to an outer edge at r = 10¹². We note that most calculations in the literature are limited to a very small radial range due to the need to always resolve the time-dependent compact object (e.g., McKinney & Narayan 2007a). Our choice of a very large range of radius allows us to study both the initial "acceleration zone" of the magnetized wind as well as the asymptotic "coasting zone."

At the polar axis, θ = 0, and at the midplane, θ = π/2, we use the usual antisymmetric boundary conditions. At the outer radial boundary (r = 10¹²) we apply an outflow condition. At the stellar surface, we set the poloidal component of the 3-velocity of the wind to a fixed value v_p directed along the poloidal magnetic field $\vec{B}_p$. We choose v_p = 1/2 in all the simulations reported here. Also, we choose the angular velocity of the star to be Ω = 3/4 (i.e., v_φ = 3/4 at the stellar equator). These boundary conditions are the same for all simulations.

We assume that the magnetic field is frozen into the star. The radial component of the field is continuous across the stellar surface, so we enforce the boundary condition B_r = 1 at r = 1. Since we inject a sub-Alfvénic flow at the surface of the rotating star, Alfvén and fast magnetosonic waves communicate from the magnetosphere back to the surface and generate self-consistent nonzero values of B_θ and B_φ. These values at the stellar surface are set by using "outflow"-type boundary condition (i.e., flowing out of the computational domain into the star) on B_θ and B_φ.

Since we fix v_p and Ω at the stellar surface, the toroidal component of the 3-velocity v_φ is determined by the condition of stationarity:

$$\begin{equation} v_\varphi = \Omega R + B_\varphi v_p/B_p. \end{equation} \tag{ 1 }$$

This equation follows by decomposing the wind velocity into rotation with the field line (the first term) plus motion parallel to the field line (the second term).

The final boundary condition at the stellar surface is the plasma density. This controls how much mass is loaded onto field lines at their footpoints. The different simulations we report in this paper correspond to different choices for ρ(θ_fp), the profile of density as a function of polar angle θ_fp of field line footpoints across the stellar surface.

We solve the time-dependent axisymmetric equations of special relativistic MHD (ignoring gravity) at zero temperature with second-order accuracy. For this we use the code HARM (Gammie et al. 2003; McKinney & Gammie 2004; McKinney 2006b) with recent improvements (Mignone & McKinney 2007; Tchekhovskoy et al. 2007, 2008). We use HARM's second-order MC limiter method for spatial interpolations and a second-order Runge–Kutta time integration. To improve the accuracy of the simulation, before each reconstruction step we interpolate the ratio of the numerical to approximate analytical solution, as described in Appendix C. To speed up the computations, we stop evolving regions of the solution where the wind has achieved a steady state (see, e.g., Komissarov et al. 2007; Tchekhovskoy et al. 2008). This technique allows a large gain in speed of up to a factor ∼10¹⁰, proportional to the maximum radius of the simulation. Also, since our interest is in cold flows in which the plasma internal energy and pressure are negligibly small, we set these quantities (and all their derivatives) identically to zero and ignore the energy evolution equation. For example, the conversion of conserved to primitive quantities is performed identically to the thermal case (Mignone & McKinney 2007), but pressure and internal energy and all their derivatives are set to zero.

HARM is a flexible code that permits the use of an arbitrary coordinate system. We employ a radial grid in which the resolution is very good near the star, but the cells become much more widely spaced at large radii. In terms of a uniformly spaced internal code coordinate x₁, we write the radial coordinate of the grid as

$$\begin{equation} r(x_1) = r_0 + \exp [x_1+c H(x_1-\tilde{x}_1) \times (x_1-\tilde{x}_1)^n], \end{equation} \tag{ 2 }$$

where H(x) is the Heaviside step function. The lower (upper) x₁ value is determined by the lower (upper) radial edge of the grid: r = 1 (r = 10¹²). The parameters r₀, $\tilde{r} = r(\tilde{x}_1)$, c ⩾ 0, and n > 2 allow us to control how resolution varies with r. For instance, at $r_0 \ll r < \tilde{r}$ the grid is logarithmic, with near-uniform relative grid cell size Δr/r ≈ const. At $r = \tilde{r}$ the grid smoothly switches to hyper-logarithmic, with Δr/r ∝ ln^1−1/nr for $r\gg \tilde{r}$. Using a hyper-logarithmic grid at large radii is sufficient since all nontrivial changes in quantities (e.g., the Lorentz factor and collimation angle) occur logarithmically slowly with r.

For the angular coordinate, we map θ to a uniform code coordinate x₂, which goes from 0 to 1, according to

$$\begin{equation} \theta (x_2) = x_2 \frac{\tilde{\theta }}{\tilde{x}_2} + H(x_2-\tilde{x}_2) \left(\frac{\pi }{2} - \frac{\tilde{\theta }}{\tilde{x}_2}\right) \times \left(\frac{x_2-\tilde{x}_2}{1-\tilde{x}_2}\right)^7. \end{equation} \tag{ 3 }$$

In this grid, a fraction $\tilde{x}_2$ (typically ∼1/2) of the cells are distributed uniformly between θ = 0 and $\theta =\tilde{\theta }$ and the remaining cells are distributed nonuniformly at larger angles.

In different models we utilize grids with different values of the parameters r₀, $\tilde{r}$, c, n, $\tilde{x}_2$, and $\tilde{\theta }$. Table 1 gives the details and estimates the effective resolution of our models in terms of the typically used uniform angular grid and exponential radial grid (i.e., r(x₁) = r₀ + exp [x₁]).

Table 1. Simulation Parameters

Name	θmax (°)	μmax	μout	r0	$\tilde{r}$	c	n	$\tilde{\theta }$ (°)	$\tilde{x}_2$	Resolution	Eff. Resolution
Constant density-on-the-star models
M90	90	460	⋅⋅⋅	0	5	0.25	4	11.25	1/4	1536 × 384	3918 × 768
1000M90	90	1000	⋅⋅⋅	0	100	0.25	3	90	1	2048 × 384	2784 × 384
Variable density-on-the-star models
M45	45	460	40	0	100	0.25	3	90	1	2048 × 384	2784 × 384
M20	20	460	40	0.55	100	0.05	3	20	1/3	3072 × 384	3216 × 576
M10	10	460	40	0	5	0.25	4	10	1/3	3072 × 384	7837 × 1152
Variable density-on-the-star + wall at θ = θmax
W10	10	460	40	0.7	100	0.25	3	90	1	2048 × 128	2786 × 128
W5	5	460	40	0.55	100	0.25	3	90	1	6144 × 128	8358 × 128

Download table as:  ASCIITypeset image

For our fiducial baseline model, called M90, we set the density ρ(θ_fp) at the surface of the star to a constant value equal to 0.001914, independent of θ, such that B²_r/(ρ₀c²) ≈ 522. The rotation of the star launches a relativistic magnetized wind that accelerates the input mass. The simulation is run for a time 10¹² in units of the light-crossing time of the star. By this time, the wind reaches a steady-state solution that is independent of initial transients out to a radius 10¹⁰, so all our results are reported out to this radius. The top two panels in Figure 1 show the results.

Zoom In Zoom Out Reset image size

Consider first the shapes of the field lines in the poloidal plane. We see that, over most of the solution, the field lines are only slightly perturbed from their initial purely radial configuration. This might be surprising since rotation causes the toroidal component of the field to grow substantially. In fact, over most of the solution, B_φ is tens of times larger than the poloidal field $B_p\equiv \sqrt{B_r^2+B_\theta ^2}$, and one might think that the hoop stress of the strong toroidal field would cause substantial collimation of the field lines. This does not happen because the relativistic wind sets up an electric field, and an associated electric force per unit volume $\rho _e\vec{E}$ (charge density ρ_e), that almost exactly cancels the magnetic hoop stress. This is a unique feature of relativistic winds. It has been explored in detail by Narayan et al. (2007) and Tchekhovskoy et al. (2008, hereafter TMN08 in connection with winds in the force-free approximation (McKinney 2006a; McKinney & Narayan 2007b).

While it is true that the distortion of field lines in the poloidal plane is small, nevertheless, there is some distortion. It is especially obvious near the rotation axis, where we see field lines converging toward the pole and bunching up. In fact, very close to the rotation axis, the distortion appears to be quite large.⁵ Thus, the results indicate that the polar regions of a rotating monopole wind are qualitatively different from the equatorial regions of the outflow.

We next consider the acceleration of the wind. The top left panel of Figure 1 shows the Lorentz factor γ as a function of position. We see that there is an inner acceleration zone extending out to a radius r ∼ 10², where the Lorentz factor increases from its initial small value (≈1.2 at the stellar surface). Beyond this we find a large coasting zone where there is very little change in γ. Over most of the coasting zone, the flow reaches ultrarelativistic Lorentz factors of γ ∼ 40; however, this value is much less than expected if all the available free magnetic energy were to be used.

An axisymmetric magnetized wind has several conserved quantities along field lines. Two of these are the enclosed magnetic flux Φ and the angular velocity Ω. Another is the ratio of poloidal magnetic flux to rest-mass flux (Chandrasekhar 1956; Mestel 1961; Li et al. 1992; Beskin 1997):

$$\begin{equation} \eta (\Phi) = \frac{\gamma \rho v_p}{B_p} = {\rm const.}{\rm \ along\ field\ line.} \end{equation} \tag{ 4 }$$

Yet another conserved quantity is the quantity μ, which is the ratio of the total energy flux to the rest-mass flux (Chandrasekhar 1956; Mestel 1961; Lovelace et al. 1986; Begelman & Li 1994; Beskin 1997; Chiueh et al. 1998):

$$\begin{equation} \mu (\Phi) = \frac{\mathcal S+\mathcal K}{\mathcal R} = \frac{E \left|B_\varphi \right| + \gamma ^2 \rho v_p}{\gamma \rho v_p} = {\rm const.\ along\ field\ line}. \end{equation} \tag{ 5 }$$

Here, ${\mathcal S}$ is the Poynting flux, ${\mathcal K}$ is the mass energy flux (rest-mass $\mathcal R$ plus kinetic energy), B_φ is the toroidal field,

$$\begin{equation} E=\Omega R B_p \end{equation} \tag{ 6 }$$

is the poloidal electric field, and ρ is the mass density in the comoving frame of the fluid. The denominator of Equation (5) is the rest-mass flux. Since we consider highly magnetized winds, we have ${\mathcal S} \gg {\mathcal K}$ at the surface of the star. Moreover, at r = 1 in the monopolar flow (Michel 1969, 1973),

$$\begin{equation} \left|B_\varphi \right| \approx E \approx \Omega \sin \theta _{\rm fp}, \end{equation} \tag{ 7 }$$

where θ_fp represents the value of θ at the footpoint of a field line. Also, at the footpoint, v_p and ρ are constant, and γ ≈ 1. We thus have⁶

$$\begin{equation} \mu (\theta _{\rm fp}) \propto \sin ^2\theta _{\rm fp}, \end{equation} \tag{ 8 }$$

for μ ≫ γ₀ (γ₀ is the initial Lorentz factor at r = 1), i.e., μ is a rapidly increasing function of θ_fp. This can be seen in the top right panel in Figure 1 and also in Figure 2.

Zoom In Zoom Out Reset image size

As mentioned above, the quantity μ is conserved along each field line. However, the two energy contributions to μ, the Poynting flux ${\mathcal S}$ and the mass energy flux ${\mathcal K}$, are not individually conserved. In fact, the outflowing wind converts Poynting flux to mass energy, thereby accelerating the wind and causing γ to increase.

If the wind were maximally efficient at accelerating the matter, we would expect ${\mathcal S}\rightarrow 0$ at large distance from the star. The Lorentz factor would then achieve its maximum value

$$\begin{equation} \gamma _{\rm max} = \mu. \end{equation} \tag{ 9 }$$

In practice, the wind falls far short of this maximum. As Figures 1 and 2 show, μ is quite large in model M90, with a value of 460 at the equator (θ_fp = π/2). However, the Lorentz factor of the wind, even at a radius of 10¹⁰, does not exceed 40. Thus, a magnetized monopole wind is very inefficient at accelerating the gas. For an equatorial field line in M90, the efficiency factor γ/μ is only about 0.09, as shown in Figure 2.

Another way of describing the efficiency of conversion of energy from electromagnetic to kinetic form is via the magnetization parameter σ, which is the ratio of Poynting to mass energy flux (Kennel & Coroniti 1984a; Li et al. 1992; Begelman & Li 1994; Vlahakis 2004; Komissarov et al. 2007),

$$\begin{equation} \sigma = \frac{\mathcal S}{\mathcal K} = \frac{E \left|B_\varphi \right|}{\gamma ^2 \rho v_p}. \end{equation} \tag{ 10 }$$

Substituting into Equation (5) we see that the conserved quantity μ is related to σ by

$$\begin{equation} \mu = \gamma (\sigma +1). \end{equation} \tag{ 11 }$$

The smallest value possible for σ is zero. Therefore, the maximum value of γ is γ_max = μ (Equation 9).

The magnetization σ is not conserved along a field line. At the surface of the star, where γ ≈ 1, we have σ ≈ μ − 1. As the magnetized wind flows out and energy is transferred from Poynting to matter energy, γ increases and σ decreases. An efficient wind would be one in which σ asymptotes to a value ⩽1, so that the outflowing material is able to convert at least half of its energy flux into matter energy. The numerical solutions shown in Figures 1 and 2 fail to satisfy this criterion by a large factor in the equatorial regions. This implies there is no ideal, axisymmetric MHD solution to the σ problem (Rees & Gunn 1974; Kennel & Coroniti 1984a, 1984b).

There is, however, one promising feature in the results: the polar regions of the wind are efficient, with γ/μ → 1 and σ ≪ 1 at large distance from the star (Figure 2). This interesting feature of the monopole problem has not been emphasized in the literature. Most previous analyses and discussions have focused on equatorial field lines where the efficiency is, indeed, too low to solve the σ problem.

Unfortunately, the actual Lorentz factor along polar field lines in M90 is only ∼20 since this is the value of μ for these lines. Would we continue to have high efficiency in the polar region even with larger values of μ? In particular, is it possible to have acceleration with high efficiency up to Lorentz factors γ>100, as observed for instance in GRBs? For this we need to study a model with larger values of μ near the pole. We describe such a model in the following subsection.

From Equation (5), we see that an obvious way to increase μ is to lower the density of the wind at the stellar surface. For instance, if we were to reduce ρ by a factor of ∼30 relative to M90, then we would have a model with μ∼  few hundred for a field line with θ_fp ∼ 10° (see Figure 3). We could then explore acceleration along this field line and determine whether or not the outflowing wind achieves a coasting γ>100.

Zoom In Zoom Out Reset image size

Since μ varies as sin²θ_fp, this approach would lead to extremely large values of μ at the equator. As a result, the model would require very large resolution to simulate accurately and would be extremely expensive. Therefore, for numerical convenience, we consider a model in which we choose the profile ρ(θ_fp) such that μ is large near the pole, reaches a maximum μ_max at a specified footpoint angle θ_fp = θ_max, and then decreases with increasing θ_fp to an outer value μ_out at θ_fp = π/2.⁷ To achieve this, we choose the density profile on the star to be

$$\begin{equation} \rho (\theta _{\rm fp}) = \rho _0 + \rho _1 \sin ^{\alpha } \theta _{\rm fp}, \end{equation} \tag{ 12 }$$

where the constants ρ₀, ρ₁, and α are adjusted so that the model has the desired values of μ_max, θ_max, and μ_out.

We have simulated a series of such models, all with μ_max = 460 and μ_out = 40,⁸ and with different values of θ_max, viz., 45°, 20°, 10°, which we refer to as models M45, M20, M10, respectively. Figure 3 shows the profiles of ρ(θ_fp) and μ(θ_fp) for the model M10. Parameters of the various models are summarized in Table 1. The models are well converged, with field-line invariants Ω(Φ), η(Φ), and μ(Φ) conserved along field lines to better than 15%, as Figure 4 shows for, e.g., models M90 and M10.

Zoom In Zoom Out Reset image size

The bottom two panels in Figure 1 show results corresponding to M10. As before, we see that the poloidal structure of the field is largely unaffected by rotation. There is of course some lateral shift of field lines, the effect being larger near the pole than near the equator, with maximum field line bunching close to the axis.

The Lorentz factor distribution again confirms the trends seen in M90. Near the equator, the asymptotic γ is only ∼12, giving an inefficient flow with γ/μ ∼ 0.3. The efficiency increases near the axis and becomes practically equal to unity close to the pole; equivalently, σ becomes much less than unity for these field lines. Most interestingly, Lorentz factors nearly as large as ∼200 are obtained for these field lines. In other words, there is no σ problem near the axis and it is possible to obtain quite large Lorentz factors in this region of the outflow.

From the results shown in Figures 1 and 2, we conclude the following.

• 1.  
  Field lines near the equator of a rotating monopole largely retain their monopolar configuration out to large radii, whereas lines near the pole tend to bunch up around the axis.
• 2.  
  The acceleration efficiency γ/μ of a rotating monopole magnetosphere is low (≪1) for field lines in the equatorial region, but quite high (∼1) for field lines near the pole.

We would like to develop an understanding of the physics behind of these effects. We would also like to know how the two effects are related to each other. This is the topic of the following two sections.

As plasma streams along field lines, relativistic effects become important near the so-called light cylinder, R = R_L = 1/Ω, where the corotation velocity ΩR equals the speed of light and E = B_p (see Equation (6)). As we show in Appendix A, far outside the light cylinder, where ΩR ≫ 1 and γ ≫ γ₀, the plasma simply drifts perpendicular to $\vec{E}$ and $\vec{B}$ at the drift velocity (Beskin et al. 1998, 2004; Vlahakis 2004),

$$\begin{equation} v \approx v_{\rm dr} = \left|\frac{\vec{E} \times \vec{B}}{B^2}\right| = \frac{E}{B}. \end{equation} \tag{ 13 }$$

The corresponding Lorentz factor is

$$\begin{equation} \gamma ^2 \approx \gamma _{\rm dr}^2 = \frac{B^2}{B^2-E^2}. \end{equation} \tag{ 14 }$$

Near the star this formula becomes inaccurate since the plasma moves at the initial Lorentz factor,

$$\begin{equation} \gamma ^2 \approx \gamma _0^2. \end{equation} \tag{ 15 }$$

In Appendix A, we show that a combination of these two formulae does a very good job of describing the Lorentz factor at all distances from the star:

$$\begin{equation} \gamma ^2 = \gamma _0^2 - 1 + \gamma _{\rm dr}^2. \end{equation} \tag{ 16 }$$

In the asymptotic region of a relativistic outflow (where ΩR ≫ 1, γ ≫ 1), we have according to Equation (14),

$$\begin{equation} \left|B_\varphi \right| \approx E. \end{equation} \tag{ 17 }$$

Further, this relation is also true at the surface of the central compact star (for the monopolar flow, see Equation (7)). This allows us to write the difference between the maximum allowed Lorentz factor γ_max = μ and the local Lorentz factor γ in the following convenient form (the numbers in parentheses refer to the equations used to derive this result):

$$\begin{equation} \mu -\gamma \mathop{\approx}^{({5}),({7}),({17})} \frac{E^2}{\gamma \rho v_p} \mathop{=}^{({4})} \frac{(\Omega R B_p)^2}{\eta B_p} = \frac{\Omega ^2(\Phi)}{\eta (\Phi)} {B_p R^2}. \end{equation} \tag{ 18 }$$

By dividing this equation by itself as evaluated at the footpoint, and approximating γ − γ_fp ≈ γ and μ − γ_fp ≈ μ, we obtain

$$\begin{equation} \frac{\gamma }{\mu } \approx 1 - \frac{B_p R^2}{[B_p R^2]_{\rm fp}} \approx 1 - \frac{a}{a_{\rm fp}}, \end{equation} \tag{ 19 }$$

where we have used the subscript "fp" to denote quantities evaluated at the field line footpoint⁹ and have defined the quantity

$$\begin{equation} a \equiv B_p R^2. \end{equation} \tag{ 20 }$$

Equation (19) demonstrates that, in order to convert an appreciable fraction of the total energy flux μ along a field line into matter energy flux (γ times the mass flux), the quantity a has to decrease appreciably from its initial value at the footpoint. This result is known (compared to Begelman & Li 1994; Chiueh et al. 1998; Vlahakis 2004), though we have not seen as simple a derivation as the one given above.

For a precisely monopole field, B_p ∝ 1/R² along each field line. Therefore, a = B_pR² is constant along a field line and so no efficient acceleration is possible. This explains why acceleration is so difficult in the monopole problem. In order to permit acceleration, field lines must move in a cooperative fashion transverse to one another so as to allow a to reduce with increasing distance from the star. We discuss how this is accomplished in the following subsection.

Meanwhile, as an aside, we describe here an improvement to the approximate result (Equation (19)) which gives the correct numerical value of efficiency in a precisely monopolar field. For such a field, the Lorentz factor at asymptotically large distances is (Michel 1969; Camenzind 1986)

$$\begin{equation} \gamma _\infty ^\mathrm{radial} = \mu ^{1/3} \ll \mu. \end{equation} \tag{ 21 }$$

This gives an extremely inefficient acceleration, γ/μ ≈ μ^−2/3 ≪ 1. However, the efficiency is not zero as Equation (19) might suggest, so we need a more accurate version of Equation (17). According to Equation (14),

$$\begin{equation} B_\varphi \approx \frac{\gamma E}{\sqrt{\gamma ^2-1}} = \frac{\gamma \Omega R B_p}{\sqrt{\gamma ^2-1}} \approx \Omega R B_p \left(1+\frac{1}{2\gamma ^2}\right). \end{equation} \tag{ 22 }$$

Using this equation, assuming B_pR² = const. along a field line, and introducing γ_∞ = lim_R→∞γ, we obtain in the limit R → ∞

$$\begin{equation} \mu (\Phi) \approx \gamma _\infty +\frac{\Omega ^2(\Phi)}{\eta (\Phi)}B_p R^2\left(1+\frac{1}{2\gamma _\infty ^2}\right). \end{equation} \tag{ 23 }$$

Now, assuming that the system settles down to a state with minimum total energy flux μ(Φ) (equivalent to the minimal torque condition of Michel 1969), we find the terminal Lorentz factor γ_∞ that minimizes the right-hand side of this equation:

$$\begin{equation} \gamma _\infty ^\mathrm{radial} =\left[\frac{\Omega ^2(\Phi)}{\eta (\Phi)}B_p R^2\right]^{1/3} \approx \mu ^{1/3}, \end{equation} \tag{ 24 }$$

where the approximate equality comes from Equation (18) evaluated at the footpoint. This reproduces Equation (21) and shows that indeed, for a precisely radial flow, only a small (but still nonzero) fraction of the Poynting flux is converted to the kinetic energy of the matter. This derivation was performed for precisely monopolar field lines. In actuality, the shape of the poloidal field lines is slightly changed from radial. The effect is small for equatorial field lines and the order of magnitude estimate (24) continues to hold. The deviations are larger for polar field lines and the Lorentz factor obtained along these lines is very different from Equation (24).

Let us first apply Formula (19) to a field line near the midplane. For acceleration to be efficient, a must decrease along the field line, i.e., B_p must decrease faster than 1/R². This can be accomplished by moving field lines away from the equator toward the axis. Consider a field line with its footpoint located at a small angle θ'_fp = π/2 − θ_fp ≪ 1 from the equator. For this field line, we can write Equation (19) as

$$\begin{equation} \frac{\gamma }{\mu }\approx 1-\frac{a}{a_{\rm fp}} \approx 1-\frac{\theta ^{\prime }_{\rm fp}}{\theta ^{\prime }}, \end{equation} \tag{ 25 }$$

where θ' is the polar angle of the field line at a large distance from the star, where the Lorentz factor is γ. Here we used the fact that Φ' ≈ 2πa θ' ≈ 2πa_fpθ'_fp, where Φ' is the amount of flux enclosed between the field line and the midplane.

For a nearly monopolar configuration of the field in which θ' ≈ θ'_fp, obviously we will not have much acceleration (γ/μ ≪ 1). In order to obtain high acceleration efficiency in the equatorial region, field lines must diverge from the equator so that the θ' values of field lines increase with distance. For this to happen, the rest of the magnetosphere must collectively move away from the equator toward the pole. For reasons that are discussed in Section 5, this does not happen, and so acceleration efficiency is at best modest near the equator.

The story is quite different for field lines close to the axis (θ ≪ 1). Consider the initial undistorted monopole configuration. At the surface of the star, B_p is constant, and [B_pR²]_fp is simply equal to Φ/π, where Φ is the flux interior to the field line. For a pure monopole, this relation is valid at any distance from the star, i.e., B_pR² = Φ/π at all radii, and therefore a/a_fp = 1 and acceleration would be inefficient.

In analogy with the previous discussion for equatorial field lines, let us now imagine uniformly expanding or contracting the field lines near the axis. That is, for each field line with a given θ_fp ≪ 1, let the polar angle far from the star become θ = kθ_fp, with the same value of k for all lines. For such a uniform expansion or contraction of the field, B_p transforms to B_p/k². However, at the same time R becomes kR, and so a = B_pR² is unaffected. In other words, there is no effect on acceleration.

The key to obtaining acceleration along polar field lines is not uniform lateral expansion (divergence) or contraction (collimation) of field lines, but differential bunching of field lines. To see this, rewrite Equation (19) as

$$\begin{equation} \frac{\gamma }{\mu } \approx 1-\frac{a}{a_{\rm fp}} \approx 1 - \frac{\pi B_p R^2}{\Phi }, \end{equation} \tag{ 26 }$$

where we have used the fact that a_fp ≈ Φ/π near the pole. Clearly, for efficient acceleration, we must make B_p substantially smaller than the mean enclosed field Φ/πR². That is, the field lines interior to the reference field line must be bunched in such a way that most of the flux has been pulled inward. As an example, consider a power-law distribution of the field strength,

$$\begin{equation} B_p(R) \propto R^{-\xi }, \end{equation} \tag{ 27 }$$

where the index ξ measures the degree of bunching. This distribution gives

$$\begin{equation} \frac{\gamma }{\mu } \approx \frac{\xi }{2}, \end{equation} \tag{ 28 }$$

which shows that the acceleration efficiency increases with increasing ξ. We reach equipartition between Poynting and matter energy flux (σ ∼ 1, γ ∼ μ/2) for ξ = 1, and we obtain arbitrarily large efficiency (σ ≪ 1) as ξ → 2.

As we have described in Section 3, polar field lines in simulation M90 are very efficient with γ/μ → 1. According to the above, this would seem to suggest that ξ(R) should be ≈2 at the axis and should decrease with increasing R. The simulations, however, show that this does not happen: as we show later, ξ ≈ 1, i.e., it stays roughly constant over a range of R. Instead, higher efficiency near the jet axis is achieved in a different way: as field lines bunch around the jet axis (Chiueh et al. 1991; Eichler 1993; Bogovalov 1995; Bogovalov & Tsinganos 1999), they form a concentrated core (Heyvaerts & Norman 1989; Bogovalov 2001; Lyubarsky & Eichler 2001; Beskin & Nokhrina 2008) that takes up a finite amount of flux Φ₀, leading roughly to the following poloidal field strength profile

$$\begin{equation} B_p(R) = \Phi _0 \delta (\pi R^2) + B_0 (R/R_0)^{-\xi }, \end{equation} \tag{ 29 }$$

where we have approximated the core profile with the Dirac delta function. This gives

$$\begin{equation} \frac{\gamma }{\mu } \approx \left\lbrace \begin{array}{l@{\quad }l} {\xi }/{2}, & \Phi _0 \ll \pi B_pR^2 / (1-\xi /2), \cr 1 - {\pi B_pR^2}/{\Phi _0}, & \Phi _0 \gg \pi B_pR^2 / (1-\xi /2). \end{array}\right. \end{equation} \tag{ 30 }$$

That is, in the limit when the flux in the concentrated core is small compared to the flux in the surrounding power-law field distribution, the efficiency is the same as in Equation (28). However, in the opposite limit, i.e., sufficiently close to the core where the flux in the core dominates, the angular profile of the poloidal field distribution (and the value of ξ) become irrelevant for determining the acceleration efficiency. Thus, (1) differential bunching and (2) the resulting development of a concentrated core, are the key requirements for efficient acceleration.

Note the following important corollary from the above discussion. It does not matter whether the particular field line of interest collimates toward the axis or diverges from it. This has no effect on the acceleration. What we need is that (1) other field lines closer to the axis must converge more, or diverge less, compared to the reference field line, and/or (2) a concentrated core at the jet axis must contain a significant amount of magnetic flux.

Figure 5 shows results for models M90 and M10. The top panels show the behavior of a/a_fp at different distances from the central star. We see that a/a_fp decreases toward the pole, exactly where γ/μ is largest and σ is smallest in Figure 2. Also, Figure 2 quantitatively confirms the validity of Equation (19) by plotting prediction (19) over the numerical solution.

Zoom In Zoom Out Reset image size

The lower panels in Figure 5 illustrate the effects described in the previous two subsections. In the equatorial regions, we see that B_pr² decreases from its footpoint value of unity. It is this decrease that allows whatever acceleration is observed in this region of the outflow. However, the decrease is modest, so the acceleration is not very large.

For angles closer to the pole, B_pr² actually increases relative to its nominal initial value of unity. Nevertheless, this does not mean that there is deceleration because, as we argued above, acceleration near the axis is associated with differential bunching, not with any overall expansion or contraction. For both M90 and M10, we see that the magnetic field sets up the required bunching so that the poloidal field is maximum at the axis and decreases with increasing distance from the pole. It is this outward decrease, coupled with the presence of magnetic flux Φ₀ in a concentrated core at the pole (see Equation (30)), that is associated with a decrease in a/a_fp for polar field lines and the reason for strong acceleration.

The effect is most clearly seen in model M90, as illustrated in Figure 5. Between r = 10² and r = 10³ the angular profile of B_p becomes a steeper function of polar angle, leading to an associated decrease in a/a_fp. Between r = 10³ and r = 10⁴ the trend reverses, and the poloidal field profile actually becomes a shallower function of polar angle; however, a/a_fp continues to decrease. Beyond r ∼ 10⁴ the magnetic field profile does not evolve further. In this regime, it can be well fitted by a broken power law,

$$\begin{equation} B_p r^2 = 0.85 + 0.075 \sin ^{-1.1}\theta. \end{equation} \tag{ 31 }$$

Therefore, at these distances, a = B_pr²sin²θ does not evolve with r either. Despite this, a/a_fp ≈ πa/Φ decreases with increasing r and the value of Φ increases. This is solely due to the increase in the amount of magnetic flux Φ₀ contained in the concentrated core.

The picture that emerges is the following. The poloidal magnetic field establishes some equilibrium angular profile B_p(r, θ)r² at intermediate latitudes that does not evolve with r. At progressively larger r, each individual magnetic field line becomes progressively more collimated toward the pole, but the same angular profile of magnetic field is maintained. This means that magnetic flux Φ flows out from the low-latitude equatorial region (B_p slightly decreases there), flows through intermediate latitudes without changing the profile of B_p there, and ends up in the concentrated core at the polar axis, thereby uniformly shifting the angular poloidal flux distribution Φ(θ) up. The effect is clearly seen in the lower left panel of Figure 5.

Note that it is not necessary to numerically resolve the concentrated core in order to accurately describe the jet structure since it is only the total amount of magnetic flux contained in the concentrated core that matters for the acceleration efficiency. In particular, our numerical method is well suited for capturing such a core, even if unresolved, since our method conserves the magnetic flux to machine precision and can accurately capture the amount of magnetic flux that enters the core and remains there. We note that within a few grid cells from the polar axis, where the unresolved magnetic flux accumulates, are least accurate but this does not affect the quality of the solution at larger angles. To verify this, we have checked convergence of our models with angular resolution. For this, we ran a version of model M90 that uses a uniform angular grid and has a factor of 2 lower effective angular resolution near the pole. Using this less-resolved model leads to a maximum relative difference in the flux function Φ and Lorentz factor of less than 15%, even near the rotation axis. This difference is less than 2% at most radii (r < 10² and r > 10⁵) and is smaller at larger θ. This confirms the accuracy of the numerical solution. Future higher resolution models that resolve smaller angles and the concentrated core are required to determine how the solution connects to the polar axis and for independent verification of results.

Beskin et al. (1998) have discussed the physical reason for inefficient acceleration in the equatorial regions. They show that it is related to the fast magnetosonic point. In the comoving frame of a cold MHD plasma, fast magnetosonic waves travel with a speed v_f given by (Gammie et al. 2003; McKinney 2006b)

$$\begin{equation} \gamma _f v_f = \left({b^2\over \rho }\right)^{1/2}, \end{equation} \tag{ 32 }$$

where b is the comoving magnetic field strength. It is straightforward to relate b to field components in the lab frame:

$$\begin{equation} b^2 = B^2-E^2 \mathop{\approx}^{({14})}\frac{B^2}{\gamma ^2} \mathop{\approx}^{({17})}\frac{E \left|B_\varphi \right|}{\gamma ^2} \mathop{\approx}^{({10})}\rho \sigma. \end{equation} \tag{ 33 }$$

We thus find

$$\begin{equation} \gamma _f v_f \approx \sigma ^{1/2}. \end{equation} \tag{ 34 }$$

Consider a streamline in the wind that moves outward with a local Lorentz factor γ. Let us first consider the limit of infinitely high magnetization, ρ → 0, i.e., the force-free limit (Goldreich & Julian 1969; Okamoto 1974; Blandford 1976; Lovelace 1976; Blandford & Znajek 1977; MacDonald & Thorne 1982; Fendt et al. 1995; Komissarov 2001, 2002a, b; McKinney 2006a; Narayan et al. 2007, TMN08). In this limit γ_f → ∞, therefore throughout the solution we have γ < γ_f. This means that the fast magnetosonic surface, defined by the condition γ = γ_f, where the wind becomes causally detached from fluid farther back along its streamline, is located at infinity. Such a force-free wind has a simple analytic solution in which the Lorentz factor increases roughly linearly with distance (Michel 1973),

$$\begin{equation} \gamma \approx \Omega R. \end{equation} \tag{ 35 }$$

In general, ρ ≠ 0, yet we might expect that the behavior of the Lorentz factor in the subfast region γ < γ_f is similar to the force-free solution (35). This has been shown to indeed be the case (Beskin et al. 1998). However, the acceleration in the super-fast region γ>γ_f has been found to become logarithmic, i.e., inefficient (Beskin et al. 1998). Once the wind has crossed the fast magnetosonic point γ = γ_f, it becomes causally detached from the fluid farther back along its streamline. We might therefore expect efficient acceleration to cease beyond this fast magnetosonic point¹⁰ for all field lines.

If we define the fast Mach number M_f by

$$\begin{equation} M_f = \frac{\gamma v}{\gamma _fv_f} \approx \frac{\gamma v}{\sigma ^{1/2}}, \end{equation} \tag{ 36 }$$

then the fast point is the location at which M_f ≈ 1 (this equality would be exact if there was only motion along the poloidal field line; however, there is also a slow rotation in the toroidal direction which introduces a negligible correction that we ignore). For a relativistic flow (v ≈ 1), Equation (11) lets us recast the above expression in a useful form:

$$\begin{equation} M_f^2 \approx \frac{\gamma ^2}{\sigma }= \frac{\gamma ^{3}}{\mu -\gamma }. \end{equation} \tag{ 37 }$$

Figure 6 shows results for a field line in M90 with θ_fp = θ = π/2. Until the flow reaches the fast magnetosonic point, we see that σ falls rapidly and γ increases rapidly. However, both trends slow down substantially once the flow crosses the fast point. Beyond this point, Beskin et al. (1998) have shown that σ and γ vary as the one-third power of log r. We confirm this dependence below.

Zoom In Zoom Out Reset image size

The relatively abrupt cessation of acceleration beyond the fast point for equatorial field lines is obvious in Figure 1, where we see that γ stops increasing once the flow crosses the fast magnetosonic surface (the middle of the three thick solid lines). However, it is also clear from Figure 1 that something else operates on polar field lines. Model M10, in particular, shows substantial continued acceleration well after polar field lines have crossed the fast surface. We discuss next the relevant physics for these field lines.

We showed in Section 4 that, for efficient acceleration along a field line, other neighboring field lines must shift laterally. At the equator, we need lines to move away from the midplane, while near the pole, we need field lines to experience differential bunching or develop a flux core. In order for any given field line to sustain efficient energy conversion, it must be able to communicate to other regions of the magnetized wind that undergo differential bunching or cause a concentrated flux core. This suggests that the fast magnetosonic point, which determines where the fluid can no longer communicate back along its motion, is perhaps not so important. A more relevant issue is whether or not the fluid can communicate with regions near the axis that have field bunching or a flux core. We refer to the point at which a fluid element loses contact with the axis as the "causal point," and call the locus of causal points over all field lines as the "causality surface." By the above arguments, we expect that this surface, rather than the fast magnetosonic surface, plays the role of the boundary for efficient acceleration.¹¹ We expect efficient acceleration inside this surface and inefficient, logarithmic acceleration outside the surface (for a related discussion, see Zakamska et al. 2008; Komissarov et al. 2009). In general fast waves propagate away from any given point toward the rotation axis through an intermediate region where the density, velocity, and magnetic field vary. Hence, one should trace the position of fast waves emitted from any given point outward over all angles and identify the "causal point" for each field line as where finally no such traces can reach the polar axis. For simplicity, we instead use only the local fast wave speed at a given point, and we identify the approximate "causal point" by where the locally emitted fast waves move away from any given point with a lab-frame local angle of θ>0 with respect to the rotation axis, so that the waves do not reach the rotation axis over a finite propagation distance. We now calculate the approximate location of the causality surface using this approach.

Consider a segment of the relativistic magnetized wind propagating with a velocity vector $\vec{v}$ and Lorentz factor γ at an angle θ_j to the rotation axis.¹² In Appendix D, we show that fast magnetosonic waves, which are emitted isotropically in the comoving frame of the fluid, in the lab frame will be collimated along $\vec{v}$ into a Mach cone with a half-opening angle

$$\begin{equation} \sin \xi _{\rm max} = \frac{\gamma _f v_f}{\gamma v} = \frac{1}{M_f}. \end{equation} \tag{ 38 }$$

By the argument given earlier, field line bunching and efficient acceleration are possible only when the fluid can communicate with the axis, i.e., only if ξ_max ≳ θ_j, i.e., only if

$$\begin{equation} \sin \theta _j \lesssim \frac{1}{M_f}. \end{equation} \tag{ 39 }$$

For an equatorial wind, i.e., θ_j = π/2, Equation (39) shows that acceleration stops when M_f = 1, i.e., γ = γ_f ≈ σ^1/2 (compare to Equation (34), assuming v_f → 1). That is, the wind stops efficient acceleration as soon as it crosses the fast magnetosonic point. However, for smaller values of θ_j, we obtain a different result.

According to Equation (38), communication with the axis and efficient acceleration are possible until

$$\begin{equation} \gamma v\approx {\gamma _fv_f\over \sin \theta _j} = {\sigma ^{1/2} \over \sin \theta _j}. \end{equation} \tag{ 40 }$$

The presence of the factor sin θ_j in the denominator means that communication extends to larger values of γ, i.e., acceleration efficiency becomes larger as θ_j decreases. In other words, polar field lines can accelerate more easily. Using the definition of the fast magnetosonic Mach number (Equation (36)), we obtain the following relation for the causality surface:

$$\begin{equation} \sin ^2\theta _{j,c} \mathop{=}^{({40})} \frac{1}{M_{f,c}^2} \mathop{\approx}^{({37})} \frac{\sigma _c}{\gamma _c^2} = \frac{\mu -\gamma _c}{\gamma _c^3}, \end{equation} \tag{ 41 }$$

where the subscript "c" indicates quantities evaluated at the causality surface. Inside the causality surface we expect the Lorentz factor to increase roughly linearly with distance, compare to Equation (35). Based on our earlier arguments, once outside the causality surface the acceleration will only be logarithmic. Using Equation (41), we can estimate the distance at which the causality surface is located:

$$\begin{equation} r_c \sim \frac{\mu ^{1/3}}{\Omega \sin ^{5/3}\theta _{j,c}}, \end{equation} \tag{ 42 }$$

where we have assumed that Equation (35) and γ ≪ μ hold for r ≲ r_c.

Figure 1 confirms that the causality surface (the outermost of the three thick solid lines) provides a better approximation to the boundary between the acceleration and coasting zones compared to the fast magnetosonic surface. Figure 7 shows detailed results for a polar field line with sin θ_fp = 0.1 in model M10. Note that efficient acceleration continues well past the fast magnetosonic point ("F," dotted line on the left); acceleration slows down only after the field line has crossed the causal point ("C," dotted line on the right). This is the reason why this particular field line is able to achieve a Lorentz factor of 170 with a high efficiency of γ/μ ≈ 0.6, which is much larger than for equatorial field lines.

Zoom In Zoom Out Reset image size

We now develop an analytical approximation to calculate the Lorentz factor as a function of distance for any field line in a monopole magnetized wind. Generally, for a force-free jet, there exist two distinct acceleration regimes, as explained in TMN08. In Appendix A, we generalize these results to MHD (finite-magnetization) jets. We summarize the results here. In the first acceleration regime, which is realized near the compact object, the Lorentz factor of the flow increases roughly linearly with distance:

$$\begin{equation} \gamma _1^2 \approx \gamma _0^2 + (\Omega R)^2, \end{equation} \tag{ 43 }$$

where γ₀ is the initial Lorentz factor at r = 1.

Based on earlier arguments, beyond the causality surface the acceleration is only logarithmic. This is the second acceleration regime in which the Lorentz factor is determined by the poloidal shape of the field lines (Beskin et al. 1998; TMN08). Using the results of Appendix B as a guide (see also Beskin et al. 1998; Lyubarsky & Eichler 2001), we expect in this region

$$\begin{equation} \gamma _2 \propto \ln ^{1/3}r. \end{equation} \tag{ 44 }$$

To make this formula quantitative, we demand that it gives the correct value of the Lorentz factor at the causality point r = r_c (see Equation (41)):

$$\begin{equation} \gamma _c \approx \mu ^{1/3} \sin ^{-2/3}\theta _{j,c}. \nonumber \end{equation} \tag{ 45 }$$

We do this by choosing the solution in the following form:

$$\begin{eqnarray} \gamma _2 &\approx & C_1 \, \gamma _c \, \ln ^{1/3}(1+C_2r/r_c) \end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray} &\approx & C_1 \, \mu ^{1/3} \frac{\ln ^{1/3}(1+C_2r/r_c)}{\sin ^{2/3}\theta _{j,c}}, \end{eqnarray} \tag{ 47 }$$

where C₁ and C₂ are numerical factors of order unity that we later determine by fitting to the numerical solution. The radial and angular scalings in Equation (47) agree with the analytic expectations (Beskin et al. 1998; Lyubarsky & Eichler 2001; Lyubarsky 2009).

Note, however, that Formulae (44)–(47) become inconsistent at low magnetization since γ cannot exceed μ, whereas the right-hand sides of these equations are unbound. Noting that the fast wave Mach number M_f is unbound and γ/μ^1/3 ≈ M^2/3_f for γ ≪ μ (Equation (37)), we empirically modify Equation (47) by replacing γ₂/μ^1/3 with M^2/3_f:

$$\begin{equation} M_f^{2/3} = C_1 \frac{\ln ^{1/3}(1+C_2r/r_c)}{\sin ^{2/3}\theta _{j,c}}, \end{equation} \tag{ 48 }$$

where C₁ and C₂ are numerical factors of order unity (see below). Substituting for the fast wave Mach number M_f using Equation (37), we obtain a cubic equation for the Lorentz factor on the field line beyond the causality surface, r ≳ r_c:

$$\begin{equation} \gamma _2 = C_1 \, (\mu -\gamma _2)^{1/3}\frac{\ln ^{1/3}(1+C_2 r/r_c)}{\sin ^{2/3}\theta _{j,c}}, \end{equation} \tag{ 49 }$$

where μ is the value of the total specific energy flux on the field line in question (Equation (5)), θ_j,c is the value of the angle θ_j that the field lines makes with the polar axis at the causality surface (see footnote 8), and C₁ ≃ C₂ ≃ 1 (see below). In the limit γ₂ ≪ μ this equation reduces to Equation (47).

We now combine the two approximations, Equations (43) and (49), to write (see Appendix A)

$$\begin{equation} \frac{1}{\gamma ^2} = \frac{1}{\gamma _1^2} +\frac{1}{\gamma _2^2}. \end{equation} \tag{ 50 }$$

Clearly, the smaller of γ₁ and γ₂ determines the total Lorentz factor: in accordance with the above discussion, near the compact object γ ≈ γ₁ and at a large distance (outside the causality surface) γ ≈ γ₂. We find that we obtain good agreement with our simulation results when we choose C₁ = 2, C₂ = 0.4. In fact, for this single set of parameters Formula (50), with γ₁ and γ₂ given by (43) and (49), does quite well for all field lines in all simulations, both in the limit of low and high magnetizations. The various dotted lines in Figures 6–10 have all been calculated using this formula, and clearly provide an excellent representation of the numerical results.

We have so far discussed in detail the representative models M90 and M10. However, we have carried out a number of other simulations. We mentioned models M20 and M45 in Section 3.2. We have also carried out models with walls: W5, W10.

Simulation W10 has the same setup as M10 but has an impenetrable perfectly conducting wall at θ = 10°. Figure 8 shows that the wall in simulation W10 keeps field lines near the wall from collapsing onto the pole and prevents the rest of the field lines from developing the lateral nonuniformity required for efficient acceleration. As a result, field lines very close to the wall accelerate more efficiently but the rest of the field lines have a suppressed efficiency: a field line with sin θ_fp = 0.1 in model W10 has a lower Lorentz factor at r = 10¹⁰, γ ≈ 150, than the corresponding field line in model M10, which has γ ≈ 170 (compare Figures 7 and 8). Simulation W5, which has the wall at θ = 5°, is the most collimated model that we have simulated. Figure 9 shows a field line for that model that makes an angle sin θ_fp = 0.06 at the surface of the central star. This field line reaches equipartition by r ∼ 10¹⁰ with γ ≈ 200, the largest Lorentz factor we have achieved among all simulations in this paper.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have also performed a simulation called 1000M90 that has a uniform density profile on the stellar surface ρ(θ_fp) = ρ₀ with the same density at θ_fp = 0 as model M45. Thus, models 1000M90 and M45 have similar profiles of density and μ near the pole, but they differ near the equator. We expect the two simulations to show nearly identical behavior for polar field lines. This is indeed confirmed, as seen in Figure 10. The point of this model is to verify that models with variable density, e.g., Equation (12), give reliable results near the axis, independent of how we modify the mass-loading of equatorial field lines. The numerical results confirm that this is indeed so.

Zoom In Zoom Out Reset image size

Figure 10 shows transversal cuts through each of our models at distances of 10⁶ and 10⁹, and compares the simulation results with the analytic approximation Equation (50) described in Section 5.3. For models with a constant density profile ρ(θ_fp) on the star, the analytic approximation works extremely well at all distances r and all polar angles. For models with a variable angular profile of density on the surface of the star, the agreement is excellent along the field lines originating in the constant-density core while for other field lines, the best-fit values of factor C₁ in Equation (49) apparently varies from one field line to the next. Therefore, adopting a single value C₁ = 2 provides only a rough description of acceleration along these field lines.

Equations (40) and (52) predict that at roughly the same value of σ ∼ 1, an asymptotically more collimated simulation reaches a larger Lorentz factor:

$$\begin{equation} \gamma _\mathrm{asym}\propto \frac{1}{\sin \theta _{j,c}}, \end{equation} \tag{ 53 }$$

where θ_j,c ∼ θ_fp ≳ θ_j is the angle at the causality surface. Figure 10 confirms this for the sequence of models M45–M20–M10, all of which have the same maximum value of μ, μ_max = 460. If all of these models could convert all of the electromagnetic energy flux into kinetic energy flux, each of them would reach the maximum energetically allowed Lorentz factor, γ = μ_max. While such a full conversion does not happen in any of these models, more collimated models reach higher Lorentz factors, in agreement with Equation (53): each subsequent model in the sequence of models M45–M20–M10 is roughly twice as collimated as the previous one and at r = 10⁹ attains approximately twice as large a Lorentz factor, 50–95–170.

The net conclusion of the previous discussion is that, whereas there is indeed a serious σ problem for pulsar winds, there is no similar problem for relativistic jets. We show that even if a flow is unconfined and has a monopolar-like shape, it can still efficiently accelerate in the polar region. For a jet angle θ_j ∼ 2°, for instance, Equation (52) predicts an asymptotic Lorentz factor γ_asym > 100 if the jet efficiently converts electromagnetic to kinetic energy and reaches σ ∼ 1. In fact, if the jet is not efficient and carries more of its energy as Poynting flux (e.g., for a larger value of μ, see Equations (52) and (51)), then σ>1, and we will have even larger values of γ_asym. These estimates are confirmed by the numerical simulations described in this paper.

The inferred total power of long GRBs is on the order of 10⁵¹ erg (Piran 2005; Meszaros 2006; Liang et al. 2008); however, much less energetic events, with total energy release as low as 10⁴⁸ erg, have also been observed (Soderberg et al. 2004, 2006). Let us compute the power output of jets in our numerical models and make sure that our jets are energetic enough to be consistent with these observations. The total power coming out from the compact object surface within the low-density core, θ < θ_max, is

$$\begin{eqnarray} P^{\rm jet} &=& \int \limits _0^{\theta _{\rm max}} 2\pi r^2d\theta \sin \theta \, S(r,\theta)|_{r=1} \approx \frac{\Omega ^2}{2}\int \limits _0^{\theta _{\rm max}} d\theta \sin ^3\theta \nonumber \\ &=& \frac{\Omega ^2}{2} \times \frac{4(2+\cos \theta _{\rm max})}{3} \times \sin ^4\frac{\theta _{\rm max}}{2} \approx \frac{\Omega ^2}{8} \,\theta _{\rm max}^4, \end{eqnarray} \tag{ 54 }$$

where the last equality is for θ_max ≪ 1. Converting the result to physical units, we obtain

$$\begin{equation} P^{\rm jet} \approx \frac{1}{8}\Omega ^2 B_r^2 r_0^4 \theta _{\rm max}^4/ c, \end{equation} \tag{ 55 }$$

where B_r is the value of radial magnetic field component on the surface of the compact object and r₀ is the compact object radius. Evaluating the jet power for a magnetar with a characteristic period P = 1 ms and a surface magnetic field of 10¹⁵ G, we get

$$\begin{eqnarray} P^{\rm jet}_{\rm NS} &\approx & 2 \times 10^{47}\;({\rm erg\;s}^{-1})\nonumber\\ &\times& \left(\frac{1\,\rm{ms}}{P}\right)^2 \left(\frac{B_r}{10^{15}\,\rm{G}}\right)^2 \left(\frac{r_0}{10\,\rm{km}}\right)^2 \left(\frac{\theta _{\rm max}}{10^\circ }\right)^4. \end{eqnarray} \tag{ 56 }$$

For a maximally spinning black hole with dimensionless spin parameter a = 1, mass 3 M_☉, and surface magnetic field strength 10¹⁶ G (McKinney 2005), we have

$$\begin{equation} P^{\rm jet}_{\rm BH} \approx 5 \times 10^{48} \, a^2 \left(\frac{B_r}{10^{16}\,\mathrm{G}}\right)^2 \left(\frac{M}{3\;M_\odot }\right)^2 \left(\frac{\theta _{\rm max}}{10^\circ }\right)^4 \left[\frac{\mathrm{erg}}{\mathrm{s}}\right]. \end{equation} \tag{ 57 }$$

Given a characteristic duration of 10–100 s for a long GRB, our numerical jets provide, for the black hole case, 10⁵⁰–10⁵¹ erg per event for model M10, 10⁵¹–10⁵² erg for M20, and 10⁵²–10⁵³ erg for M45. Therefore, the simulated jets from black holes in models M10 and M20 are energetic enough and move at sufficiently high Lorentz factors (γ ≳ 100) to account for most long GRBs and can certainly account for low luminosity events. The energetics is also right for short GRBs: the simulated jets output 10⁴⁹–10⁵⁰ erg during a characteristic event duration of 1 s (Nakar 2007). For the magnetar case the energetics is lower: 10⁴⁷–10⁴⁸–10⁵⁰ erg s⁻¹ for a sequence of models M10–M20–M45. Therefore, the simulated jets from magnetars can account for less-luminous long GRB events and most short GRBs.

While the energy fraction in the polar jet is small as compared to the total energy extracted by magnetic fields from the black hole, P^jet/P^tot ∼ θ⁴_max (Equation (55), which assumes a uniform magnetic field distribution at the BH horizon), the absolute value of jet power P^jet is sufficiently large to account for long and short GRBs. We point out that the magnetic flux in accreting black hole systems is nonuniformly concentrated in the polar region of the BH (e.g., McKinney 2005, due to ambient pressure of the accretion flow), and therefore the total energy losses of the spinning black hole are actually dominated by the losses from the polar region rather than from the midplane region, meaning a larger fraction of power in the jet than given by Equation (55).

A very interesting question is whether the models suggest any characteristic value for the quantity γ_jθ_j: is this quantity generally smaller or larger than unity For a jet with γ_jθ_j ≫ 1 (say, ∼10) only part of the jet within the beaming angle θ_b ≈ 1/γ_j ≪ θ_j is visible to a remote observer. As the interaction with the ambient medium decelerates the jet and the beaming angle becomes comparable to the jet opening angle, the edges of the jet come into sight and the light curve steepens achromatically, displaying a "jet break" (Piran 2005; Meszaros 2006). In the other limit, γ_jθ_j ≪ 1, a jet would be incapable of producing achromatic breaks in GRB light curves.

Achromatic breaks have been found in pre-Swift times (Frail et al. 2001; Bloom et al. 2003; Zeh et al. 2006) and have yielded θ_j ∼ 0.03–1 rad. The situation is quite different in the post-Swift times for which a large amount of data available, and no jet breaks have been found to fully satisfy closure relations in all bands (Liang et al. 2008). However, if one or more of the closure relations are relaxed, some of the breaks may be interpreted as "achromatic" and be used to derive jet opening angles that span a similar range as pre-Swift GRBs (Liang et al. 2008). Overall, it appears that the Lorentz factor of most GRBs is γ_j ≳ 100 (Piran 2005; Meszaros 2006, up to ∼400, Lithwick & Sari 2001) which, with the above estimates for θ_j, gives γ_jθ_j ≳ 3 and indicates that in principle achromatic jet breaks are possible.

Our numerical models have γ_jθ_j ∼ 5–15, where the index j denotes quantities evaluated at the jet boundary which we define as the boundary of matter-dominated region σ < 1. In this respect, it is particularly fruitful to compare simulations M10 and W10. As we discussed in Section 5.4, the wall in model W10 prevents field lines from collapsing onto the pole as much as they do in model M10. Comparison of field lines with σ ∼ 1 for these models (see Figures 7 and 8), reveals a difference in the Lorentz factor of at most 20% and a much larger difference in the value of γ_jθ_j (caused by a large difference in θ_j due to the effect of the wall): γ_jθ_j ≈ 5 for M10 and ≈12 for W10. The most collimated model W5 produces an even larger value, γ_jθ_j ≈ 15.

We now analytically confirm that in general the quantity γ_jθ_j ∼ 5–10. For the jet boundary, using Equation (52) and characteristic values Ω ∼ 1, θ_j,c ∼ 0.04–0.4, μ ∼ 10³, r ∼ 10⁶–10⁹, θ_j/θ_j,c ∼ 0.5–1 (see Figures 7–9), we get

$$\begin{equation} \gamma _j\theta _j \sim \frac{4\theta _j}{\theta _{j,c}} \sigma ^{1/2}\log _{10}^{1/2} \left(\frac{r}{10^2\hbox{--}10^3}\right) \sim (4\hbox{--}10)\sigma ^{1/2} \end{equation} \tag{ 58 }$$

which is in good agreement with the simulation results. The large value of γ_jθ_j in this analysis arises solely due to the logarithmic factor that appears because a significant fraction of the acceleration occurs after crossing the causality surface, in the inefficient acceleration region: this is the case for all unconfined flows studied in our paper. This should be contrasted with confined outflows, collimated by walls with prescribed shapes, for which most of the acceleration tends to complete before crossing the causality surface and which have γ_jθ_j ∼ 1 (Komissarov et al. 2009). We note that if we do not require that jets are matter-dominated, e.g., we allow σ>1 as in Lyutikov & Blandford's (2003) model of GRBs, then the value of γ_jθ_j will be even higher (see Equation (58)).

According to Equation (58), we can attribute the fact that some post-Swift GRBs show quasi-achromatic jet breaks, while many do not, by associating the former with jets that have γ_jθ_j > 1 and the latter with those that have γ_jθ_j < 1. Such a scatter in γ_jθ_j might be naturally produced by differences in GRB environment (affecting θ_j/θ_j,c, see Equation (58)) or the properties of the central engine (affecting σ). Indeed, according to Equation (58), a low value of σ (≪1) or θ_j (≪θ_j,c) in our jets would mean γ_jθ_j < 1 and so the absence of a jet break.

One could use Equation (51), which is based on our analytical model of the simulations, to obtain the Lorentz factor of any GRB jet. If the black hole or neutron star is nearly maximally spinning with Ω ∼ 0.25 and has a polar region with the reasonable value of μ ∼ 1000 (TMN08), and if we consider an opening angle θ_j ∼ 4°, then by r ∼ 10⁸ ∼ 10¹⁴ cm one obtains γ_j ∼ 250 and γ_jθ_j ∼ 17. For a range of opening angles with sufficient luminosity, one obtains a range of Lorentz factors consistent with both short- and long-duration GRB jets (Piran 2005; Meszaros 2006). Further, the product γ_jθ_j ≳ 1, indicating an afterglow can exhibit the so-called "achromatic jet breaks," where observations imply γ_jθ_j ≳ 3 for long-duration GRB jets (Piran 2005; Meszaros 2006). Our simulations and analytical models have γ_jθ_j ∼ 5–15, which is a proof of principle that magnetically driven jets can produce jet breaks.