Model for how an accretion disk drives astrophysical jets and sheds angular momentum

Observational evidence strongly suggests that accretion disks and astrophysical jets are closely coupled and that disks act as the driver for jets. There is a large, rich literature on the topic of astrophysical jets and disks. We will not attempt to review this literature here other than mentioning a few relevant papers; a more extensive discussion of relevant literature is given in the introductions of [1] and [3]. There exists a variety of astrophysical jets which are distinguished by size and velocity. The smaller jets, called protostellar or Herbig–Haro jets, have a size of the order of the solar system, are non-relativistic and are associated with the early life of a star and the formation of planets. The larger jets are highly relativistic, have dimensions many orders of magnitude greater than protostellar jets, and are associated with black holes and neutron stars. A common feature of all jets is that the jet is associated with an accretion disk and this feature suggests that existence of an accretion disk is a necessary condition for existence of a jet.

Most of the astrophysical literature addresses either disks or jets with relatively few papers discussing the coupling of disks to jets. In one of the earliest jet papers Blandford and Zjanek [4] proposed a model for the highly relativistic situation of a jet coming from a black hole. They assumed that magnetic flux was frozen into accreting material so that accretion would concentrate flux and that the disk was force-free. They noted the need for a radial electric potential drop in the disk plane and suggested there must be an effective finite resistivity to enable accreting material to cross magnetic field lines (i.e., not be frozen to the magnetic field). Their mechanism involved extracting energy from the rotation of the black hole to drive the jet.

The disk papers have primarily focussed on the magneto rotational instability (MRI). The original version of the MRI was presented by Balbus and Hawley [5] who assumed ideal MHD (i.e., Ohm's law given by ${\bf{E}}+{\bf{U}}\times {\bf{B}}=0$ so there is zero resistivity) and a 100% ionized plasma undergoing Kepler rotation as the initial condition. Accreting particles contain angular momentum which must be shed because otherwise centrifugal force will reflect the particles at small radius. Classic viscosity in principle could transport angular momentum from small to large radius but the quantitative value of classical viscosity is negligible and so there has been a search for some non-viscous mechanism for transporting angular momentum outwards. The MRI papers argue that MRI turbulence provides the desired mechanism for outward transport of angular momentum in an accretion disk but generally do not address how disks drive jets. Blaes and Balbus [6] noted that protostellar accretion disks are weakly ionized so that the interaction between ions and neutrals would modify the MRI. Gammie [7] argued that if the ionization is very small, the effective resistivity can become so large as to suppress the MRI and denoted regions where the ionization is negligible as 'dead zones'. Sano and Stone [8, 9] and Lesur [10] noted that the Hall term (which is missing from the ideal MHD equations used in [5–7]) becomes important when the disk is weakly ionized and affects the behavior of the MRI instability. The MRI investigations of the disk typically do not use the actual geometry but instead model a small portion of the disk where this portion, known as a shearing box, has an imposed velocity shear that mimics the radial dependence of the Kepler orbital velocity.

Jet models typically use ideal MHD and use the disk region as a boundary condition (this was done for example by Zhai et al [3]). However, a few papers such as those by Zanni et al [11], by Sheikhnezami et al [12], and by Stepanovs and Fendt [13] present models that couple disks and jets with the aim of showing how disks launch jets. These disk-jet coupling models invoked the ideal MHD equations for the jet and resistive MHD equations for the disk. The resistivity was made anomalously large and this enhancement in magnetic diffusivity was presumed to result from MRI turbulence. The value of the invoked anomalous resistivity was not calculated from first principles, but rather was chosen to provide plausible numerical results in numerical calculations. Furthermore, [11–13] did not take into account weak ionization or Hall physics.

Pandey and Wardle [14] considered the effect of the Hall term in a weakly ionized plasma from a more general point of view and noted that the effective ion cyclotron frequency is reduced by the fractional ionization which means that the effective ion cyclotron frequency can become very small. Since the cross-over from the MHD regime to the whistler regime occurs in the vicinity of the ion cyclotron frequency, they noted that weak ionization has the consequence of greatly lowering the frequency at which whistler waves occur.

Because the MHD equations have no intrinsic scale, these equations characterize situations ranging over many orders of magnitude and, in particular, the same dimensionless equations can characterize both a lab experiment and an astrophysical situation. A system for scaling from lab experiments to astrophysical situations was given by Remington et al [15] while Staff et al [16] noted 'the many similarities between laboratory jets and astrophysical jets despite the immense difference in scale speaks to the fundamental nature of the MHD jets'. This similarity means that important insight into how jets work can be gained from studies of laboratory jets since, in principle, all aspects of laboratory jets can be controlled and all parameters can be measured. Jet lab experiments have been reported by Hsu and Bellan [17], You et al [18], Ciardi et al [19], Kumar et al [20], Suzuki-Vidal et al [21], Moser and Bellan [22], and Li et al [23].

The purpose of this paper is to present an integrated disk/jet model in the context of protostellar jets. This model does not depend on the MRI instability and so is fundamentally different from the MRI-based models because there is no invocation of turbulence or anomalous magnetic diffusivity. Furthermore, unlike the shearing box models typically used in MRI models of disks, the geometry here is global and contains both the disk and the jet regions extending to infinity. The model depends on previously unrecognized particle kinetics and involves the new concept that a select group of particles have trajectories substantially different from Kepler or cyclotron motion. These trajectories are at the cross-over from Kepler to cyclotron orbits. A thought experiment indicating that such orbits must exist consists of imagining that a spacecraft slowly becomes electrically charged until its charge-to-mass ratio becomes the same as that of an ion. The spacecraft orbit would then have to change from being a Kepler orbit to being a cyclotron orbit and at the change-over there would have to be a trajectory that is neither Kepler nor cyclotron. Analysis shows that these special trajectories result from fundamental Hamiltonian orbit considerations missed when using the much coarser ideal MHD description. The existence of the group of particles having non-Kepler, non-cyclotron trajectories can equivalently be derived from the Hall MHD equations for a weakly ionized plasma.

The jet dynamics are governed by the ideal MHD equations

$$\begin{eqnarray}&&\rho \left(\displaystyle \frac{\partial {\bf{U}}\ }{\partial t}+{\bf{U}}\cdot {\rm{\nabla }}{\bf{U}}\right)={\bf{J}}\times {\bf{B}}-{\rm{\nabla }}P,\end{eqnarray} \tag{ 3a }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\bf{U}})=0,\end{eqnarray} \tag{ 3b }$$

$$\begin{eqnarray}&&{\bf{E}}+{\bf{U}}\times {\bf{B}}=0,\end{eqnarray} \tag{ 3c }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\times {\bf{E}}=-\displaystyle \frac{\partial {\bf{B}}}{\partial t},\end{eqnarray} \tag{ 3d }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\times {\bf{B}}={\mu }_{0}{\bf{J}}.\end{eqnarray} \tag{ 3e }$$

Since the components of the convective derivative are

$$\begin{eqnarray}&&{({\bf{U}}\cdot {\rm{\nabla }}{\bf{U}})}_{r}={U}_{r}\displaystyle \frac{\partial {U}_{r}}{\partial r}+{U}_{z}\displaystyle \frac{\partial {U}_{r}}{\partial z}-\displaystyle \frac{{U}_{\phi }^{2}}{r},\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&{({\bf{U}}\cdot {\rm{\nabla }}{\bf{U}})}_{\phi }={U}_{r}\displaystyle \frac{\partial {U}_{\phi }}{\partial r}+{U}_{z}\displaystyle \frac{\partial {U}_{\phi }}{\partial z}+\displaystyle \frac{{U}_{\phi }{U}_{r}}{r},\end{eqnarray} \tag{ 4b }$$

$$\begin{eqnarray}&&{({\bf{U}}\cdot {\rm{\nabla }}{\bf{U}})}_{z}={U}_{r}\displaystyle \frac{\partial {U}_{z}}{\partial r}+{U}_{z}\displaystyle \frac{\partial {U}_{z}}{\partial z},\end{eqnarray} \tag{ 4c }$$

on defining the total force ${\bf{F}}={\bf{J}}\times {\bf{B}}-{\rm{\nabla }}P$, the MHD equation of motion can be expressed as

$$\begin{eqnarray}&&\rho \left(\displaystyle \frac{\partial {U}_{r}}{\partial t}+{U}_{r}\displaystyle \frac{\partial {U}_{r}}{\partial r}+{U}_{z}\displaystyle \frac{\partial {U}_{r}}{\partial z}-\displaystyle \frac{{U}_{\phi }^{2}}{r}\right)={F}_{r},\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}&&\rho \left(\displaystyle \frac{\partial {U}_{\phi }}{\partial t}+{U}_{r}\displaystyle \frac{\partial {U}_{\phi }}{\partial r}+{U}_{z}\displaystyle \frac{\partial {U}_{\phi }}{\partial z}+\displaystyle \frac{{U}_{\phi }{U}_{r}}{r}\right)={F}_{\phi },\end{eqnarray} \tag{ 5b }$$

$$\begin{eqnarray}&&\rho \left(\displaystyle \frac{\partial {U}_{z}}{\partial t}+{U}_{r}\displaystyle \frac{\partial {U}_{z}}{\partial r}+{U}_{z}\displaystyle \frac{\partial {U}_{z}}{\partial z}\right)={F}_{z}.\end{eqnarray} \tag{ 5c }$$

The MHD continuity equation is

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}(\rho {{rU}}_{r})+\displaystyle \frac{\partial }{\partial z}(\rho {U}_{z})=0.\end{eqnarray} \tag{ 6 }$$

Adding ${U}_{r},{U}_{\phi },$ or U_z times equation (6) to the respective $r,\phi ,$ or $\ z$ component of equation (5) gives

$$\begin{eqnarray}&&\displaystyle \frac{\partial \ }{\partial t}(\rho {U}_{r})+\ {\rm{\nabla }}\cdot (\rho {U}_{r}{\bf{U}})={F}_{r}+\displaystyle \frac{\rho {U}_{\phi }^{2}}{r},\end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \ }{\partial t}(\rho {{rU}}_{\phi })+{\rm{\nabla }}\cdot (\rho {{rU}}_{\phi }{\bf{U}})={{rF}}_{\phi },\end{eqnarray} \tag{ 7b }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \ }{\partial t}(\rho {U}_{z})+\ {\rm{\nabla }}\cdot (\rho {U}_{z}{\bf{U}})={F}_{z}.\end{eqnarray} \tag{ 7c }$$

Using ${\bf{J}}\times {\bf{B}}={{\bf{J}}}_{{\rm{tor}}}\times {{\bf{B}}}_{{\rm{pol}}}+{{\bf{J}}}_{{\rm{pol}}}\times {{\bf{B}}}_{{\rm{tor}}}+{{\bf{J}}}_{{\rm{pol}}}\times {{\bf{B}}}_{{\rm{pol}}}$ the force vector can be expressed as

$$\begin{eqnarray}{\bf{F}} & = & \displaystyle \frac{1}{4{\pi }^{2}}\left[\displaystyle \frac{1}{{\mu }_{0}}{\rm{\nabla }}\cdot \left(\displaystyle \frac{1}{{r}^{2}}{\rm{\nabla }}\psi \right){\rm{\nabla }}\psi -\displaystyle \frac{{\mu }_{0}I}{{r}^{2}}{\rm{\nabla }}I\right.\\ & & -{\rm{\nabla }}\phi {\rm{\nabla }}\psi \cdot ({\rm{\nabla }}I\times {\rm{\nabla }}\phi )]-{\rm{\nabla }}P\end{eqnarray} \tag{ 8 }$$

so equations (7a)–(7c) become

$$\begin{eqnarray}&&\displaystyle \frac{\partial \ }{\partial t}(\rho {U}_{r})+\ {\rm{\nabla }}\cdot (\rho {U}_{r}{\bf{U}})\\ &&\quad =\,\displaystyle \frac{1}{4{\pi }^{2}}\left(\displaystyle \frac{1}{{\mu }_{0}}{\rm{\nabla }}\cdot \left(\displaystyle \frac{1}{{r}^{2}}{\rm{\nabla }}\psi \right)\ \displaystyle \frac{\partial \psi }{\partial r}-\displaystyle \frac{{\mu }_{0}I}{{r}^{2}}\displaystyle \frac{\partial I}{\partial r}\right)\\ &&\quad -\,\displaystyle \frac{\partial P}{\partial r}+\displaystyle \frac{\rho {U}_{\phi }^{2}}{r},\end{eqnarray} \tag{ 9a }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \ }{\partial t}(\rho {{rU}}_{\phi })+{\rm{\nabla }}\cdot (\rho {{rU}}_{\phi }{\bf{U}})=\displaystyle \frac{1}{4{\pi }^{2}}({\rm{\nabla }}I\times {\rm{\nabla }}\psi \cdot {\rm{\nabla }}\phi ),\end{eqnarray} \tag{ 9b }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \ }{\partial t}(\rho {U}_{z})+\ {\rm{\nabla }}\cdot (\rho {U}_{z}{\bf{U}})\\ &&\quad =\,\displaystyle \frac{1}{4{\pi }^{2}}\left(\displaystyle \frac{1}{{\mu }_{0}}{\rm{\nabla }}\cdot \left(\displaystyle \frac{1}{{r}^{2}}{\rm{\nabla }}\psi \right)\ \displaystyle \frac{\partial \psi }{\partial z}-\displaystyle \frac{{\mu }_{0}I}{{r}^{2}}\displaystyle \frac{\partial I}{\partial z}\right)\\ &&\quad -\,\displaystyle \frac{\partial P}{\partial z}.\end{eqnarray} \tag{ 9c }$$

We now make the following simplifying, but realistic and relevant assumptions:

• 1.  
  As in the derivation of the Grad–Shafranov equation [24, 25] we assume that $I=I(\psi )\$ in which case ${\rm{\nabla }}I\times {\rm{\nabla }}\psi =0$ so the current flows on and not across poloidal flux surfaces. This means that the right hand side of equation (9b) vanishes so there is no force in the ϕ direction. As discussed in [1] we assume ${U}_{\phi }=0$ in the jet initially and since there is no force in the ϕ direction in the jet, ${U}_{\phi }$ will remain zero in the jet. The simplest nontrivial dependence is to have

  $$\begin{eqnarray}&&{\mu }_{0}I=\lambda \psi \end{eqnarray} \tag{ 10 }$$

  which is closely related to force-free equilibria but not exactly the same because there is no assumption that the toroidal current is proportional to the toroidal flux.
• 2.  
  We consider the geometry in figure 2 and in particular, assume that near the z axis, toroidal field pressure has stretched the poloidal magnetic field axially away from its original vacuum profile. The flux function ψ is assumed to be highly distended relative to its vacuum state so poloidal field lines near the z axis are nearly parallel to the z axis. The axial magnetic field ${B}_{z}={(2\pi r)}^{-1}\partial \psi /\partial r$ is thus nearly uniform in the vicinity of the z axis; this uniformity corresponds to the requirement that near the z axis the poloidal flux must scale as r² to satisfy mathematical regularity. Because the axial field is nearly uniform near the z axis, if $a(z)$ is the assumed current channel radius, the poloidal flux for $r\leqslant a$ can be approximated as

  $$\begin{eqnarray}&&\psi (r,z)={\psi }_{0}{\left(\displaystyle \frac{r}{a(z)}\right)}^{2},\end{eqnarray} \tag{ 11 }$$

  where the dependence of a on z is weak so $\partial /\partial z\ll {a}^{-1}$. The poloidal field is assumed to have an umbrella- or mushroom-like shape with equation (11) describing the stalk only. At z = 0 this stalk intercepts what will be called the center electrode and so the center electrode has radius $a(0)$. At large distances the poloidal field bends around to form the cap of the mushroom and then returns back to the z = 0 plane at radii much larger than $a(0)$. This large radius part of the z = 0 plane will be called the outer electrode.
• 3.  
  We assume ${U}_{r}\ll {U}_{z}$ and also that U_r is so small that it can be neglected.

Equations (9a)–(9c) are an exact re-statement of equation (3a) for an axisymmetric system. We seek a steady-state solution, i.e., a solution where $\partial /\partial t=0$. and first note that for the poloidal flux specified by equation (11),

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \left(\displaystyle \frac{1}{{r}^{2}}{\rm{\nabla }}\psi \right)=\displaystyle \frac{{\psi }_{0}}{{r}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {z}^{2}}{\left(\displaystyle \frac{1}{a(z)}\right)}^{2}.\end{eqnarray} \tag{ 12 }$$

This shows that ${\rm{\nabla }}\cdot ({r}^{-2}{\rm{\nabla }}\psi )$ can be ignored compared to other terms in equations (9a) and (9c) because $\partial /\partial z$ is small compared to other terms In equation (9c) the other terms are first order in $\partial /\partial z$ but the term involving ${\rm{\nabla }}\cdot ({r}^{-2}{\rm{\nabla }}\psi )$ is second order in $\partial /\partial z$ and so still can be dropped. The system is solved by considering the radial, azimuthal, and axial components of the equation of motion separately and then combining the results.

Radial component: Taking into account the assumptions ${U}_{r}\simeq 0,$ ${U}_{\phi }=0,$ and using equation (10), equation (9a) becomes

$$\begin{eqnarray}&&0=-\displaystyle \frac{1}{4{\pi }^{2}}\displaystyle \frac{{\lambda }^{2}\psi }{{\mu }_{0}{r}^{2}}\displaystyle \frac{\partial \psi }{\partial r}-\displaystyle \frac{\partial P}{\partial r}\end{eqnarray} \tag{ 13 }$$

which is essentially a balancing between the outward radial force from the pressure and the inward radial magnetic pinch force from $-{J}_{z}{B}_{\phi }$ since ${J}_{z}\sim {r}^{-1}\partial \psi /\partial r$ and ${B}_{\phi }\sim \psi /r$. Using equation (11), equation (13) can be expressed as

$$\begin{eqnarray}&&\displaystyle \frac{{\lambda }^{2}{\psi }_{0}^{2}}{2{\pi }^{2}{\mu }_{0}{(a(z))}^{4}}\ r+\displaystyle \frac{\partial P}{\partial r}=0\end{eqnarray} \tag{ 14 }$$

which can be integrated using the boundary condition P = 0 at r = a to give

$$\begin{eqnarray}&&P(r,z)=\ \displaystyle \frac{{\lambda }^{2}{\psi }_{0}^{2}}{4{\pi }^{2}{\mu }_{0}{(a(z))}^{2}}\left(1-\displaystyle \frac{{r}^{2}}{{(a(z))}^{2}}\right).\end{eqnarray} \tag{ 15 }$$

It should be noted here that the pressure given by equation (15) at r = 0 scales as $1/{(a(z))}^{2}$ and so is higher when a(z) is smaller. This implies existence of an axial pressure gradient that will tend to drive flows from regions of small a(z) towards regions of large $a(z)$.

Azimuthal component: Equation (9b) is satisfied because ${U}_{\phi }$ is assumed zero and $I=I(\psi )$ so ${\rm{\nabla }}I\times {\rm{\nabla }}\psi =0$.

Axial component: Using ${U}_{r}\simeq 0\$ and ${U}_{\phi }=0$ equation (9c) becomes

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial z}(\rho {U}_{z}^{2})=-\displaystyle \frac{{\mu }_{0}I}{4{\pi }^{2}{r}^{2}}\displaystyle \frac{\partial I}{\partial z}-\displaystyle \frac{\partial P}{\partial z}\end{eqnarray} \tag{ 16 }$$

which can be rearranged as a Bernoulli-like relation

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial z}\left(\rho {U}_{z}^{2}+\displaystyle \frac{{\mu }_{0}{\lambda }^{2}{\psi }^{2}}{8{\pi }^{2}{\mu }_{0}{r}^{2}}+P\right)=0.\end{eqnarray} \tag{ 17 }$$

On the z-axis ${\psi }^{2}/{r}^{2}=0$ and so on on the z-axis equation (17) can be integrated to give the conventional Bernoulli relation

$$\begin{eqnarray}&&{(\rho {U}_{z}^{2}+P)}_{r=0}={\rm{const}}.\end{eqnarray} \tag{ 18 }$$

Evaluation of $\rho (r,z),$ ${U}_{z}(r,z)$ and $P(r,z)$ on the z-axis axis at the two positions z = 0 and z = L gives

$$\begin{eqnarray}&&\rho (0,0){U}_{z}^{2}(0,0)+P(0,0)=\rho (0,L){U}_{z}^{2}(0,L)+P(0,L).\end{eqnarray} \tag{ 19 }$$

Next, assume that the axial position L is so large that $a(L)\gg a(0)$ and also assume that ${U}_{z}(0,0)$ is nearly zero. This means that $\rho (0,0){U}_{z}^{2}(0,0)\ll P(0,0)$ and, since P scales as ${a}^{-2},$ that $P(0,L)\ll P(0,0)$. Thus, equation (19) becomes

$$\begin{eqnarray}&&P(0,0)=\rho (0,L){U}_{z}^{2}(0,L).\end{eqnarray} \tag{ 20 }$$

From equation (15) it is seen that

$$\begin{eqnarray}&&P(0,0)=\displaystyle \frac{{\lambda }^{2}{\psi }_{0}^{2}}{4{\pi }^{2}{\mu }_{0}{(a(0))}^{2}}=\displaystyle \frac{{\mu }_{0}{I}^{2}}{{(2\pi a(0))}^{2}}\end{eqnarray} \tag{ 21 }$$

so equation (20) can be solved to give

$$\begin{eqnarray}&&{U}_{z}(0,L)=\displaystyle \frac{1}{\sqrt{{\mu }_{0}\rho (0,L)}}\displaystyle \frac{{\mu }_{0}I}{2\pi a(0)}.\end{eqnarray} \tag{ 22 }$$

Equation (22) can also be expressed as

$$\begin{eqnarray}&&{U}_{z}(0,L)=\displaystyle \frac{{B}_{\phi }(0,a(0))}{\sqrt{{\mu }_{0}\rho (0,L)}}\end{eqnarray} \tag{ 23 }$$

which is formally like the velocity of an Alfvén wave, but involves only ${B}_{\phi }$ rather than the total magnetic field $B=\ \sqrt{{B}_{\phi }^{2}+{B}_{z}^{2}}$ and, furthermore, involves the magnetic field and the mass density evaluated at different locations.

The middle term in equation (17) is finite off the z axis and corresponds to the ${J}_{r}{B}_{\phi }\sim -\partial {B}_{\phi }^{2}/\partial z\$ axial force since ${J}_{r}\sim {r}^{-1}\partial \psi /\partial z$ and ${B}_{\phi }\sim {r}^{-1}\psi$. Equation (17) also shows that the generation of axial kinetic energy ultimately comes from $\partial a/\partial z$ being finite. This is because if a were independent of z then both ${\psi }^{2}$ and P would be independent of z and so there would be no axial force both on and off the z-axis.

Equation (22) has been verified experimentally at Caltech by Kumar and Bellan [20] who measured the time of flight of an MHD jet as a function of both current I and mass density $\rho$. Zhai et al [3] solved the MHD equations numerically for plasma parameters corresponding to the Caltech experiment and observed jet behavior consistent with the model presented here. The magnetically generalized Bernoulli relation given by equation (17) has been verified in the 3D numerical simulation by Zhai et al [3] who observed that the local values at $r=a(z)$ of $\rho {U}_{z}^{2}+{B}_{\phi }^{2}/(2{\mu }_{0})\$ remain constant with respect to change in z. The solutions in Zhai et al [3] were shown to scale to astrophysical jet parameters.

The jet velocity increases with distance from the electrode so long as a(z) is a weak function of z as discussed above. However, at the jet tip a(z) is no longer a weak function of z because at the tip the poloidal magnetic field originating from the inner electrode bends over and goes back to intercept the outer electrode. Thus, the acceleration argument leading to equation (17) is no longer applicable. The bending over of the poloidal magnetic field results in a decelerating force that slows down the jet at the jet tip. This slowing down at the tip can be understood intuitively by recalling that the magnetic force can be decomposed as

$$\begin{eqnarray}&&{\bf{J}}\times {\bf{B}}=-{{\rm{\nabla }}}_{\perp }\left(\displaystyle \frac{{B}^{2}}{2{\mu }_{0}}\right)-\displaystyle \frac{{B}^{2}}{{\mu }_{0}}\displaystyle \frac{\hat{R}}{R}.\end{eqnarray} \tag{ 24 }$$

The first term, $-{{\rm{\nabla }}}_{\perp }{B}^{2}/2{\mu }_{0}$, behaves like a pressure in the direction perpendicular to the magnetic field while the second term, $-{\mu }_{0}^{-1}{B}^{2}\hat{R}/R,$ points in the direction opposite to the radius of curvature vector ${\bf{R}}$ of the magnetic field and has magnitude inversely proportional to the radius of curvature. Suppose that the jet is extremely long and collimated in which case the poloidal field has a near hair-pin shape at the jet tip. Thus, a typical poloidal field line is approximately straight going from the inner electrode up to the jet tip and then reverses direction to return to the outer electrode. This reversal of poloidal field direction near the jet tip means there is a sharp 180° turn of the poloidal field at the jet tip. The radius of curvature of the poloidal field will thus be small at the jet tip but large elsewhere. Since the associated force scales inversely with ${R}^{-1}$ and is directed in the $-\hat{R}$ direction, there will be a strong retarding force at the jet tip but not elsewhere.

An equivalent way of looking at this development of a retarding force at the jet tip is to revisit the argument that $\sim {\mu }_{0}^{-1}[{\rm{\nabla }}\cdot ({r}^{-2}{\rm{\nabla }}\psi )]\partial \psi /\partial z$ can be dropped compared to $-{r}^{-2}{\mu }_{0}I\partial I/\partial z$ in the axial force specified by the right hand side of equation (9c). Let us assume temporarily that $a={a}_{0}\exp ({\int }_{0}^{z}\kappa ({z}^{{\prime} }){{\rm{d}}{z}}^{{\prime} })$ where $\kappa (z)$ is small and nearly independent of z for most of the jet length but then suddenly becomes large at the jet tip so as to provide the flaring out of the poloidal flux at the jet tip (the mushroom cap). With this assumption we express

$$\begin{eqnarray}\displaystyle \frac{1}{{\mu }_{0}}{\rm{\nabla }}\cdot \left(\displaystyle \frac{1}{{r}^{2}}{\rm{\nabla }}\psi \right) & = & \displaystyle \frac{{\psi }_{0}}{{\mu }_{0}{r}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {z}^{2}}{\left(\displaystyle \frac{1}{a(z)}\right)}^{2}\\ & = & \displaystyle \frac{{\psi }_{0}}{{\mu }_{0}{a}_{0}^{2}{r}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {z}^{2}}\exp \left(-2{\displaystyle \int }_{0}^{z}\kappa ({z}^{{\prime} }){{\rm{d}}{z}}^{{\prime} }\right)\\ & \sim & 4\displaystyle \frac{{\kappa }^{2}}{{\mu }_{0}{r}^{2}}\psi +...\end{eqnarray} \tag{ 25 }$$

so

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\mu }_{0}}\left[{\rm{\nabla }}\cdot \left(\displaystyle \frac{1}{{r}^{2}}{\rm{\nabla }}\psi \right)\right]\displaystyle \frac{\partial \psi }{\partial z}\ \sim 4\displaystyle \frac{{\kappa }^{2}}{{\mu }_{0}{r}^{2}}\psi \displaystyle \frac{\partial \psi }{\partial z}.\end{eqnarray} \tag{ 26 }$$

By comparison, the term providing acceleration scales as

$$\begin{eqnarray}&&-\displaystyle \frac{{\mu }_{0}I}{{r}^{2}}\displaystyle \frac{\partial I}{\partial z}=-\displaystyle \frac{{\lambda }^{2}\psi }{{\mu }_{0}{r}^{2}}\displaystyle \frac{\partial \psi }{\partial z}\end{eqnarray} \tag{ 27 }$$

so the retardation given by equation (26) will overcome the acceleration given by equation (27) when the flaring out of the poloidal flux is such that $\kappa (z)\gt \lambda /2$.

Being concentrated in the vicinity of the jet tip, the retarding force slows down the jet in the vicinity of the jet tip. An observer sitting in the jet frame near the tip would consequently see a converging axial velocity, i.e., a negative ${\rm{\nabla }}\cdot {\bf{U}}$. The equation of continuity in the jet frame expressed as $\partial \rho /\partial t+{\bf{U}}\cdot {\boldsymbol{\nabla }}\rho =-\rho {\rm{\nabla }}\cdot {\bf{U}}$ implies that the plasma density seen by the observer in the jet frame increases where this local slowing occurs, much like traffic on a highway bunches up when a group of fast cars approaches a group of slower cars traveling in front of the fast cars. The axial compression of the jet plasma compresses the toroidal magnetic flux embedded in the jet plasma and this compression corresponds to increasing ${B}_{\phi }$ because ${B}_{\phi }$ is the density of toroidal flux. Since ${B}_{\phi }$ is what causes the pinch force, the bunching up of the fast jet approaching the slower jet at the tip increases pinching and so the jet becomes collimated. The radial pinching by axial deceleration can be seen quantitatively by first dotting the induction equation with ${\rm{\nabla }}\phi$ to obtain

$$\begin{eqnarray}&&{\rm{\nabla }}\phi \cdot \displaystyle \frac{\partial {\bf{B}}}{\partial t}={\rm{\nabla }}\phi \cdot {\rm{\nabla }}\times ({\bf{U}}\times {\bf{B}})\end{eqnarray} \tag{ 28 }$$

and then re-arranging equation (28) as

$$\begin{eqnarray}\displaystyle \frac{\partial }{\partial t}\left(\displaystyle \frac{{B}_{\phi }}{r}\right) & = & {\rm{\nabla }}\cdot [({\bf{U}}\times {\bf{B}})\times {\rm{\nabla }}\phi ]\\ & = & {\rm{\nabla }}\cdot \left[{\bf{B}}\displaystyle \frac{{U}_{\phi }}{r}-{\bf{U}}\displaystyle \frac{{B}_{\phi }}{r}\right]\\ & = & {\bf{B}}\cdot {\rm{\nabla }}\displaystyle \frac{{U}_{\phi }}{r}-{\bf{U}}\cdot {\rm{\nabla }}\displaystyle \frac{{B}_{\phi }}{r}-\displaystyle \frac{{B}_{\phi }}{r}{\rm{\nabla }}\cdot {\bf{U}}.\end{eqnarray} \tag{ 29 }$$

Equation (29) is then recast to be in the jet frame as

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\left(\displaystyle \frac{{B}_{\phi }}{r}\right)={\bf{B}}\cdot {\rm{\nabla }}\displaystyle \frac{{U}_{\phi }}{r}-\displaystyle \frac{{B}_{\phi }}{r}{\rm{\nabla }}\cdot {\bf{U}}\ ,\end{eqnarray} \tag{ 30 }$$

where d/dt denotes the convective derivative $\partial /\partial t+{\bf{U}}\cdot {\boldsymbol{\nabla }}$, i.e., the time derivative seen by an observer in the jet frame. Since ${U}_{\phi }=0$ was stipulated, equation (30) reduces to

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\left(\displaystyle \frac{{B}_{\phi }}{r}\right)=-\displaystyle \frac{{B}_{\phi }}{r}{\rm{\nabla }}\cdot {\bf{U}}.\end{eqnarray} \tag{ 31 }$$

However, the continuity equation can be written as

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\rho }{{\rm{d}}{t}}+\rho {\rm{\nabla }}\cdot {\bf{U}}=0\end{eqnarray} \tag{ 32 }$$

so elimination of ${\rm{\nabla }}\cdot {\bf{U}}$ from equations (31) and (32) gives

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{B}_{\phi }}{r}\right)}^{-1}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\left(\displaystyle \frac{{B}_{\phi }}{r}\right)=\displaystyle \frac{1}{\rho \ }\displaystyle \frac{{\rm{d}}\rho }{{\rm{d}}{t}}\end{eqnarray} \tag{ 33 }$$

or

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{\rm{ln}}\left(\displaystyle \frac{{B}_{\phi }}{r}\right)=\displaystyle \frac{{\rm{d}}\ }{{\rm{d}}{t}}{\rm{ln}}\rho \end{eqnarray} \tag{ 34 }$$

which leads to

$$\begin{eqnarray}&&\displaystyle \frac{{B}_{\phi }}{\rho r}={\rm{const}}.\ \mathrm{in}\ \mathrm{the}\ \mathrm{jet}\ \mathrm{frame}.\end{eqnarray} \tag{ 35 }$$

Thus axial compression increases ρ at a given r which causes a corresponding increase in ${B}_{\phi }\$ and so an increase in pinching. Since the total current $I=2\pi {{aB}}_{\phi }/{\mu }_{0}$ is a constant, the current channel radius a must decrease in the jet frame as ${B}_{\phi }$ increases. This reduction of a at the jet tip collimates the jet. The jet tip thus self-collimates and this self-collimation acts like a zipper continuously squeezing together poloidal field lines as the jet propagates. Because the flux tube is current-carrying, the magnetic field is helical and can be represented as ${\bf{B}}={B}_{\phi }\hat{\phi }+{B}_{z}\hat{z}$. The jet flow is away from the region where the flux tube diameter is small and the jet flow direction remains the same when the current polarity is reversed. This is because the jet mechanism ultimately depends on the axial gradient of ${B}_{\phi }^{2}$ which clearly is independent of the sign of ${B}_{\phi }$ and hence independent of the sign of the current.

The jet has a constant poloidal flux, but this flux is being axially distended since it is frozen into the jet. However, the toroidal flux in the jet is increasing because while ${B}_{\phi }$ is constant in the jet, the toroidal flux scales as the length of the jet. Thus, the linked toroidal and poloidal flux, i.e., the helicity, is increasing with time. The rate of helicity injection is proportional to the radial voltage drop at the electrodes driving the jet since this voltage corresponds to the rate of injection of toroidal flux. The jet requires a radial electric field to drive the poloidal current I and create the ever-increasing toroidal flux resulting from the lengthening of the jet. This radial electric field has been provided by a capacitor in laboratory experiments [17, 18, 22, 26] but obviously must be provided by other means in the situation of an accretion disk.

I and ψ are defined for both positive and negative z and, as sketched in figure 2, $\ I$ is assumed to be an odd function of z while ψ is assumed to be an even function of $z$. Both I and ψ depend on r and both vanish at infinity. The disk lies in the z = 0 plane and is thin, weakly ionized, and in a Kepler orbit around a star. ${B}_{\phi }$, ${B}_{r,}$ and J_z are antisymmetric with respect to z since ${B}_{\phi }\sim I,$ ${B}_{r}\sim \partial \psi /\partial z$ and ${J}_{z}\sim \partial I/\partial r$ whereas ${B}_{z}\sim \partial \psi /\partial r$, ${J}_{r}\sim \partial I/\partial z$ and ${J}_{\phi }$ are symmetric with respect to $z$. The axial component of the magnetic force is thus antisymmetric with respect to z since

$$\begin{eqnarray}&&{({\bf{J}}\times {\bf{B}})}_{z}=\mathop{\underbrace{{J}_{r}}}\limits_{{\rm{symmetric}}}\mathop{\underbrace{{B}_{\phi }}}\limits_{{\rm{antisymmetric}}}-\mathop{\underbrace{{J}_{\phi }}}\limits_{{\rm{symmetric}}}\mathop{\underbrace{{B}_{r}}}\limits_{{\rm{antisymmetric}}},\end{eqnarray} \tag{ 36 }$$

whereas the radial component of the magnetic force is symmetric with respect to z since

$$\begin{eqnarray}&&{({\bf{J}}\times {\bf{B}})}_{r}=\mathop{\underbrace{{J}_{\phi }}}\limits_{{\rm{symmetric}}}\mathop{\underbrace{{B}_{z}}}\limits_{{\rm{symmetric}}}-\mathop{\underbrace{{J}_{z}}}\limits_{{\rm{antisymmetric}}}\mathop{\underbrace{{B}_{\phi }}}\limits_{{\rm{antisymmetric}}}.\end{eqnarray} \tag{ 37 }$$

Because ${B}_{\phi }$ is antisymmetric with respect to z while B_z is symmetric, the $z\gt 0$ and $z\lt 0$ regions have opposite helicity magnetic fields so the total helicity of the entire system remains zero. The radial electric field E_r and the radial electric current J_r are oppositely directed in the z = 0 plane which indicates there is a power source in the z = 0 plane. This negative value of ${\bf{E}}\cdot {\bf{J}}$ in the disk (corresponding to the existence of a power source) is missing from models invoking anomalous resistivity (anomalously large magnetic diffusion) since increasing the resistivity corresponds to creating a sink for power and not a source. Since the disk is presumed to drive the jets, the disk should be a power source and the jets a power sink.

The accretion disk is assumed to be in a Kepler orbit, i.e., at each radius the outward centrifugal force $\rho {U}_{\phi }^{2}/r$ is balanced by the inward gravitational force $\rho {MG}/{r}^{2}$ where M is the mass of the star and G is the gravitational constant so ${U}_{\phi }=\sqrt{{MG}/r}$, the Kepler orbital velocity. Previous astrophysical jet models have presumed that the accretion disk acts as a homopolar generator [27] in the presence of some pre-existing constant poloidal magnetic field, i.e., the B_z field threading the disk plane. According to this presumption, the Kepler rotation of the disk cuts the B_z field lines to create a radial electric field as prescribed by the ideal MHD Ohm's law equation (3c). The radial component of equation (3c) is

$$\begin{eqnarray}&&{E}_{r}+{U}_{\phi }{B}_{z}-{U}_{z}{B}_{\phi }=0\end{eqnarray} \tag{ 38 }$$

and because the anti-symmetry of the bidirectional jets mandates ${U}_{z}=0$ in the z = 0 plane, equation (38) reduces to ${E}_{r}=-{U}_{\phi }{B}_{z}$ where ${U}_{\phi }$ is the disk's Kepler rotational velocity. Despite the intuitive appeal of this argument, it cannot be correct because if equation (3c) is to be invoked, then all components of equation (3c) must be satisfied, not just the radial component. The failure of equation (3c) becomes evident upon consideration of its azimuthal component, namely

$$\begin{eqnarray}&&{E}_{\phi }+{U}_{z}{B}_{r}-{U}_{r}{B}_{z}=0.\end{eqnarray} \tag{ 39 }$$

The B_z magnetic field threading the z = 0 plane is assumed to be both constant in time and axisymmetric. The azimuthal electric field is ${E}_{\phi }=-{r}^{-1}\partial V/\partial \phi -\partial {A}_{\phi }/\partial t$. It is seen that ${E}_{\phi }$ must vanish in the z = 0 plane since axisymmetry implies $\partial V/\partial \phi =0$ while B_z being constant in time implies $\partial {A}_{\phi }/\partial t=0$. The anti-symmetry with respect to z of the jet axial velocity implies ${U}_{z}=0$ in the z = 0 plane. Because the jets are a bidirectional outflow of matter from the disk plane and because the configuration is in steady-state in the z = 0 plane, there must be a continuous radially inward flow of matter in the z = 0 plane to supply the matter continuously flowing out of the z = 0 plane in the bidirectional jets. This implies that U_r must be both finite and radially inwards in the z = 0 plane; i.e., there must be accretion. Because B_z is finite in the z = 0 plane, the term $-{U}_{r}{B}_{z}$ in equation (39) must be finite and yet no other term exists in equation (39) to balance this term (recall that ${E}_{\phi }=0$ and ${U}_{z}=0$ in the z = 0 plane). Thus, the disk cannot be governed by equation (3c) and so ideal MHD is an inadequate model for the accretion disk in contrast to ideal MHD being adequate for the jet region external to the accretion disk.

We claim that a plasma cannot simultaneously satisfy the ideal MHD equations and be in a Kepler orbit [28]. To prove this claim we examine the distinction between Kepler and cyclotron motion in detail. In order to present the essential idea in the clearest possible fashion, we first restrict consideration to motion in the z = 0 plane and also assume the magnetic field is constant, uniform, and normal to the z = 0 plane. The magnetic field is then given by ${\bf{B}}={B}_{z}\hat{z}$ with an associated azimuthal vector potential ${A}_{\phi }={B}_{z}r/2$ and an associated poloidal flux $\psi ={B}_{z}\pi {r}^{2}$. We consider a particle having mass m and charge $q\$ in the gravitational field of a central object of mass $M$. The respective Lagrangian and Hamiltonian for the particle in the combined gravitational and magnetic field are

$$\begin{eqnarray}&&{ \mathcal L }=\displaystyle \frac{1}{2}m({v}_{r}^{2}+{r}^{2}{\dot{\phi }}^{2})+{qr}\dot{\phi }{A}_{\phi }+\displaystyle \frac{{mMG}}{r}\end{eqnarray} \tag{ 40 }$$

and

$$\begin{eqnarray}&&H=\displaystyle \frac{1}{2}m({v}_{r}^{2}+{r}^{2}{\dot{\phi }}^{2})-\displaystyle \frac{{mMG}}{r}.\end{eqnarray} \tag{ 41 }$$

The canonical angular momentum defined as ${P}_{\phi }=\partial { \mathcal L }/\partial \dot{\phi }$ is

$$\begin{eqnarray}{P}_{\phi } & = & {{mr}}^{2}\dot{\phi }+{{qrA}}_{\phi }(r)\\ & = & {{mr}}^{2}\dot{\phi }+\displaystyle \frac{q}{2\pi }\psi (r).\end{eqnarray} \tag{ 42 }$$

Lagrangian mechanics shows that because $\partial { \mathcal L }/\partial \phi =0,$ ${P}_{\phi }$ is a constant of the motion. An important consequence of the term containing $q\psi (r)$ in equation (42) is that a charged particle behaves very differently from a neutral particle. Equation (42) can be rearranged to give the angular velocity

$$\begin{eqnarray}&&\dot{\phi }=\displaystyle \frac{{P}_{\phi }-\tfrac{q}{2\pi }\psi (r)}{{{mr}}^{2}}.\end{eqnarray} \tag{ 43 }$$

On substituting equation (43) into (41) the Hamiltonian can be expressed as

$$\begin{eqnarray}&&H=\displaystyle \frac{1}{2}{{mv}}_{r}^{2}+\chi (r),\end{eqnarray} \tag{ 44 }$$

where the effective potential $\chi (r)$ can be decomposed into two terms as

$$\begin{eqnarray}&&\chi (r)={\chi }_{\phi }(r)+{\chi }_{{\rm{g}}}(r).\end{eqnarray} \tag{ 45 }$$

The first term is the azimuthal kinetic energy ${{mv}}_{\phi }^{2}/2$ expressed as [28, 29]

$$\begin{eqnarray}&&{\chi }_{\phi }(r)=\displaystyle \frac{{\left({P}_{\phi }-\tfrac{q}{2\pi }\psi (r)\right)}^{2}}{2{{mr}}^{2}}\end{eqnarray} \tag{ 46 }$$

while the second term

$$\begin{eqnarray}&&{\chi }_{{\rm{g}}}(r)=-\displaystyle \frac{{mMG}}{r}\end{eqnarray} \tag{ 47 }$$

is the gravitational potential energy. If the particle is neutral so q = 0 then χ reverts to the classic Kepler form $\chi (r)={P}_{\phi }^{2}/(2{{mr}}^{2})-{mMG}/r$ which leads to circular Kepler orbits when r is at the minimum of $\chi (r)$. Using $\psi \sim {B}_{z}\pi {r}^{2}$, comparison of the orders of magnitudes of ${\chi }_{\phi }$ and ${\chi }_{{\rm{g}}}$ gives

$$\begin{eqnarray}&&\displaystyle \frac{{\chi }_{\phi }}{{\chi }_{g}}\sim \displaystyle \frac{{q}^{2}{B}_{z}^{2}{r}^{3}}{{m}^{2}{MG}}=\displaystyle \frac{{\omega }_{{\rm{c}}}^{2}}{{\omega }_{{\rm{K}}}^{2}},\end{eqnarray} \tag{ 48 }$$

where ${\omega }_{{\rm{c}}}={{qB}}_{z}/m$ is the cyclotron frequency and ${\omega }_{{\rm{K}}}\,=\sqrt{{MG}/{r}^{3}}$ is the Kepler frequency. For typical astrophysical parameters, cyclotron frequencies are many orders of magnitude larger than the Kepler frequency so ${\chi }_{\phi }$ is orders of magnitude larger than ${\chi }_{{\rm{g}}}$. This means that gravity is of negligible importance for charged particles and so charged particles will not execute Kepler orbits. If a plasma is fully ionized, all particles are either electrons or ions, no particle executes a Kepler orbit [28], and so the plasma center of mass cannot execute a Kepler orbit. Hence, it is not possible for a fully ionized ideal MHD plasma to be in a Kepler orbit. However, if the magnitude of ${\omega }_{{\rm{c}}}$ were somehow reduced to become comparable to ${\omega }_{{\rm{K}}}$, the two effective potentials would then have associated forces of similar magnitude in which case the motion would differ substantively from either strict cyclotron or strict Kepler motion. This reduction of ${\omega }_{{\rm{c}}}^{2}/{\omega }_{{\rm{K}}}^{2}$ could be accomplished either by having a very weak magnetic field B_z or by having a very small charge to mass ratio $q/m$.

Reduction of ${\omega }_{{\rm{c}}}^{2}/{\omega }_{{\rm{K}}}^{2}$ by having charged dust grains was considered in [30] and while this appears physically possible, it seems unlikely to be the typical situation. Instead, the effect of weak ionization appears to be a more likely means by which this reduction could occur because accretion disks rotating around new-born stars are very weakly ionized [7, 14] with fractional ionization ratios in the range $\alpha \sim {10}^{-13}\mbox{--}{10}^{-8}$. The neutral particle density is quite high, nominally ${n}_{n}\sim {10}^{16}$ m⁻³. Because of the weak ionization, to a first approximation one can ignore the charged particles and simply presume that the neutrals are in Kepler orbits around the star. Also, because there are ${\alpha }^{-1}$ more neutrals than charged particles, charged particles collide mainly with neutrals rather than with other charged particles. Assuming a nominal temperature of 100 K, the thermal velocity of neutrals is ${v}_{{\rm{T}}{n}}\sim {10}^{3}$ m s⁻¹ and so assuming a nominal atomic cross section $\sigma =3\times {10}^{-20}$ m², the nominal neutral–neutral collision frequency will be ${\nu }_{{nn}}=\sigma {n}_{n}{v}_{{\rm{T}}{n}}=\ 0.3$ s⁻¹ and collisions of an ion with a neutral will have a similar frequency. It is assumed that ion-neutral collisions are so frequent that ions collide with neutrals and scatter before completing a cyclotron orbit; this means that ions can be considered unmagnetized. On the other hand, because the electron cyclotron frequency is ${m}_{i}/{m}_{e}$ times higher than the ion cyclotron frequency whereas the electron thermal velocity is only $\sqrt{{m}_{i}/{m}_{e}}$ times larger than the ion thermal velocity for equal electron and ion temperatures, regimes exist where the electrons complete a cyclotron orbit before being scattered by a collision. Thus, we are considering a regime [8, 31] where the electrons are magnetized but the ions are not. By 'magnetized' we mean that the particle completes a cyclotron orbit.

We now argue that despite the ions being unmagnetized, the combination of ions and neutrals behaves in the manner of a fictitious magnetized particle having a charge to mass ratio smaller than that of an ion by the fractional ionization $\alpha$. To see this, consider a hypothetical situation where a transient electromagnetic force imparts a certain momentum to the ions during an interval shorter than the ion-neutral collision time. After acquiring this momentum, the ions collide with the neutrals and share this momentum with the neutrals. Because of this sharing, the combination of the ions and neutrals have the momentum the ions received from the electromagnetic force and the total mass of the ions and the neutrals. This shared momentum is identical to what would have been gained by a fictitious particle having the ion charge and a mass equal to the sum of the masses of the ions and neutrals under consideration and so a charge to mass ratio $\alpha {q}_{i}/{m}_{i}$ where α is the fractional ionization. Thus, the center of mass of the combination of ions and neutrals behaves like a fictitious 'metaparticle' having cyclotron frequency ${\omega }_{{\rm{c}}}=\alpha {q}_{i}/{m}_{i};$ this reduction of effective cyclotron frequency was previously noted by Pandey and Wardle [14] who showed this would greatly lower the frequency at which whistler wave physics occurs since the cross-over from MHD physics to whistler physics occurs in the vicinity of the ion cyclotron frequency. Since $\alpha \sim {10}^{-13}\mbox{--}{10}^{-8}$ and ion cyclotron frequencies are typically 8–13 orders of magnitude larger than Kepler frequencies, metaparticles with a cyclotron frequency of the order of the Kepler frequency should exist in a typical accretion disk. Ionization is typically stratified [7] because α depends on ionization via x-ray or ultra-violet radiation, and so typically α varies from being essentially zero in the disk equatorial plane (as this plane is shielded by the disk particles from x-rays and UV) to being near unity at the interface between the disk face and the low-density, fully ionized plasma external to the disk. Thus, there will always be a disk stratum where α is such that the metaparticle cyclotron frequency is of the order of the Kepler frequency.

The time required for the neutrals to become collisionally attached to an ion can be estimated by considering a single ion and its associated ${\alpha }^{-1}$ neutrals (e.g. one ion and 10¹² neutrals) and supposing that the ion receives some momentum from an electromagnetic field. The ion can be considered to constitute a delta-function momentum input which will collisionally diffuse through the surrounding ${\alpha }^{-1}$ neutrals. The volume occupied by the ${\alpha }^{-1}$ neutrals is $V={(\alpha {n}_{n})}^{-1}$ so the characteristic linear dimension of the clump of the ${\alpha }^{-1}$ neutrals is $l={V}^{1/3}={(\alpha {n}_{n})}^{-1/3}$. The momentum diffusion coefficient is $D\sim {v}_{{\rm{T}}}^{2}/{\nu }_{{nn}}={v}_{{\rm{T}}}/(\sigma {n}_{n})$. The diffusion equation gives the time t_d to diffuse a length l as ${l}^{2}\simeq 4{{Dt}}_{{\rm{d}}}$ so

$$\begin{eqnarray}&&{t}_{{\rm{d}}}=\displaystyle \frac{{l}^{2}}{4D}=\displaystyle \frac{\sigma {n}_{n}^{1/3}{\alpha }^{-2/3}}{4{v}_{{\rm{T}}}}\end{eqnarray} \tag{ 49 }$$

which indicates times of the order of nanoseconds or less for the electromagnetic momentum gained by an ion to be shared via collisional diffusion with its associated ${\alpha }^{-1}$ neutrals assuming $\alpha \sim {10}^{-13}\mbox{--}{10}^{-8}$. Thus, it is very reasonable to consider an ion and its associated ${\alpha }^{-1}$ neutrals to be collisionally attached to each other.

The metaparticle equation of motion can be derived in an alternative, equivalent, and more formal way by consideration of the MHD equation of motion for a weakly ionized plasma in a gravitational field, namely

$$\begin{eqnarray}&&\rho \displaystyle \frac{{\rm{d}}{\bf{U}}}{{\rm{d}}{t}}={\bf{J}}\times {\bf{B}}-{\rm{\nabla }}({P}_{i}+{P}_{e}+{P}_{n})-\rho {\bf{g}}\end{eqnarray} \tag{ 50 }$$

and the Hall Ohm's law expressed in the form of the electron fluid equation of motion,

$$\begin{eqnarray}&&0={n}_{e}{q}_{e}({\bf{E}}+{{\bf{u}}}_{e}\times {\bf{B}})-{\rm{\nabla }}{P}_{e}.\end{eqnarray} \tag{ 51 }$$

Electron collisions have been dropped from equation (51) since electrons are assumed magnetized so ${\omega }_{{\rm{c}}{e}}$ exceeds the electron-neutral collision frequency. The gravitational force on electrons has also been dropped from this equation because the electron cyclotron frequency exceeds the Kepler frequency by many orders of magnitude. Using ${\bf{J}}={n}_{i}{q}_{i}{{\bf{u}}}_{i}+{n}_{e}{q}_{e}{{\bf{u}}}_{e}$, subtraction of equation (51) from (50) yields

$$\begin{eqnarray}&&\rho \displaystyle \frac{{\rm{d}}{\bf{U}}}{{\rm{d}}{t}}=-{n}_{e}{q}_{e}{\bf{E}}+{n}_{i}{q}_{i}{{\bf{u}}}_{i}\times {\bf{B}}-{\rm{\nabla }}({P}_{i}+{P}_{n})-\rho {\bf{g}}.\end{eqnarray} \tag{ 52 }$$

Since the plasma is weakly ionized, the overall mass density is essentially the neutral mass density, i.e., $\rho ={n}_{n}{m}_{n}$ where m_n is the mass of a neutral particle and $\ {n}_{n}$ is the neutral particle number density. Dividing equation (52) by ρ and invoking the quasi-neutrality condition ${n}_{e}{q}_{e}=-{n}_{i}{q}_{i}$ yields

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\bf{U}}}{{\rm{d}}{t}}=\alpha \displaystyle \frac{{q}_{i}}{{m}_{n}}\ ({\bf{E}}+{{\bf{u}}}_{i}\times {\bf{B}})-\displaystyle \frac{{\rm{\nabla }}({P}_{i}+{P}_{n})}{{n}_{n}{m}_{n}}-{\bf{g}}.\end{eqnarray} \tag{ 53 }$$

Because the ions collide frequently with the neutrals, the mean ion velocity ${{\bf{u}}}_{i}$ is nearly the center of mass velocity ${\bf{U}}$ which is essentially the neutral center of mass velocity. The mean ion and center of mass velocities are presumed to be of order the Kepler velocity which at 1 a.u. is 3 × 10⁴ m s⁻¹. A hydrogen gas with nominal temperature of 100 K has a thermal velocity of the order of 10³ m s⁻¹ which is significantly lower than the Kepler velocity. Also, the direction of rotation is assumed without loss of generality to be in the positive ϕ direction; or equivalently, the direction of Kepler rotation defines the positive ϕ direction. Finally, if the pressure term can be neglected on the basis that the random thermal velocities are small compared to the Kepler-like center of mass velocity, equation (53) reduces to

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\bf{U}}}{{\rm{d}}{t}}=\alpha \displaystyle \frac{{q}_{i}}{{m}_{n}}\ ({\bf{E}}+{\bf{U}}\times {\bf{B}})-{\bf{g}}\end{eqnarray} \tag{ 54 }$$

which is the equation of motion for a metaparticle with charge to mass ratio $\alpha {q}_{i}/{m}_{n}$. Assuming the neutral mass is similar to the ion mass, the effective cyclotron frequency is thus reduced from the ion cyclotron frequency by $\alpha ,$ i.e., by ${10}^{-13}\mbox{--}{10}^{-8}$. Ion-neutral collisions do not exist in equation (54) because the momentum ions lose from colliding with neutrals is precisely the momentum neutrals gain from these collisions and equation (54) describes the evolution of the sum of the ion and neutral momenta. Thus, ion-neutral and neutral-ion collisions have no effect on the metaparticle momentum. Instead, these collisions are the means by which the ions and neutrals are bound together to form the metaparticle. The metaparticle point of view is thus a type of two-fluid system where instead of having the two fluids being ions and electrons, now the two fluids are metaparticles and electrons.

A metaparticle is thus a clump of particles having effective charge to mass ratio $\alpha {q}_{i}/{m}_{i}$. Just as ion motion for phenomena slow compared to ${\omega }_{{\rm{c}}{i}}$ can be described by the guiding center approximation (i.e., ${\bf{E}}\times {\bf{B}}$ drifts, polarization drifts, curvature drifts, and grad B drifts), the metaparticle motion can similarly be described by the guiding center approximation for phenomena slow compared to $\alpha {\omega }_{{\rm{c}}{i}}$. However, if the Kepler frequency ${\omega }_{{\rm{K}}}$ is comparable to the metaparticle cyclotron frequency $\alpha {\omega }_{{\rm{c}}{i}}$, then the guiding center approximation cannot be used for metaparticles and a qualitatively different behavior occurs. An essential result of the guiding center approximation is that particles remain on poloidal magnetic flux surfaces to the extent that the ${{mr}}^{2}\dot{\phi }$ term in equation (42) is small compared to the $q\psi /2\pi$ term. However if ${\omega }_{{\rm{K}}}\simeq \alpha {\omega }_{{\rm{c}}{i}}$, the two terms in equation (42) become comparable, the guiding center approximation fails, and metaparticles can move across flux surfaces. In contrast to metaparticles, electrons remain governed by the guiding center approximation and so cannot move across magnetic flux surfaces. This distinction between metaparticle motion and electron motion leads to the possibility of a steady radial electric current composed of metaparticles, but not electrons, flowing across magnetic flux surfaces. The metaparticle constitutes the carrier of the radial current.

Equation (54) is the equation of motion for a charged particle with mass m_n and charge $\alpha {q}_{i}$ in a gravitational field and a magnetic field. The discussion can now be generalized by removing the four previously imposed restrictions that (i) motion is confined the z = 0 plane, (ii) the magnetic field is constant, (iii) the magnetic field is uniform, and (iv) there is no electric field. With the removal of these restrictions the metaparticle Hamiltonian H generalizes to

$$\begin{eqnarray}&&\displaystyle \frac{H}{{m}_{n}}=\displaystyle \frac{1}{2}({\dot{r}}^{2}+{r}^{2}{\dot{\phi }}^{2}+{\dot{z}}^{2})-\displaystyle \frac{{MG}}{\sqrt{{r}^{2}+{z}^{2}}}+\displaystyle \frac{\alpha {q}_{i}}{{m}_{n}}V(r,z,t),\end{eqnarray} \tag{ 55 }$$

where $V(r,z,t)$ is an electrostatic potential. Because the system remains axisymmetric, the canonical angular momentum is still a constant of the motion, i.e.,

$$\begin{eqnarray}&&{P}_{\phi }={m}_{n}{r}^{2}\dot{\phi }+\displaystyle \frac{\alpha {q}_{i}}{2\pi }\psi (r,z,t)={\rm{const}}.\end{eqnarray} \tag{ 56 }$$

On substituting for $\dot{\phi }$ using equation (56) the Hamiltonian assumes the form

$$\begin{eqnarray}&&\displaystyle \frac{H}{{m}_{n}}=\displaystyle \frac{1}{2}({\dot{r}}^{2}+{\dot{z}}^{2})+\chi (r,z,t),\end{eqnarray} \tag{ 57 }$$

where the effective potential χ now is

$$\begin{eqnarray}\chi (r,z,t) & = & \displaystyle \frac{{\left(\tfrac{{P}_{\phi }}{{m}_{n}}-\tfrac{1}{2\pi }\tfrac{\alpha {q}_{i}}{{m}_{n}}\psi (r,z,t)\right)}^{2}}{2{r}^{2}}\\ & & -\displaystyle \frac{{MG}}{\sqrt{{r}^{2}+{z}^{2}}}+\displaystyle \frac{\alpha {q}_{i}}{{m}_{n}}V(r,z,t).\end{eqnarray} \tag{ 58 }$$

Now consider the special case where ${P}_{\phi }=0$ and assume temporarily that V = 0 so the effective potential becomes

$$\begin{eqnarray}&&\chi (r,z,t)=\displaystyle \frac{{\alpha }^{2}{q}_{i}^{2}}{{m}_{n}^{2}}\displaystyle \frac{{(\psi (r,z,t))}^{2}}{8{\pi }^{2}{r}^{2}}-\displaystyle \frac{{MG}}{\sqrt{{r}^{2}+{z}^{2}}}.\end{eqnarray} \tag{ 59 }$$

Since $\ {B}_{z}$ is approximately uniform at small $r,$ the poloidal flux for small r is

$$\begin{eqnarray}&&\psi (r,x,t)=\pi {B}_{z}(0,z,t){r}^{2}\end{eqnarray} \tag{ 60 }$$

in which case

$$\begin{eqnarray}&&\chi (r,z,t)=\displaystyle \frac{{\alpha }^{2}{q}_{i}^{2}}{{m}_{n}^{2}}\displaystyle \frac{{B}_{z}^{2}(0,z,t)}{8}{r}^{2}-\displaystyle \frac{{MG}}{\sqrt{{r}^{2}+{z}^{2}}}.\end{eqnarray} \tag{ 61 }$$

The radial force on the particle is ${\bf{F}}=-\partial \chi /\partial r$ and this force is always radially inwards because the radial partial derivative of both terms on the right hand side of equation (61) is positive. Thus, a particle is subjected to an inward radial force at any $z$. If the particle reaches r = 0 at finite $z,$ the gravitational term provides an axial force towards the midplane z = 0. Thus, a particle will always experience a combination of forces directing it towards $r=0,$ z = 0, i.e., towards the central mass. The ${P}_{\phi }=0$ assumption can be used in equation (56) to give

$$\begin{eqnarray}&&\dot{\phi }=-\displaystyle \frac{1}{2\pi {r}^{2}}\displaystyle \frac{\alpha {q}_{i}}{{m}_{n}}\psi (r,z,t)\end{eqnarray} \tag{ 62 }$$

so using equation (60)

$$\begin{eqnarray}&&\dot{\phi }=-\displaystyle \frac{1}{2}\displaystyle \frac{\alpha {q}_{i}}{{m}_{n}}{B}_{z}(0,z,t)\end{eqnarray} \tag{ 63 }$$

i.e., $\dot{\phi }$ is minus half the metaparticle cyclotron frequency. Suppose a metaparticle was initially orbiting at a radius r₀ with a Kepler frequency ${\omega }_{{\rm{K}}}=\sqrt{{MG}/{r}_{0}^{3}}$ and that ${P}_{\phi }=0$ was satisfied at this initial radius. Because ${P}_{\phi }$ is a constant of the motion, ${P}_{\phi }$ will remain zero throughout the particle's entire motion. The combination of radially inward velocity and angular velocity $\dot{\phi }=-\alpha {\omega }_{{\rm{c}}{i}}/2$ means the ${P}_{\phi }=0$ metaparticle makes an inward spiral trajectory. This can be easily verified by solving equation (54) numerically with ${P}_{\phi }=0$ set as an initial condition. Note that B_z must be negative at small r since we have assumed that $\dot{\phi }$ is positive; this corresponds to the downward direction of the magnetic field at small r in figure 2.

Since metaparticles have positive charge and spiral inwards, their accumulation near $r=0,$ z = 0 implies a build-up of positive charge near the origin and so the electrostatic potential V will become finite and there will be a radially outward electric field ${E}_{r}$. Metaparticle accumulation then is the mechanism for producing the radial electric field discussed in section 3. The concentration near r = 0 of positive charge from the accumulated metaparticles cannot be neutralized by electron motion in the z = 0 plane because electrons cannot move across poloidal flux surfaces. This constraint on electron motion results from conservation of the electron canonical angular momentum ${P}_{\phi }^{e}={m}_{e}{r}^{2}\dot{\phi }\,-e\psi (r,z,t)/2\pi ;$ noting that m_e is negligible ${P}_{\phi }^{e}={\rm{const}}.$ implies that $\psi (r,z,t)$ is conserved for an electron. Electrons therefore cannot deviate from their initial poloidal flux surface and so cannot move in the z = 0 plane to neutralize the accumulated positive charge of the inward spiraled metaparticles. However, as sketched in figure 2, electrons can move along constant ψ contours outside of the z = 0 plane and so can make their way from large r in the disk plane to small r in the disk plane by moving along a constant ψ surface. This electron motion out of the disk plane along constant ψ surfaces corresponds to a clockwise electric current in the region above the z = 0 plane and a counter-clockwise electric current in the region below the z = 0 plane. At small r, this poloidal electric current is axially outwards from the z = 0 midplane above and below the midplane. These poloidal currents produce a toroidal magnetic field into the page above the z = 0 plane and out of the page below the z = 0 plane. The topology of the current is precisely what is required to drive bidirectional jets away from the midplane.

Also, a portion of the inward spiraling metaparticles can serve as the feedstock for the jet mass flux. Since ${U}_{\phi }=r\dot{\phi }$ and $\dot{\phi }$ is constant for the inward spiraling metaparticles, the metaparticles have vanishing ${U}_{\phi }$ at small radius and so the assumption that the jet starts with ${U}_{\phi }=0$ is seen to be reasonable.

When the Hall term is included, the ideal MHD Ohm's law generalizes to be the Hall Ohm's law

$$\begin{eqnarray}&&{\bf{E}}+{\bf{U}}\times {\bf{B}}-\displaystyle \frac{1}{{n}_{e}e}{\bf{J}}\times {\bf{B}}=0.\end{eqnarray} \tag{ 64 }$$

The addition of the Hall term $-{\bf{J}}\times {\bf{B}}/({n}_{e}e)$ does not correspond to adding 'new' physics to the model but instead corresponds to correcting an over-simplification inherent in ideal MHD. This over-simplification is the assumption in ideal MHD that the Hall term can be dropped because the difference between each component of the perpendicular electron and ion velocities is much less than the ion velocity component, i.e., ideal MHD is based on the assumption that

$$\begin{eqnarray}&&\left|\displaystyle \frac{{J}_{k}}{{n}_{e}e}\right|=| {u}_{{ik}}-{u}_{{ek}}| \ll | {u}_{{ik}}| \ \approx | {U}_{k}| \end{eqnarray} \tag{ 65 }$$

for each perpendicular component $k$. Specifically, for a magnetic field in the z direction, k would be the r or ϕ directions. The assumption prescribed by equation (65) is true if both the ion and electron perpendicular motions are dominantly ${\bf{E}}\times {\bf{B}}$ drifts. However the assumption prescribed by equation (65) is clearly false when there are inward spiraling metaparticles since ions, being constituents of the metaparticles, are spiraling in whereas the electrons, frozen to flux surfaces, are not. The entire radial electric current is thus produced by the radially inward ion velocity and so equation (65) is false because ${u}_{{er}}=0\$ causes equation (65) to imply $| {u}_{{ir}}| \ll | {u}_{{ir}}|$ which is clearly incorrect. An equivalent point of view is to consider the azimuthal component of equation (64), namely

$$\begin{eqnarray}&&-{U}_{r}{B}_{z}+\displaystyle \frac{1}{{n}_{e}e}{J}_{r}{B}_{z}=0;\end{eqnarray} \tag{ 66 }$$

on dividing by B_z this becomes

$$\begin{eqnarray}&&{J}_{r}={n}_{e}{{eU}}_{r}=\alpha {n}_{n}{{eU}}_{r}\end{eqnarray} \tag{ 67 }$$

which shows that the radial current comes from the metaparticles.

Because the Hall term in equation (66) is not contained in ideal MHD there is no term to balance the $-{U}_{r}{B}_{z}$ term in the azimuthal component. Inclusion of the Hall Ohm's law is effectively a tautology as it contains no physical information, but when the Hall term is omitted, the tautology is falsified. This is seen by noting that the Hall Ohm's law in the z = 0 plane is just the electron equation of motion for the situation where the electrons have no radial motion. The assumption in equation (65) which is the basis for ideal MHD, falsely assigns the electrons a radial motion comparable to the ion radial motion.

An essential requirement for any accretion model is a mechanism for removing angular momentum of accreting mass. Suppose for the moment that the accreting mass is neutral so the effective potential is

$$\begin{eqnarray}&&\chi (r)=\displaystyle \frac{{P}_{\phi }^{2}}{2{{mr}}^{2}}-\displaystyle \frac{{mMG}}{r}.\end{eqnarray} \tag{ 68 }$$

If angular momentum were not removed, then as $r\to 0$ the ${P}_{\phi }^{2}/2{{mr}}^{2}$ contribution to the effective potential energy would diverge making it energetically impossible for accreting neutral material to attain arbitrarily small r. Because torque is the rate of change of angular momentum, removing angular momentum of accreting matter at some radial location r means there must be a negative or 'braking' torque at this r. Magnetic forces can provide such a braking torque and the detailed mechanism can be understood in terms of conservation of canonical angular momentum or, equivalently, by consideration of the azimuthal component of the MHD equation of motion.

Consider an axisymmetric magnetic field and its associated vector potential ${{rA}}_{\phi }(r,z,t)=\psi (r,z,t)/2\pi \$ so conservation of canonical angular momentum is expressed as

$$\begin{eqnarray}&&{P}_{\phi }={{mr}}^{2}\dot{\phi }+\displaystyle \frac{q}{2\pi }\psi (r,z,t)={\rm{const}}.\end{eqnarray} \tag{ 69 }$$

The time derivative of equation (69) for each species σ gives

$$\begin{eqnarray}&&0=\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}({m}_{\sigma }{r}^{2}\dot{\phi })+\displaystyle \frac{{q}_{\sigma }}{2\pi }\left(\displaystyle \frac{\partial \psi }{\partial t}+\displaystyle \frac{\partial \psi }{\partial r}\displaystyle \frac{{\rm{d}}{r}}{{\rm{d}}{t}}+\displaystyle \frac{\partial \psi }{\partial z}\displaystyle \frac{{\rm{d}}{z}}{{\rm{d}}{t}}\right),\end{eqnarray} \tag{ 70 }$$

where ${\rm{d}}/{\rm{d}}{t}=\partial /\partial t+{\bf{v}}\cdot {\rm{\nabla }}$ is the time derivative experienced by an observer moving with a particle having velocity ${\bf{v}}$. Recalling that

$$\begin{eqnarray}&&{B}_{r}=-\displaystyle \frac{1}{2\pi r}\displaystyle \frac{\partial \psi }{\partial z},\ \ {B}_{z}=\displaystyle \frac{1}{2\pi r}\displaystyle \frac{\partial \psi }{\partial r},\end{eqnarray} \tag{ 71 }$$

equation (70) can be expressed as

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}({m}_{\sigma }{r}^{2}\dot{\phi }) & = & -\displaystyle \frac{{q}_{\sigma }}{2\pi }\left(\displaystyle \frac{\partial \psi }{\partial t}+\displaystyle \frac{\partial \psi }{\partial r}\displaystyle \frac{{\rm{d}}{r}}{{\rm{d}}{t}}+\displaystyle \frac{\partial \psi }{\partial z}\displaystyle \frac{{\rm{d}}{z}}{{\rm{d}}{t}}\right)\\ & = & -\displaystyle \frac{{q}_{\sigma }}{2\pi }\displaystyle \frac{\partial \psi }{\partial t}-{q}_{\sigma }r({B}_{z}{v}_{r\sigma }-{B}_{r}{v}_{z\sigma }).\end{eqnarray} \tag{ 72 }$$

Multiplying equation (72) by ${n}_{\sigma }$ and summing over species gives

$$\begin{eqnarray}&&\displaystyle \sum _{\sigma }\left(\displaystyle \frac{{\rm{d}}\ }{{\rm{d}}{t}}({n}_{\sigma }{m}_{\sigma }{{ru}}_{\sigma \phi })-{m}_{\sigma }{{ru}}_{\sigma \phi }\displaystyle \frac{{{\rm{d}}{n}}_{\sigma }\ }{{\rm{d}}{t}}\right)\\ &&\quad =\,-{{rB}}_{z}{J}_{r}+{{rB}}_{r}{J}_{z}\\ &&\quad =\,r\hat{\phi }\cdot ({{\bf{J}}}_{{\rm{pol}}}\times {{\bf{B}}}_{{\rm{pol}}}).\end{eqnarray} \tag{ 73 }$$

Summing ${n}_{\sigma }{q}_{\sigma }{v}_{r\sigma }$ and ${n}_{\sigma }{q}_{\sigma }{v}_{z\sigma }$ on the right hand side of equation (72) gives the current densities J_r and J_z. Summing ${n}_{\sigma }{q}_{\sigma }$ over species gives zero because of quasi-neutrality so no term involving $\ \partial \psi /\partial t\$ survives the summation of equation (72) over species. The continuity equation for each species can be written as

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}{n}}_{\sigma }}{{\rm{d}}{t}}+{n}_{\sigma }{\rm{\nabla }}\cdot {{\bf{u}}}_{\sigma }=0\end{eqnarray} \tag{ 74 }$$

so equation (73) becomes

$$\begin{eqnarray}&&\displaystyle \sum _{\sigma }\left(\begin{array}{c}\displaystyle \frac{\partial }{\partial t}({n}_{\sigma }{m}_{\sigma }{{ru}}_{\sigma \phi })\\ +{{\bf{u}}}_{\sigma }\cdot {\rm{\nabla }}({n}_{\sigma }{m}_{\sigma }{{ru}}_{\sigma \phi })\\ +{m}_{\sigma }{{ru}}_{\sigma \phi }{n}_{\sigma }{\rm{\nabla }}\cdot {{\bf{u}}}_{\sigma }\end{array}\right)=r\hat{\phi }\cdot ({{\bf{J}}}_{{\rm{pol}}}\times {{\bf{B}}}_{{\rm{pol}}}).\end{eqnarray} \tag{ 75 }$$

When summing over species in the first term on the left-hand side, the random velocities sum to zero so what remains is the mean velocity of the species. equation (75) then becomes

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}(\rho {{rU}}_{\phi })+{\rm{\nabla }}\cdot \displaystyle \sum _{\sigma }({m}_{\sigma }{{ru}}_{\sigma \phi }{n}_{\sigma }{{\bf{u}}}_{\sigma })=r\hat{\phi }\cdot ({{\bf{J}}}_{{\rm{pol}}}\times {{\bf{B}}}_{{\rm{pol}}}).\end{eqnarray} \tag{ 76 }$$

We may ignore electron mass since m_e is much smaller than the ion or neutral masses and note that the strong collisional coupling between ions and neutrals causes ${{\bf{u}}}_{i}\simeq {{\bf{u}}}_{n}\simeq {\bf{U}}$ so

$$\begin{eqnarray}&&\displaystyle \sum _{\sigma }({m}_{\sigma }{{ru}}_{\sigma \phi }{n}_{\sigma }{{\bf{u}}}_{\sigma })\simeq \displaystyle \sum _{i,e}({m}_{\sigma }{{rU}}_{\phi }{n}_{\sigma }{\bf{U}})=\rho {{rU}}_{\phi }{\bf{U}}.\end{eqnarray} \tag{ 77 }$$

Here any terms involving the squares of random velocities have been dropped because these would eventually give a pressure gradient in the ϕ direction, but because of the assumed azimuthal symmetry $\partial P/\partial \phi =0$.

Because ${{\bf{B}}}_{{\rm{pol}}}\,=$ ${(2\pi )}^{-1}$ ${\rm{\nabla }}\psi \times {\rm{\nabla }}\phi$, ${{\bf{J}}}_{{\rm{pol}}}\,=$ ${(2\pi )}^{-1}$ ${\rm{\nabla }}I\times {\rm{\nabla }}\phi ,$ and $r\hat{\phi }={r}^{2}{\rm{\nabla }}\phi$ the right hand side of equation (76) can be expressed as

$$\begin{eqnarray}r\hat{\phi }\cdot ({{\bf{J}}}_{{\rm{pol}}}\times {{\bf{B}}}_{{\rm{pol}}}) & = & r(\hat{\phi }\times {{\bf{J}}}_{{\rm{pol}}})\cdot {{\bf{B}}}_{{\rm{pol}}}\\ & = & {(2\pi )}^{-1}\ r(\hat{\phi }\times ({\rm{\nabla }}I\times {\rm{\nabla }}\phi ))\cdot {{\bf{B}}}_{{\rm{pol}}}\\ & = & {(2\pi )}^{-1}\ {r}^{2}\ ({\rm{\nabla }}\phi \times ({\rm{\nabla }}I\times {\rm{\nabla }}\phi ))\cdot {{\bf{B}}}_{{\rm{pol}}}\\ & = & {(2\pi )}^{-1}{\rm{\nabla }}I\cdot {{\bf{B}}}_{{\rm{pol}}}\\ & = & {(2\pi )}^{-1}{\rm{\nabla }}\cdot (I{{\bf{B}}}_{{\rm{pol}}})\\ & = & {(2\pi )}^{-2}{\rm{\nabla }}\cdot (I{\rm{\nabla }}\psi \times {\rm{\nabla }}\phi ).\end{eqnarray} \tag{ 78 }$$

Finally, using equations (77) and (78), equation (76) assumes the form

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}(\rho {{rU}}_{\phi })+{\rm{\nabla }}\cdot (\rho {{rU}}_{\phi }{\bf{U}})=\displaystyle \frac{1}{4{\pi }^{2}}{\rm{\nabla }}\cdot (I{\rm{\nabla }}\psi \times {\rm{\nabla }}\phi )\end{eqnarray} \tag{ 79 }$$

which is identical to equation (9b).

If the right hand side of equation (79) were zero, then this equation would be a conservation equation for the angular momentum density $\rho {{rU}}_{\phi }$. However, the right hand side can be finite and so act as either a source or sink of angular momentum depending on whether it is positive or negative. Because the right hand side is in the form of a divergence and because both ψ and I must vanish at infinity, integration of equation (79) over volume up to infinity will result in there being no contribution from the right hand side. Thus, if there is a source of angular momentum over some subvolume, the divergence form of the right hand side of equation (79) shows that some other subvolume must absorb all the angular momentum produced by the source. Because the flux term ${\rm{\nabla }}\cdot (\rho {{rU}}_{\phi }{\bf{U}})$ on the left-hand side of equation (79) also has the form of a divergence, integration over the entire volume up to infinity shows that all contributions from angular momentum flux will similarly cancel. The result is that the volume integral up to infinity of equation (79) gives

$$\begin{eqnarray}&&\int \rho {{rU}}_{\phi }{{\rm{d}}}^{3}r=\mathrm{const}.\end{eqnarray} \tag{ 80 }$$

so the total angular momentum is exactly conserved even though there can be local fluxes of angular momentum and local sources or sinks, i.e., local torques from the ${\bf{J}}\times {\bf{B}}$ force.

The right hand side of equation (79) describes the torque associated with magnetic forces. This magnetic torque, a source/sink for angular momentum, provides a non-mechanical means for transporting angular momentum from one subvolume to another spatially separated subvolume. This transport of angular momentum to a distant location can be understood by writing the right hand side as

$$\begin{eqnarray}&&\displaystyle \frac{1}{4{\pi }^{2}}{\rm{\nabla }}\cdot (I{\rm{\nabla }}\psi \times {\rm{\nabla }}\phi )=\displaystyle \frac{1}{4{\pi }^{2}}({\rm{\nabla }}I\times {\rm{\nabla }}\psi )\cdot {\rm{\nabla }}\phi \end{eqnarray} \tag{ 81 }$$

and examining ${\rm{\nabla }}I\times {\rm{\nabla }}\psi$. If surfaces of constant I coincide with surfaces of constant ψ then, as in the Grad–Shafranov equation, this means $I=I(\psi )$ so ${\rm{\nabla }}I={I}^{{\prime} }{\rm{\nabla }}\psi$ and ${\rm{\nabla }}I\times {\rm{\nabla }}\psi =0$ so no magnetic torque exists, or equivalently there is no source or sink of angular momentum.

In figure 2 the constant I surfaces approximately coincide with the constant ψ surfaces in the regions above and below the z = 0 plane. This corresponds to the poloidal electric current flowing along the poloidal magnetic flux surfaces in the regions above and below the z = 0 plane. Thus, above and below the z = 0 plane the poloidal flux surfaces can be considered to behave as set of concentric insulated conductors carrying the poloidal electric current.

However, if surfaces of constant I do not coincide with surfaces of constant ψ there will be a source or sink of angular momentum. This happens in the $z\simeq 0$ disk region where the current is radial whereas the magnetic field is vertical. This means ${\rm{\nabla }}\psi$ is horizontal while ${\rm{\nabla }}I$ is vertical so ${\rm{\nabla }}I\times {\rm{\nabla }}\psi \ne 0$. Thus $-{{rJ}}_{r}{B}_{z}$ is finite in the $z\simeq 0$ disk region, i.e., torques exist in the disk region. However, while J_r is always radially inwards, B_z is negative at small r and positive at large r so the torque has opposite polarities at small and large $r$. This torque reversal can also be seen from the point of view of the poloidal flux function. In figure 2 ψ is negative everywhere and has a minimum at the black circle. In the z = 0 plane and to the left of this minimum, $\partial \psi /\partial r$ is negative corresponding to negative B_z while to the right of the black circle, $\partial \psi /\partial r$ is positive corresponding to positive ${B}_{z};$ see the vertical arrows indicating B_z in figure 2.

The inward spiraling metaparticles constitute a radially inward electric current so J_r is negative in the $z\simeq 0$ disk region. Since both J_r and B_z are negative in the $z\simeq 0$ disk region to the left of the black circle, the first term in equation (73) causes $-{{rB}}_{z}{J}_{r}$ to be negative corresponding to a removal of angular momentum in this region. On the other hand B_z is positive to the right of the black circle in the $z\simeq 0$ disk region while J_r remains negative corresponding to a creation of angular momentum in the $z\simeq 0$ disk region to the right of the black circle. Angular momentum is consequently being transported from the left to the right. Unlike viscosity which transports angular momentum between adjacent regions, magnetic torque transports angular momentum to a spatially separated region. In particular, angular momentum is extracted by the magnetic torque from a small-r region and then deposited into a non-adjacent large-r region that is connected, not by geometric proximity, but instead by having the same ψ flux surface. The connectivity occurs via electric current flowing along the flux surfaces in the region external to the $z\simeq 0$ disk region and there is no need for the angular momentum ${{mr}}^{2}\dot{\phi }$ to be finite in this external region.

To see this in more detail, consider two nearby flux surfaces labeled ${\psi }_{1}$ and ${\psi }_{2}$, let V₁₂ denote the volume of the toroidal shell between these two flux surfaces, and let ${S}_{\mathrm{1,2}}$ denote the geometric surfaces associated with ${\psi }_{\mathrm{1,2}}$. The rate of change of angular momentum due to magnetic torque in this toroidal shell is prescribed by integration over the volume of the shell, i.e.,

$$\begin{eqnarray}&&{\displaystyle \int }_{{V}_{12}}\displaystyle \frac{1}{4{\pi }^{2}}{\rm{\nabla }}\cdot (I{\rm{\nabla }}\psi \times {\rm{\nabla }}\phi ){{\rm{d}}}^{3}r\\ &&\quad =\,\displaystyle \frac{1}{4{\pi }^{2}}\left({\displaystyle \int }_{{S}_{2}}(I{\rm{\nabla }}\psi \times {\rm{\nabla }}\phi )\cdot {\rm{d}}{\bf{s}}-{\displaystyle \int }_{{S}_{1}}(I{\rm{\nabla }}\psi \times {\rm{\nabla }}\phi )\cdot {\rm{d}}{\bf{s}}\right)\\ &&\quad =\,0,\end{eqnarray} \tag{ 82 }$$

where the last line is zero because ${\rm{d}}{\bf{s}}$ is parallel to ${\rm{\nabla }}\psi$ on the respective geometric surfaces ${S}_{\mathrm{1,2}}$ and because $I{\rm{\nabla }}\psi \times {\rm{\nabla }}\phi$ has no ϕ component so Gauss's law can be used without concern about the shell being doubly-connected. Equation (82) states that the angular momentum removed magnetically from the disk region between ${\psi }_{1}$ and ${\psi }_{2}$ at small radius (left of black circle) is deposited into the disk region between ${\psi }_{1}$ and ${\psi }_{2}$ at large radius (right of black circle). Thus, angular momentum is transported electrically along the flux surfaces from the inner region to the distant outer region.

This magnetic transport of angular momentum has an important consequence regarding energy. A flux surface passing through small r to the left of the minimum of ψ in the z = 0 plane (left of black circle) and passes through large r to the right of this minimum. Furthermore, the smaller the value of r is to the left for this intersection of the flux surface with the z = 0 plane, the larger the value of r is for the intersection to the right. Because the kinetic energy of angular momentum is $\ {{mv}}_{\phi }^{2}/2={({{mr}}^{2}\dot{\phi })}^{2}/(2{r}^{2})$, the kinetic energy of a certain amount of angular momentum scales as $1/{r}^{2}$. Thus, the energy required to create angular momentum at extremely large r is negligible compared to the energy absorbed on removing angular momentum at small $r$. This means that there is negligible energy cost in transporting angular momentum from small r to large $r$. While both angular momentum and azimuthal kinetic energy are removed from the same small r location, the creation of angular momentum and of energy do not occur at the same location. As will be argued below, the removed energy goes into driving the jet while the removed angular momentum is shed in the disk plane at extremely large $r$.

This transfer of angular momentum can be visualized by imagining the accretion disk to have a set of 'gear teeth' at small radius and another set of gear teeth at large radius. A DC electric generator engages the gear teeth at small radius and a DC electric motor engages the gear teeth at large radius. The rotating disk at small radius drives the generator and the load provided by the generator acts as a braking torque. At the same time, the motor at large radius spins up the plasma by an infinitesimal amount at large radius to shed the angular momentum absorbed at small radius. Because angular momentum scales as ${{mr}}^{2}\dot{\phi }$ a tiny amount of angular velocity at large radius absorbs the angular momentum associated with finite $\dot{\phi }$ at small radius. The electrical circuit connecting the generator to the motor has one wire in the z = 0 plane and another wire arching out of the disk into the region above (or below) the z = 0 plane following the poloidal flux surfaces. Thus, angular momentum is transferred from the inner region (generator) to the outer region (motor) with no mechanical contact between the two regions.

The energy budget is seen by considering the poloidal flux surfaces as equivalent to an azimuthally symmetric set of wires rising out of the z = 0 plane and carrying electric currents flowing in the poloidal direction. At small radius, these currents will attract each other and thus squeeze the plasma between the wires to produce an axial pressure gradient. This axial pressure gradient accelerates plasma in the z direction and this acceleration of matter constitutes an energy sink; off the z-axis there is additional axial acceleration provided by the magnetic force $\hat{z}\cdot {\bf{J}}\times {\bf{B}}={J}_{r}{B}_{\phi }-{J}_{\phi }{B}_{r}$. This force is dominantly given by the ${J}_{r}{B}_{\phi }$ term when the jet is nearly collimated so B_z is approximately uniform in which case ${J}_{\phi }$ is small. Because magnetic field is frozen into the plasma in the region external to the disk, the poloidal field is frozen into the jet and so is distended by the jet. This stretching of the poloidal field plus the increasing volume of toroidal field in the lengthening jet correspond to an increase of the magnetic energy in the region external to the disk. The energy extracted from the disk goes into the jet and its frozen-in magnetic field. Thus the angular momentum and gravitational potential energy of the accreting disk particles are absorbed at a single location (small radius region of the disk) but are not deposited at some other single location. The extracted energy is deposited at the jet location which is in the vicinity of the z-axis and away from the z = 0 plane; this deposited energy becomes the jet kinetic and magnetic energy. In contrast, the extracted angular momentum is deposited at near infinite radius in the vicinity of the z = 0 plane.

Maximum power coupling occurs when the load resistance equals the battery internal resistance and if this is so, the voltage at the battery terminals is half the open-circuit voltage. Bellan [1] showed that this situation was a stable operating point, i.e., any deviation would cause a heating or cooling of the disk that would change the disk internal resistance in such a way as to drive the system back to this operating point. Because the electrical resistivity in the disk results from collisions with neutrals rather than from Spitzer resistivity, the resistivity increases with temperature. Stability occurs because if the disk is at the operating point where the load resistance matches the internal resistance and there is a perturbation that increases the current, this will heat the disk and so increase the disk electrical resistance which will then reduce the current.

The open-circuit voltage occurs when the outward force of the electric field resulting from the accumulation of metaparticles at small radius balances the inward gravitational force. This is because if the electric field balances gravity, there is no more radial infall of metaparticles and so no current, i.e., an open-circuit. Thus, the open-circuit electric field is given by

$$\begin{eqnarray}&&{q}_{i}{E}_{r}^{{\rm{open}}}=\displaystyle \frac{{\alpha }^{-1}{m}_{n}{MG}}{{r}^{2}}\end{eqnarray} \tag{ 83 }$$

since a metaparticle has charge q_i and mass ${\alpha }^{-1}{m}_{n}$ where m_n is the mass of a neutral. Since the voltage under matched load is half the open-circuit voltage, the loaded electric field will be

$$\begin{eqnarray}&&{E}_{r}^{{\rm{loaded}}}=\displaystyle \frac{{\alpha }^{-1}{m}_{n}{MG}}{2{q}_{i}{r}^{2}}.\end{eqnarray} \tag{ 84 }$$

Now consider electrons and metaparticles at $r=a,$ z = 0, i.e., at the center of the black circle in figure 2 where there is a magnetic field null. These particles will experience the inward pull of gravity, outward centrifugal force, and a force from the electric field. Thus, radial force balance for the electrons at the black circle is given by

$$\begin{eqnarray}&&\displaystyle \frac{{m}_{e}{u}_{\phi e}^{2}}{a}+{q}_{e}{E}_{r}^{{\rm{loaded}}}=\displaystyle \frac{{m}_{e}{MG}}{{a}^{2}},\end{eqnarray} \tag{ 85 }$$

whereas radial force balance for metaparticles at the black circle is given by

$$\begin{eqnarray}&&\displaystyle \frac{{\alpha }^{-1}{m}_{n}{u}_{\phi {\rm{mp}}}^{2}}{a}+{q}_{i}{E}_{r}^{{\rm{loaded}}}=\displaystyle \frac{{\alpha }^{-1}{m}_{n}{MG}}{{a}^{2}}.\end{eqnarray} \tag{ 86 }$$

On substituting for E_r^loaded the electron and metaparticle toroidal velocities are prescribed by

$$\begin{eqnarray}&&{u}_{\phi e}^{2}=\displaystyle \frac{{MG}}{a}\left(1+\displaystyle \frac{{\alpha }^{-1}{m}_{n}}{2{m}_{e}}\right),\end{eqnarray} \tag{ 87 }$$

$$\begin{eqnarray}&&{u}_{\phi {\rm{mp}}}^{2}=\displaystyle \frac{{MG}}{2a}\end{eqnarray} \tag{ 88 }$$

which shows that the toroidal velocity of the metaparticles is slower than the Kepler velocity whereas the toroidal velocity of the electrons is much faster than the Kepler velocity. This is because gravity is aided by the electric force for electrons and opposed by the electric force for metaparticles. Assuming that both electrons and metaparticles are moving in the positive ϕ direction, the net toroidal current would then be negative and this would produce a poloidal magnetic field with the sense shown in figure 2. The electron motion is essentially a balance between the radially inward force from the electric field and the outward centrifugal force. The value of α at $r=a,$ z = 0 will in general not be the same as for the metaparticles having zero-canonical angular momentum and the situation might be complicated by electron-neutral collisions that tend to dissipate the toroidal current and by capture of new particles that would tend to enhance the current. However, despite these complications there would tend to be a difference between electron and metaparticle toroidal velocities and hence a current that produces the poloidal magnetic field.