Hyper-Resistive Model of Ultra High Energy Cosmic Ray Acceleration by Magnetically Collimated Jets Created by Active Galactic Nuclei
T. Kenneth Fowler, Hui Li, Richard Anantua
TL;DR
This paper presents a hyper-resistive MHD model linking dynamo processes in AGN accretion disks to the formation of magnetically collimated jets, explaining their role in accelerating ultra high energy cosmic rays beyond 10^19 eV.
Contribution
It introduces a novel hyper-resistive MHD framework that connects dynamo formation to large-scale jet collimation and cosmic ray acceleration, bridging simulations and observations.
Findings
Model predicts jet length, opening angle, and power consistent with observations.
Explains how jets accelerate particles to ultra high energies (>10^19 eV).
Provides nine observable predictions matching astrophysical data.
Abstract
This is the fourth in a series of companion papers showing that, when an efficient dynamo can be maintained by accretion disks around supermassive black holes in Active Galactic Nuclei (AGNs), it will lead to the formation of a powerful, magnetically-collimated helix that could explain both the observed jet/radiolobe structures on very large scales and ultimately the enormous power inferred from the observed ultra high energy cosmic rays (UHECRs) with energies > 10^19 eV. Many timescales are involved in this process. Our hyper-resistive magnetohydrodynamic (MHD) model provides a bridge between General Relativistic MHD simulations of dynamo formation, on the short accretion timescale, and observational evidence of magnetic collimation of large-scale jets on astrophysical timescales. Given the final magnetic structure, we apply hyper-resistive kinetic theory to show how instability causes…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 2
Figure 3
Figure 4
Figure 5
Figure 5
Figure 5
Figure 5
Figure 5
Figure 6Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Hyper-Resistive Model of Ultra High Energy Cosmic Ray Acceleration by Magnetically Collimated Jets Created by Active Galactic Nuclei
T. Kenneth Fowler
University of California, Berkeley, CA, 94720 USA
Hui Li
Los Alamos National Laboratory, Los Alamos, NM, 87545 USA
Richard Anantua
University of California, Berkeley, CA, 94720 USA
Harvard-Smithsonian Center for Astrophysics, Cambridge, MA, 02138 USA
(Received ____ ; Revised ____ ; Accepted ____ )
Abstract
This is the fourth in a series of companion papers showing that, when an efficient dynamo can be maintained by accretion disks around supermassive black holes in Active Galactic Nuclei (AGNs), it will lead to the formation of a powerful, magnetically-collimated helix that could explain both the observed jet/radiolobe structures on very large scales and ultimately the enormous power inferred from the observed ultra high energy cosmic rays (UHECRs) with energies . Many timescales are involved in this process. Our hyper-resistive magnetohydrodynamic (MHD) model provides a bridge between General Relativistic MHD simulations of dynamo formation, on the short accretion timescale, and observational evidence of magnetic collimation of large-scale jets on astrophysical timescales. Given the final magnetic structure, we apply hyper-resistive kinetic theory to show how instability causes slowly-evolving magnetically-collimated jets to become the most powerful relativistic accelerators in the Universe. The model yields nine observables in reasonable agreement with observations: the jet length, radiolobe radius and apparent opening angle as observed by synchrotron radiation; the synchrotron total power, synchrotron wavelengths and maximum electron energy (TeVs); and the maximum UHECR energy, the cosmic ray energy spectrum and the cosmic ray intensity on Earth.
cosmic ray acceleration, accretion disks, jets, hyper-resistivity
††journal: ApJ
1 Introduction
The importance of understanding the origin of ultra high energy cosmic rays (UHECRs), with energies up to , has long been recognized (Cronin, 1999). Many explanations have been offered (Bierman, 1997), none yet fully satisfactory (Blandford, Meier, & Redhead, 2019).
This paper builds on Colgate & Li (2004) in which it was hypothesized that UHECRs arise as electric currents in magnetic jets created by massive black holes inside Active Galactic Nuclei (AGNs), still a viable hypothesis (Blandford & Anantua, 2017). Colgate and Li had already concluded that only AGNs with black hole masses to could account for the total energy in UHECRs. The discovery by Balbus & Hawley (1998) that magneto-rotational (MRI) instability in rotating accretion disks around AGNs can create powerful dynamos strongly suggests that AGNs eject the magnetic jets observed by synchrotron radiation (Krolik, 1999); and General Relativistic Magnetohydrodynamic (GRMHD) simulations have produced these dynamos (McKinney, et al., 2012). Direct attempts to correlate UHECRs with known AGNs are limited by statistics but may yet pin down UHECR origins (Pierre Auger Collaboration, 2007, 2014, 2017).
It has long been appreciated that AGN dynamos could produce the volts needed to accelerate UHECRs (Lovelace, 1976; Lynden-Bell, 2006). The uncertainty has concerned what mechanism transfers dynamo voltage to ion acceleration. Unlike the transient mechanisms mentioned in Colgate & Li (2004), the mechanism proposed in this paper steadily accelerates ions ejected from the accretion disk over the full length of the jet, including the lobe regions. This acceleration mechanism is based on known plasma phenomena. Acceleration occurs in two stages. First is a weaker precursor stage due to current-driven magnetohydrodynamic (MHD) kink instability known to accelerate ions in the laboratory (Rusbridge, et al., 1997). This is followed by a powerful, purely kinetic stage based on a known ion cyclotron resonance instability that could be further explored in the laboratory, as discussed in Fowler & Li (2016).
Confidence in our accelerator model depends on confidence in the jet structure undergirding this model. We begin this paper with a brief review of our model of the jet accelerator structure which is supported by the recent new measurements from the RadioAstron mission (Giovannini, et al., 2018). This review appears in Section 2, together with a discussion of MRI-driven dynamos in Appendix A. This is followed by a discussion of ion acceleration, in Section 3 together with Appendix B; the predicted UHECR energy spectrum and cosmic ray intensity on Earth, in Section 4; synchrotron radiation as a signature of our model, in Section 5 and Appendix C; and a summary of this paper and all papers in this series, in Section 6.
Our jet model is axisymmetric, though analogous WKB solutions would apply when jets are bent by encounters with the ambient (Begelman, et al., 1984). We use a stationary system of cylindrical coordinates in which the disk spins about a fixed -axis with angular frequency pointing along the -direction in the inner region of the disk, giving positive toroidal magnetic field and negative in the same region. Except as noted, units are cgs, often introducing , the speed of light.
2 THE MAGNETIC ACCELERATOR STRUCTURE
We begin with a review of our disk/jet model in Colgate, et al. (2014)—hereafter Paper I— and jet structure and stability in Colgate, et al. (2015)—hereafter Paper II.
2.1 Jets as Current Loops
Together, Papers I and II provide the basis for the final quasi-static magnetic structure yielding steady-state acceleration of cosmic rays.
This final structure is shown in Figure 1, taken from Paper II, showing an - cross-section of vertically-expanding magnetic flux surfaces that is one way to describe the jet produced by the spinning accretion disk at the bottom of the figure, and often the description provided by MHD simulations of disk/jet systems. For our purposes it is more useful to note that the structure in Figure 1 is mainly produced by electric current flowing vertically up the Central Column, then bending to form the jet “nose” where the majority of cosmic ray acceleration will take place. Following Lynden-Bell (2003), we will interchangeably refer to jets as the Central Column or as a “magnetic tower.” We postpone to Section 2.8 to discuss the role of the weaker but important current surrounding the Central Column, labeled Diffuse Pinch in Figure 1.
As is shown in Figure 2, the Central Column current finally loops through the disk, serving as an electric circuit extracting power from the disk. In our model, the Central Column current is ultimately self-collimated by the magnetic pinch force, , creating a magnetic tower (or “pinch”). As noted in the Abstract, we have shown in Fowler & Li (2016) – hereafter Paper III – that the existence of a self-collimated jet current automatically produces the electric fields that accelerate cosmic rays, as is further elaborated here in Section 3.
But besides , a spinning accretion disk produces an electric field that actually turns out to be the force that launches jets from an accretion disk (McKinney, et al., 2012). Thus how jets launched by evolve to magnetically-collimated “towers” becomes a vital part of the story.
2.2 How Magnetically-Collimated Jets Are Launched by an Accretion Disk
Magnetically-collimated jets produced in the laboratory employ plasma guns driven by a capacitor bank (Hooper, et al., 2012; Zhai et al., 2014; Fowler & Li, 2016). Accretion disk jets are driven by disk rotation , creating a dynamo similar to a Faraday disk but with self-excited growing magnetic flux (Balbus & Hawley, 1998), given any pre-existing seed field near an evolving black hole (Appendix A.3). We find that jets created by rotation can be fundamentally different from laboratory jets, yet they evolve to the same thing, but in a manner not easily accessible to GRMHD computer simulations of accretion disks.
The main features of jet evolution occur in three steps:
- (1)
First the dynamo electric field launches a jet at velocity when the dynamo voltage drives the dynamo current up to a value given by with free-space impedance in our units (Blandford & Anantua, 2017). This jet is not collimated but instead has the conical shape sketched in Figure 3 (Tchekhovskoy, et al., 2008). 2. (2)
Next, as the dynamo current continues to rise, the jet current remains fixed as long as the jet velocity persists at . The growth in dynamo current greater than the jet current forces a portion of the dynamo current to short-circuit through the disk corona as shown at the bottom of Figure 3. The short-circuit is produced by acceleration of the jet to speed . 3. (3)
Finally, anything that slows down eliminates the short-circuit, allowing the conical jet to evolve to a collimated magnetic tower with current equal to the full current produced by the dynamo. Why this last stage of jet evolution is not observed in GRMHD simulations is discussed in Section 2.7.
Figure 3 again calls attention to the utility of describing jets as current rather than flux. The current in Step (1) above is the total current produced by the accretion disk acting as a self-excited dynamo creating its own magnetic field (Balbus & Hawley, 1998). The current is the maximum current in the system, referred to interchangeably as the disk current or the dynamo current.
Loops in Figure 3 represent the main current paths emerging near the axis of symmetry and returning to the disk. Figure 3 shows a branching of the dynamo current as it emerges to form a jet. The initial branch forms the jet labeled Cone; followed by the short-circuit as the dynamo current grows while the jet current stays constant; and finally the Tower branch grows as the Cone dies so that the tower jet becomes the total current produced by the dynamo.
2.3 Dynamics of Accretion Disk/Jet Structures: Mathematical Model
The conclusions in Section 2.2 can be derived from the following relativistic current and momentum equations, with gravitational potential to represent the black hole. General relativity included in GRMHD codes is not required except very near the black hole, leaving then special-relativistic fluid equations in three dimensions (3D) derived from Vlasov distribution functions with position , momentum and time in our reference frame. We obtain:
[TABLE]
Here is the current density where sums over ions (hydrogen) and electrons with charge , rest mass ; is the bulk flow momentum; is the charge density; K.E. represents small kinetic terms; and is the particle velocity giving for example a relativistic correction arising from in an integration by parts (Montgomery & Tidman, 1964).
Gravity in Equation (1b) divides the system into a denser accretion disk and lower density jet ejected in the disk corona. We obtain mean-field (2D) equations for the jet by averaging over and fluctuation scales in time and space. Ambient pressure represents ionized gas outside the accretion environment. Dropping the ambient pressure for the moment, we obtain:
[TABLE]
Equation (2a) is Ohm’s Law obtained by dropping terms in Equation (1a) as discussed in Section 3, with fluid velocity , vector potential , electrostatic potential and hyper-resistivity representing any kind of turbulence, MHD or kinetic (Fowler & Gatto, 2007). Equation (2b) is the momentum equation, with , and it includes the electric force due to space charge density , given in Equation (2c). Equation (2d) is the solution to Equation (2b) giving as acceleration-driven current flow between flux surfaces moving at velocity in our reference frame. Equation (2e) is the Force Free Degenerate Electrodynamics (FFDE) approximation to Equation (2b) (Meier, 2012), assuming that hyper-resistivity .
2.4 Step 1: The Formation of Conical Jets
We now apply Equations (2a) - (2e) to see how theory corresponds to the three steps of jet evolution described in Section 2.2, from Cone to Short Circuit to Magnetic Tower. Three distinct regions of the disk/jet system will emerge associated with three different timescales within the overall lifetime of the jet ( yrs):
- (1)
The disk itself: where the dynamo current reaches quasi-steady state in a few accretion times (only years near the black hole); 2. (2)
The disk corona: home of a persisting short-circuit that allows continuation of a conical jet for a duration of yrs (see Section 2.7); 3. (3)
The jet as magnetic tower: enduring for the remaining .
It is the great disparity in timescales that makes it difficult to span the entire evolution of a magnetic tower in GRMHD disk+jet simulations which focused on the early phases of dynamo formation; or in special relativistic MHD simulations which focused on jet propagation.
2.4.1 Initial Conditions
As in many GRMHD simulations, we consider an initial state in which a pre-existing black hole is embedded in an accretion disk with rotation frequency , threaded by a vacuum poloidal seed field created by distant currents in the ambient. Both dipole and quadrupole seed fields eject jets upward and downward, the upward jets being similar except very near the disk midplane (see Paper I, Appendix B). Figures 1 - 3 depict an upward jet.
In the absence of MRI, the ideal Ohm’s Law only allows gravity to compress flux, sometimes observed in simulations. Here we are interested in the formation of a self-excited dynamo that grows its own magnetic field to the kilogauss levels inferred from AGN observations (Paper I, Balbus & Hawley 1998). Rotation immediately causes an accretion disk to charge up to (positive for our sign convention) with a charge density . Because the ordinary resistive conductivity is negligible in disks, no dynamo current flows through the disk until the MRI instability sets in, creating hyper-resistive current flow radially across flux surfaces inside the disk (as ordinary resistivity allows inside a Faraday disk). Closure of this MRI-driven current in or above the corona allows a jet to be ejected.
2.4.2 Two Quasi-Steady States
We note that, for the quasi-steady state described by Equation (2e), there are two jet solutions, one with and ; and another with giving . Since initially both and are zero in the vacuum seed field, before a jet is launched is zero, giving, by Equation (2b), jet acceleration that saturates at . In Equation (2e), this corresponds to , possible if with solution as the threshold to launch a jet. This is the conical jet solution in GRMHD simulations (McKinney, et al., 2012), conical rather than collimated because cancels the pinch force. A second solution exists, however, when . In this limit, a jet can be launched by which gives jets as magnetic towers in the laboratory (Hooper et al. 2012; Paper II). Finally, we call attention to the kinetic term in that is dropped in the FFDE approximation but plays a vital role in sustaining a short circuit in the disk corona, as discussed in Section 2.5.
2.4.3 Which Comes First - Cone or Tower?
As the dynamo current builds up, the conical jet comes first because more dynamo current is required to launch a tower jet than a cone. This is shown in Figure 4, which plots the cumulative jet current in units of , where with (Schwarzschild radius). The shape of is only approximate at . But for , the curves are calculated from the radial force equation with angular momentum as a constraint, given by with Keplarian rotation, in Equation (6c) below. The lower curve for the conical jet gives with from force balance in Equation (2e) with , giving . The upper curve for the tower is taken from Figure 2 in Paper I, obtained by an exact solution of with the angular momentum constraint. The total current asymptotes to a constant value. Note that for the cone peaks at , indicating a reversal in sign of the jet current density noted in simulations, while, for the tower, is merely the radius of the Central Column. The value of and the ratio of cone and tower currents at are derived in Section 2.8.
2.4.4 Continued Growth of the Disk Current
Next we consider the timescale of processes inside the disk. The growth of the disk current in terms of the vector potential is described by Equation (2a) with equal to that for MRI. Near the black hole, the growth rate will turn out to be of the order of the accretion rate , giving a rapid buildup of dynamo current to balance gravitational energy input by accretion (only years, accessible by GRMHD, as noted above). We obtain inside the disk:
[TABLE]
We note that, as the disk develops MRI, these equations show that gradually grows from zero, until reaching which launches the conical jet, but then continuing smoothly to the value needed to sustain a magnetic tower. It is this gradual growth of inside the disk, independent of the jet, that led us to postulate the creation of a short-circuit in the corona isolating the growing current in the disk from a constant current in the jet. Though initially , it grows smoothly to equal and surpass , as discussed below.
In Equations (3a) and (3b), we apply Equation (3c) to eliminate zeroth-order terms, leaving as the driver of the self-excited dynamo magnetic field. Then Equation (3d) gives the steady state with obtained as follows. We set and drop compared to inside the disk, giving then and in quasi-steady state. And most importantly, we apply Equation (3e), derived from an ordering scheme in Appendix A of Paper I, and reaffirmed for an MRI-driven in Appendix A of this paper (giving in Paper I). Equation (3d) shows why in the disk grows to exceed in the disk, enough so to drive the jet current up to the magnetic tower threshold in Figure 4.
2.5 Step 2: The Short Circuit
Why a magnetic tower does not form immediately concerns the persistence of the initial conical jet. This is due to the jet structure whereby the jet current returns close to the outgoing jet so that jet magnetic energy only changes as the length changes, at a rate (Paper II, Appendix C; Blandford & Anantua (2017)). As was noted in Step (2), Section 2.2, a short-circuit is required to sustain a constant jet current as the dynamo current rises. The required short-circuit is automatically accommodated by in Equation (2d). This occurs through the kinetic term of that allows a quasi-steady conical jet to persist even as the disk current continues to grow, giving then the short circuit current and a corresponding parallel current as follows:
[TABLE]
Here is the width of the current path and the factor is the final disk current with factor fitted to Figure 4 at . Equation (4a) is obtained from Equation (2d) with from Equation (2b) with . The cross-product of with gives ; and maximizing on the divergence gives . The final expression approximates the integral over by the fluid velocity divided by the Alfvén-like velocity in Equation (4d). Equation (4c) is Ohm’s Law in the conical jet, whereby jet acceleration to maintains in the jet.
2.5.1 Parallel Current, Circuit Notation
The distinction between and in Equations (4a) and (4b) reflects the branching of the total jet current in Figure 3. While the MRI-driven current loop for a magnetic tower must close along closed poloidal field lines, the current loop of the initial conical jet closes via a current flowing directly across poloidal flux surfaces. This flows across the nose of the growing jet. But because the conical jet current remains constant as the disk current continues to increase, a short-circuiting must also flow in the corona, giving the disk current as the sum of three branches in Equation (4e). We take the conical jet current equal to its return current at the nose so that Equation (4e) gives in Equation (4b).
In the next Section, we explicitly describe Figure 3 as an electric circuit. Here and hereafter, the circuit variables and are defined as follows. The symbol will be the disk or dynamo current, the largest current in the system, ultimately equal to the jet current in a self-collimated magnetic tower. As in Figure 4, , the current inside the Central Column . The symbol will denote the dynamo voltage, defined as the voltage across the Central Column where it intercepts the disk, as in Step (1), Section 2.2, giving .
2.6 Step 3: Slowing Down of the Jet, Cones to Towers
As in Paper II, we describe jet propagation by an electric circuit satisfying energy conservation. We obtain:
[TABLE]
[TABLE]
Equation (5a) is obtained in the usual way by adding the results from dotting into Equation (2b) and dotting into Maxwell’s equation, using Equation (2a) giving if we neglect in the jet. Integrating Equation (5a) over the jet volume gives Equation (5b), with dynamo power on the left hand side (from the Poynting vector in Equation (5a)), and in Equation (5e) a dissipation factor with as the efficiency of converting magnetic energy to ion acceleration and a factor representing ambient shocks, as in Paper II. As in Figure 3, Equation (5b) approximates the jet as a conical jet with current embedded inside a magnetic tower with current . We drop and divide into three parts: (a) the Central Column giving in the bracket ; (b) the conical jet giving the term ; and (c) the enveloping magnetic tower gives the term , valid when , causing current to twist around field lines. As in Paper II, we use the fact that, as falls below , Alfvén waves at speed (for the jet density in Section 2.8) spread flux radially to produce a blunt nose like Figure 1, whereby the jet can be approximated as a cylinder of fixed radius and expanding length approximately independent of so that can be removed from the integral to give Equation (5c).
Equation (5c) derives the jet velocity from Equation (5b) in quasi-steady state (constant ), using Equation (4a) for and Equation (4b) for . The slow growth of the jet current to match the final dynamo current is approximated by giving the correct limits both as the jet slows down and when the jet is first launched with jet current .
Formulating jet propagation as an electric circuit has the advantages that the black hole region omitted in Equation (1b) is included in the circuit (McDonald & Thorne, 1982; Frank, et al., 2002), and the magnetic tower can be solved analytically at , giving yielding the logarithmic factor in Equation (5b) (Papers I & II). That also holds inside the conical jet, giving for a cone with outer and inner radii and , can be derived from the -component of Ohm’s Law in Equation (2a), with in the jet. We obtain:
[TABLE]
The right hand side of Equation (6a) is the square of Equation (4c), with Keplerian rotation for black hole mass and Newtonian gravitational constant . Solving Equation (6a) for gives exactly the middle expression in Equation (6b) using angular momentum conservation discussed in Paper II, giving Equation (6c) with fields defined at the disk corona () and accretion rate . We approximate the relativistic Lorentz factor by with fluid velocities and . Then solving for and expanding gives Equation (6d), valid if , true at all with for numbers below. Substituting Equation (6d) into Equation (6b) gives for the conical jet.
That viscosity can be neglected in Equation (6c) is key to our model, whereby disk/jet angular momentum conservation governs all model predictions of jet properties and UHE cosmic rays and synchrotron radiation in Table 1 below. In Appendix A of Paper I, we derive viscosity as giving for our consistent ordering and equipartition of turbulence confirmed for MRI in Equation (A2.b) as noted earlier. (These conclusions are true at all despite an error in Paper I in estimating as a projection perpendicular to that underestimates at large .)
2.6.1 How Cones Become Towers
As anticipated in Appendix C of Paper II, it is the slowing down of the jet by Equation (5c) that yields the magnetic tower jet structure in Figure 1 giving the accelerator in Section 3.
In Paper II, this was shown by following the dynamical history of , given here by Equation (4c) yielding . Then anything that slows down the jet—dissipation and/or induction—implies in Equation (2e). Extended to two dimensions, Equation (2e) with becomes the Grad-Shafranov magnetic tower solution in Figure 1, as derived in Paper II.
2.6.2 Eliminating the Short-Circuit
The new result in this paper is the explicit origin of the short circuit current by jet acceleration, in Equation (2d), as is required for the conical jet current to close on itself, in Figure 3, while the parallel current producing a magnetic tower must follow poloidal flux surfaces.
That slowing down dissipates can be shown by re-writing Equation (2b), approximating with from Ohm’s Law. We obtain:
[TABLE]
where we now retain the effect of ambient back pressure in in Equation (5e). Equation (7b) is obtained by factoring out with poloidal flux , with for jet densities in Section 2.8. Then slowing down eliminates for typical ambient pressures over most of the nose surface, even though the integral effect of applying a weak back pressure over the stiff, blunt nose in Figure 1 does slow down by 50% in Equation (5e), as discussed in Paper II.
The purely inductive rate of slowing down a conical jet when ambient pressure and dissipation are negligible is difficult to estimate. If the conical jet were a straight cylinder in a constant field, would never change; would remain zero; and the conical jet solution of Equation (2e) giving would persist. Given , both the radial pinch force creating a magnetic tower and the hoop force causing growth of the radiolobe radius are zero. The actual growth of inductance is associated with the 2D structure of the conical jet in Figure 3. We have not been able to calculate this 2D structure analytically, and there is reason to believe its radial profile cannot be in steady state. For example, the base of the jet should grow at the accretion rate, while GRMHD solutions in length run for only a few accretion times.
2.7 Evidence of Magnetic Collimation: Simulations and Observations
As already noted, attempts to simulate the entire astrophysical jet formation cycle are limited by the extreme range of time and space scales, from accretion disk dynamos concentrated near the black hole, where conical jets are created, to the Mpc dimensions of fully developed jets. In particular, GRMHD codes yielding MRI-driven dynamos producing conical jets in a low density ambient do not yet exhibit the slowing down required to produce a tower (Tchekhovskoy, 2017). That this is perhaps understandable given the limited timescale of these simulations follows from Equation (5c) dropping ln(R/a) to obtain, for zero dissipation for typical cone dimensions in GRMHD simulations. On the other hand, assuming continued inductive slowing down of the jet for reasons discussed at the end of Section 2.6.2, as dissipation develops eventually giving in later sections, yielding as in Equation (5d). Equation (5d) predicts early slow-down to at 1 of the final jet length if we approximate in Equation (5d), indicating magnetic collimation over most of the life of the jet.
The slowing down by ln that produces a magnetically collimated tower is consistent with MHD simulations of magnetically-dominated jets, relativistic in Guan, et al. (2014), non-relativistic in Nakamura, et al. (2006); Carey (2009); and Carey, et al. (2011). Disk-like creation of jets by rotation of a conducting sphere in a dense ambient is featured in relativistic MHD simulations in Bromberg & Tchekhovskoy (2016), exhibiting jet slowing down by the ambient and the two classes of MHD kink modes used in our accelerator model (Nakamura, et al., 2006, 2007; Tchekhovskoy & Bromberg, 2016).
Future work might disclose other connections to our model, for example, the short-circuit current in the disk/corona in GRMHD simulations. Preliminary evidence of a short-circuit is shown in Figure 5 analyzing data from a GRMHD simulation in McKinney, et al. (2012). Shown are trajectories of (to add contrast) in Cartesian x-z planes cutting through 3D current structures ( being at ). The predicted short-circuiting current loops are indicated by overlays in Panels (c) and (d), occurring within a radius , this being the zone of greatest gravity-driven MRI activity as calculated in Section 2.8. The tilting geometry of actual loops gives evidence of looping both in in Panel (c) and in in Panel (d).
The creation and propagation of magnetic towers (or pinches) is well understood from laboratory experiments on spheromaks (Hooper, et al., 2012; Zhai et al., 2014). Spheromaks using plasma guns to replace the accretion disk produce sub-Alfvenic, collimated, blunt-nosed structures like Figure 1 that justify our calculation of independent of in Equation (5c) (see Paper II). This collimation persists the full length of the jet no matter how long the jet length grows. That jets produced by AGNs are also magnetically collimated over long distances is now supported by new radio telescopic observations (Giovannini, et al., 2018) which corroborate earlier evidence from Faraday rotation measurements indicating the presence of collimated current channels with radii of order 10 (Schwarzschild radii) far from black holes (Owen, et al., 1989; Kronberg, et al., 2011; Lovelace & Kronberg, 2013); and the corresponding evidence near the black hole (Zamaninasab, et al., 2014; Kim, et al., 2018). In Section 5 and Paper II, we show that contrary evidence based on conical synchrotron images represented by the dashed line in Figure 1 can be explained as a few wandering field lines that spread the radiated power but do not spoil the overall collimation of the current. Field lines wandering around magnetically-collimated jets have been observed in three-dimensional (3D) magnetic tower simulations, both non-relativistic (Nakamura, et al., 2006) and relativistic (Guan, et al., 2014). The large dimensions of the nose agree with data on radio lobe dimensions, represented by the closed flux in Figure 1 (Diehl, et al., 2008). And evidence that this huge jet/lobe boundary pushes away the ambient as predicted is shown by bubbles seen in galaxy clusters (McNamara & Nulsen, 2007).
2.8 Jet Parameters
We conclude this Section with a review of numbers needed for our accelerator model in later sections. Following Paper III, we note that the key quantities needed concern the final jet velocity and numbers associated with the Central Column where most of the gravitational power is deposited. As in Papers I and II, the Central Column radius is defined as the inner radius of a Diffuse Pinch characterized by Keplerian rotation in the disk, while the Central Column at is described by the electric circuit in Equation (5b). The Diffuse Pinch helps guide and stabilize the Central Column flow, as discussed in Paper II. We obtain:
[TABLE]
where Equation (8d) is the angular momentum conservation equation, Equation (6c), evaluated at . Here the subscript denotes quantities at , and is the efficiency of converting accretion power into Central Column jet power , with current and voltage with fitting factor . The accretion rate is taken to be with typical jet lifetime years (Colgate & Li, 2004); (Beskin, 2010). Simultaneous solution of Equations (8a) and (8d) gives:
[TABLE]
Here is in units of solar masses; Keplerian and ; (half from shocks, half from ion acceleration, as in Paper II); and we take from Paper I and as discussed in Paper II, Appendix B. Equation (9d) is obtained from Equation (9c) using for magnetic towers, from Paper I; to be compared with that when a conical jet is launched in Equation (9e). The ratio accounts for current ratio at in Figure 4, using in Equation (8b). Equation (9g) applies Equations (9a) and (9b) to Equations (8b) and (8c). Note that and the total jet power IV are unchanged (except as changes) during the transition from a conical jet to a magnetic tower. Note also the Diffuse Pinch current omitted here [() at in Figure 4], though of the total, only accounts for about of the total dynamo power.
The other number we need for accelerator calculations is the density in the Central Column, given in Equation (9f), where and . As is discussed in some detail in Paper II, Appendix A, the fact that the magnetic tower jet is ejected vertically with negligible requires that an electrostatic sheath form to eject ions against gravity (for our sign convention).
Why the jet current becomes the constant value in Equation (8b) is due to macrostability of the system, as discussed in Paper II, based on earlier results in Fowler, et al. (2009). That Equation (9f) is always correct, despite mass loading along the jet (not expected, Paper II) or pair-production, follows from current suppression by two-stream instability. How two-stream scattering of electrons (and positrons) kills their current in our reference frame is shown in Equation (C5). That the net ion density cannot exceed that in Equation (9f) follows from acceleration of all ions to speed c, requiring that radial hyper-resistive transport expel any excess of ions beyond emerging from an electrostatic sheath in the disk corona.
Another feature of magnetic towers is their radial confinement by ambient pressure that finally limits the growth of even as continues to grow. Following Lynden-Bell (2003), the final and other jet parameters are given by:
[TABLE]
where again years is the typical AGN jet lifetime; and in is the radius of giant radiolobes bounded by the ambient pressure (Diehl, et al., 2008). Earlier in time, the ambient pressure in Equation (10a) is replaced by ram pressure due to jet inertia. Following Paper II, we take for typical radiolobe dimensions. The final jet length in Equation (10c) is consistent with observations (Krolik, 1999).
2.8.1 Importance of Large R and Current Loop Closure
Our accelerator model will require that current loops begin to close at the nose, and that the nose radius be very large as predicted in Equation (10d). That current loops will try to close follows when slowing down gives , equivalent to Equation (7b) when . Then the procedure giving Equation (4a), applied to the nose using Equation (10a) to eliminate , gives , showing that for large R there is essentially no short-circuit over most of the nose of a slowly-evolving jet. Absent a short-circuit at the nose, in the jet return so that, to access dynamo power inside the disk, the return current must follow flux surfaces encircling an O-point inside the disk, located at where . A quasi-static poloidal seed field necessarily has such an O-point which, however, may be inaccessible at speed c. Even so, Equations (3a) and (3b) giving simultaneous growth of exceeding the seed field in the disk would create a peak in near giving an O-point where , growing eventually to as derived in Paper II.
3 ION ACCELERATION IN AGN JETS
We now begin our discussion of how magnetic energy is converted to UHE cosmic rays, using known physics of plasma turbulence discussed in Appendix B. For simplicity, we assume ions to be protons, known to be constituents of the most energetic UHECRs (Cronin, 1999), though the model applies to any ion species. As noted in Section 5, we neglect positron production. The model also describes electron acceleration yielding synchrotron radiation, discussed in Section 5 and Appendix C.
We begin with two points distinguishing our model of UHE cosmic rays. First, while AGN jets are mainly observed by synchrotron radiation, the synchrotron power is known to be a negligible fraction of overall AGN luminosity and electrons play a secondary role in our model of ion acceleration. Second, as discussed in our Paper III, our cosmic ray accelerator model is a two-stage ion accelerator. The first stage occurs in the Central Column, limited by ion radiation to energies well below UHECRs. The second stage occurs in the nose-end of the jet. In both stages, acceleration is due to plasma turbulence (hyper-resistivity in Ohm’s Law). What distinguishes these two regions is the size of the ion Larmor radius . In the Central Column, giving MHD current-driven kink modes as the main source of hyper-resistivity. In the nose with magnetic flux width , the fall of finally yields whereupon ions resonating with electron drift waves are known to produce a powerful turbulence (specifically the DCLC instability explained in Section 3.3.2) that could both accelerate ions to eV energies and eject enough ions as cosmic rays to account for the UHECR intensity on Earth (see Section 4).
The electrons behave passively in DCLC turbulence so they are not accelerated in the nose, while ions resonant with the electron waves in the nose are strongly accelerated. For a different reason, electrons play a passive role in kink instability in the Central Column. Again this is due to kinetic effects not included in MHD; namely, the two-stream instability between counter-flowing ions and electrons accelerated by the kink mode turbulence. Because of their difference in rest mass, two-stream instability scatters electrons sufficiently to eliminate the electron current in our reference frame, while ions are not much affected. Thus the kink mode accelerator in the Central Column produces current as a mono-energetic ion beam with sufficient energy to produce DCLC instability in the nose.
The mathematics supporting the above picture is developed below. A key feature is the role of turbulence-driven hyper-resistivity in Ohm’s Law. Like earlier cosmic ray models, hyper-resistivity depends on correlations in turbulent fluctuations to give additive additions to particle energy. Self-correlated hyper-resistive ion acceleration by kink modes was observed in the SPHEX experiment, discussed below (see Paper III). Two-stream instability is discussed in Appendix C.2.
3.1 Acceleration by Turbulence
Our accelerator model begins with the relativistic acceleration equation, with acceleration of momentum of a particle given by:
[TABLE]
where represents radiation loss. While taking confirms that only can accelerate ions parallel to in our reference frame fixed in the disk, by now a number of magnetic acceleration mechanisms have been identified in which invokes in moving structures (shocks, clouds, etc.). Acceleration by turbulence can be thought of this way, whereby accelerated particles encounter fluctuating electric fields in a self-organized way in a plasma.
Formally, in our model energy conversion comes from where the integral is over the entire structure of the magnetic jet and represents an average over the toroidal angle and an average over fluctuations. The integral is dominated by the Central Column which carries most of the power to the nose in Figure 1. The average yields a turbulence-generated that could serve as a quasi-steady accelerator. Most of the kinetic power is contained in , accounting for the observation of quasi-steady acceleration of ions in the SPHEX spheromak mentioned above. In this Section we will focus on the quasi-steady acceleration of ions, and defer electron acceleration to Section 5. To calculate ion acceleration, we determine from the relativistic form of Ohm’s Law derived from Equation (1a), now retaining the relativistic correction and the Hall term. We drop terms on the left hand side to obtain:
[TABLE]
where . The Hall term applies only to giving MHD. From the radiation discussion below, we will learn that always the ion Lorentz factor , which allows us to order all ion terms in Ohm’s Law as the ratio of electron rest mass to ion rest mass, hence negligible. Dropping the Hall term, we obtain for relativistic electrons:
[TABLE]
where we use and by two-stream instability discussed in Section C.2. Here, we write out explicitly giving hyper-resistivity for MHD fluctuations. The bracket indicates an axisymmetric smooth average over fluctuations giving Equation (3a). Equation (13a) yields MHD jet propagation in Section 2 while is the relativistic accelerator in the Central Column of the jet. Note the crucial role of two-stream instability in defeating a relativistic cancellation of the accelerator parallel to in our reference frame.
3.2 Ion Acceleration in the Central Column
We make the assumption, justified later, that ions and electrons in the Central Column can be described by orbits consisting of circles with small Larmor radii , gyrating around a “guiding center” momentum mainly directed along magnetic field lines. Using Ohm’s Law, the acceleration equation, Equation (11), becomes:
[TABLE]
where , and are 3D fluctuations around the mean fields due to kink modes driven by the jet current, and for MHD we drop the term in Equation (13c), leaving with . Equation (14b) can be rewritten as:
[TABLE]
The main issue is whether actual 3D fluctuations correlate to produce a finite , especially for ideal kink modes argued to dominate behavior in the Central Column of astrophysical jets in Paper II. That Equation (15) with ideal MHD perturbations does produce a finite mean electric field was demonstrated by careful measurement during kink-mode instability in the SPHEX experiment under conditions when resistivity was negligible (al-Karkhy, et al., 1993). Acceleration of ions by this electric field has been verified directly (Rusbridge, et al., 1997).
Theoretical evidence that ideal MHD fluctuations can correlate is shown in Figure 6 giving just the inductive contribution to produced by ideal MHD kink modes (McClenaghan, et al., 2014), using the non-linear non-relativistic GTC PIC code (Deng, et al., 2012).
3.2.1 Hyper-Resistivity as Diffusion
Given favorable correlations, we estimate the magnitude of as follows. We write as:
[TABLE]
using and by Equation (9g). Equation (16) represents magnetic relaxation by hyper-resistive current diffusion (Fowler & Gatto, 2007), with in Equation (13b) giving symmetrically averaged over fluctuations and . Note that has units of resistivity while (which we are calling hyper-resistivity) has units of the electric field. The term with relates hyper-resistivity to diffusion of the current, by analogy with classical diffusion via ordinary collisional resistivity. This characteristic relationship of hyper-resistivity in any direction producing diffusion perpendicular to will recur in considering ion acceleration by diffusion in Equation (21b).
We expect fluctuation amplitudes determining to saturate at levels not exceeding those required to flatten the known stability parameter over a radius where internal kink modes are localized, giving then the local source of free energy driving fluctuations with little effect on the current profile at (Fowler, 1968). To extrapolate to astrophysical dimensions, while avoiding the need to calculate correlated fluctuations in detail, we approximate , giving , where is the time required to reach quasi-steady state. Taking gives the result for the accelerating electric field on the far right hand side of Equation (16).
3.2.2 Voltage Drop in the Central Column, Curvature Radiation
We apply Equation (16) to averaged over the duration of the jet, giving, at any elapsed time , from Equation (10b). Using also Equation (16), we obtain:
[TABLE]
That this is independent of the elapsed time is a result of localization of internal kink modes to mentioned above. Then, even if the disk current profile itself does not relax, within a few wavelengths the local effect of internal kink modes can cause the jet current profiles to relax toward a stable state. As the system evolves non-linearly, instability growth rates and diminish as the current density approaches this stable profile. Then, the longer the duration, the weaker is the time-averaged needed to spread the current toward a stable profile. Implicit in this argument is the slow evolution of the profile driving kink instability, on the timescale , which is also the lifetime of the jet.
The in Equation (17) represents dissipation in the Central Column that depletes the power available to accelerate ions to UHECR energies, but only slightly because of the large inductance giving the slow evolution of the jet length, . The dissipation is in the form of acceleration of electrons giving the correct power for observed synchrotron radiation, in Section 5, and ion acceleration that is also mostly dissipated as radiation due to charged particles following the curvature of twisted magnetic field lines. By Equation (14a), the maximum energy allowed by curvature radiation is given by with (Jackson, 1998, Equation (14.31)).
[TABLE]
Here the right hand side gives the radiation with for relativistic ions or electrons following field lines with curvature radius for which the Lienard factor . The maximum Lorentz factor allowed by is independent of particle mass, given by:
[TABLE]
To get , we have used Equation (18) with by Equation (17), line curvature , and all the values from Section 2.8. A few ions at have . Note that representing the maximum acceleration energy in the Central Column is far below () required to explain UHECRs, hence the need for additional acceleration in the nose. On the other hand, that the Central Column is otherwise an excellent accelerator follows from a calculation of the ion Larmor radius and pressure parameter in the Central Column, giving:
[TABLE]
where in Equation (20c) is the ion velocity. In Equation (20a) the numerator is the ion energy and we use the density in Equation (9f) with , and in Equation (8b), and in Equation (20b), we use with at and in Equation (19b). We see that the low beta giving strong magnetic collimation is coincident with small giving good ion confinement.
That ions in the Central Column should run away to the energy in Equation (19b) and that these runaway ions carry most of the current is shown as follows. Briefly, for from Equation (17) with numbers from Section 2.8, both ions and electrons reach a maximum at a distance . But two-stream instability between the oppositely-accelerated relativistic ions and electrons spreads electron momentum between . As shown in Appendix C.2, this kills the electron current, while by momentum conservation the corresponding ion momentum spread is less by a factor so that the ion distribution is approximately a delta-function around carrying current at speed . Details of two-stream instability are discussed in Appendix C.
3.3 Ion Acceleration in the Nose
Internal kink instability due to occurring in the Central Column, as analyzed in Paper II and used in Section 3.2, should persist for some distance into the nose driven by as the current turns radially in the nose, now less affected by curvature radiation as falls like giving curvature parameter . But a new acceleration mechanism arises due to conditions in the nose region different from those in the cylindrical Central Column. One new condition concerns the sensitivity of the curvature radiation, Equation (19a), to that allows ions to be accelerated to the full voltage in the nose region where becomes large. Another new condition concerns the escape of ions from the nose, escape itself turning out to be a necessary ingredient for the powerful acceleration mechanism not present in the Central Column.
There are two ways that ions can escape from the nose. First, while ion orbits are magnetically confined in the Central Column, in the absence of an electric field ion orbits in the nose are not confined due to drifts of guiding centers in due to curvature of field lines and gradients in the magnetic field intensity. Opposite drifts by electrons establish an electric charge giving , the electron temperature, which can cause ions to escape by motion. In addition, ions can also escape by hyper-resistive diffusion.
3.3.1 Hyper-Resistive Diffusion in the Nose
As in Section 3.2.1, we represent hyper-resistivity in the nose as diffusion, noting that an electric field larger than that produced by charge separation is created if turbulence in the nose causes ions to escape by diffusion. To see this, consider the radial component of the relativistic Ohm’s Law in Equations (13a) and (13b). In the nose, with unit vector , and electron flow is either that due to drifts in the -direction or that parallel to giving so that we can drop the relativistic correction to . Applying Equation (13a) in the -direction (perpendicular to ) gives:
[TABLE]
Equation (21b) giving the cosmic ray accelerator is obtained as follows. Early in jet propagation when in the jet and nose, all of the voltage drop in the jet (giving net zero voltage around the loop) can be approximated as which appears as an inductive electric field in our disk-centered reference frame. It is this inductive field that serves to build up behind the nose advancing at velocity . Even as shocks and cosmic ray acceleration slow down by Equation (5b), it continues to be true that it is that builds up , giving Equation (21b) with in our reference frame serving as the cosmic ray accelerator in the nose. The final step relates the radial hyper-resistivity to vertical hyper-resistive diffusion that ejects ions in the -direction, analogous to radial kink mode diffusion producing in Section 3.2.1. Here is averaged over the energies of escaping ions, and is the flux width in the nose. The disk-averaged vertical velocity is dominated by disk ions giving and , which shows that it is the escape of cosmic rays vertically that drives cosmic ray acceleration radially (as weak radial escape drives vertical acceleration in the Central Column). The in Equation (21b) adds to any residual due to kink modes in the nose as field lines bend radially to create the nose, giving a continuous while turns into . Equation (21c) relates in the nose to the kinetic conversion (acceleration) efficiency in Equation (5e), aside from due to kink acceleration, in Equation (17), which we neglect here.
3.3.2 DCLC Instability
The origin of the cosmic ray accelerator field is not yet specified. Strong diffusion giving a large would inevitably arise due to the Drift Cyclotron Loss Cone (DCLC) instability (Post & Rosenbluth, 1966), not included in MHD but known from fusion research and characteristic of non-Maxwellian ion distributions such as the runaway ion beam created by kink mode acceleration in the Central Column. This runaway beam injecting energetic ions into the nose serves as the first stage of cosmic ray acceleration, DCLC in the nose being the second stage. Both the two-stream instability mentioned at the end of Section 3.2 and the DCLC instability are caused by ion excitation of electrostatic electron waves satisfying the following dispersion relation derived in Appendix B, applicable to DCLC as discussed in Appendix B.2. We obtain:
[TABLE]
where and is independent of electron mass (see Appendix B and C); and are components parallel and perpendicular to ; and is the ion term derived in Appendix B. In the Central Column, initially the dominant electron waves are plasma oscillations with on the right hand side. This is the two-stream instability that we found not to affect ions very much, in Section 3.2. But as ions enter the nose, new conditions allow the ions to excite electron “drift waves” giving DCLC instability resonant at the ion cyclotron frequency. Equation (22b) describes drift waves, with . Drift waves are caused by the electron density gradient given in the nose by . Why and when drift waves produce DCLC instability can be understood as follows.
There are two conditions yielding DCLC instability. First, DCLC is a self-driven cyclotron in which the oscillating drift wave acts as the cyclotron accelerator but the energy comes from the ions themselves. It is for this reason that the drift wave frequency must resonate with ion rotation. The second, more important, condition concerns when ions can transfer energy to the drift wave. This aspect of DCLC physics is well-known from mirror fusion devices, where DCLC is the remaining electrostatic instability when Landau damping stabilizes all modes with finite along . The DCLC with requires that the distribution function averaged over be non-Maxwellian (Paper III), as is the case in mirror devices due to the ejection of low energy ions by an electrostatic potential needed to confine the electrons (Fowler, 1981). The potential creates a hole in momentum space around , yielding free energy producing DCLC turbulence. An analogous situation occurs when ions accelerated parallel to B in the Central Column encounter the nose. The DCLC instability is a mode with driven by a hole in the averaged distribution (Paper III). As we show below, ions following bending field lines into the nose experience a centrifugal force that causes them to acquire a minimum , giving then a hole in around that can drive DCLC.
The DCLC instability occurs for a sufficiently large value of , found by choosing and to minimize the value required to satisfy Equation (22b) for real . The derivation in the relativistic limit is reviewed in Appendix B.2, giving the same criterion for the onset of instability as the well-documented result for the non-relativistic case if we simply replace rest masses by relativistic masses, giving a relativistic ion plasma frequency and a relativistic ion cyclotron frequency . Then the DCLC instability occurs if (Fowler, 1981; Post & Rosenbluth, 1966).
We rewrite the DCLC instability condition by eliminating using with dominant field , and by eliminating using from Equation (9f), giving for relativistic current carriers with for return flux width . Substituting these quantities into the DCLC instability condition, we obtain after some algebra:
[TABLE]
where we use from Equation (9f) with . Substituting Equation (23b) into (23a) gives Equation (23c), with as ions enter the nose.
Equation (23b), with in Equation (20b), is derived from the perpendicular momentum equation , as follows. Since the two-stream instability does not scatter ions very much, the orbital spin velocity in Equation (23b) is determined by balancing the centrifugal force, , due to magnetic curvature radius , against the restraining magnetic force, with . This ignores the ion synchrotron radiative force, which is much weaker, so that in Equation (23b) is always maintained. This is the minimum required for ion confinement in a twisting magnetic field, giving then the hole in momentum space causing DCLC instability to occur.
3.3.3 Onset of DCLC Instability
Initially, the low carried forward from the Central Column preserves a force free field as ions enter the nose, giving as field lines bend radially in the nose with flux width , obtained from force balance , which gives using and flux conservation . Also, initially the DCLC instability condition is not satisfied at the entry to the nose. However, for parameters above and numbers from Section 2.8, even with no further acceleration by kink modes, the condition for DCLC instability would already be satisfied at a radius where the DCLC condition in Equation (23c) is first satisfied.
We calculate the radius where DCLC commences as follows. We take to be at the margin of DCLC instability, satisfying Equation (23c) with the equality sign with and in the nose. We obtain:
[TABLE]
Again, is in units of solar masses. Equation (24a) uses Equation (23b), and Equation (24b) is Equation (23c) with the equality and . Equation (24c) comes from Equation (24b) using Equation (24a) with . Numerical values use Equation (20b). These results follow from the fact that, even though as ions flow into the nose scaling gives constant at its value in the Central Column until further acceleration occurs, the cyclotron frequency giving that must eventually exceed a fixed in the DCLC instability condition, Equation (23c).
3.3.4 Ion Acceleration to UHECR Energies
After the onset of DCLC instability at , two things happen. First, acceleration causes the ion energy to begin to increase. Secondly, DCLC scattering of ions gives in about one cyclotron period, too fast for ion synchrotron radiation to prevent this as it did in the Central Column. Then, using from Equation (9g), the relativistic ion Larmor radius in the nose becomes with an energy distribution f(E) due to downward scattering by DCLC. (Notation: note that here eV means multiply the accelerator voltage V by the electron charge, not to be confused with eV as a unit of energy.)
At the onset of DCLC, Equation (24d) shows that ions are still well confined within the initial flux width of the force free field. As the energy increases by DCLC acceleration, ions would remain confined only if the flux width widens to satisfy for all confined ions. The flux width also determines a pressure parameter giving the vertical pressure balance. We calculate in three steps:
[TABLE]
Equation (25b) uses . Dividing Equation (25b) by gives Equation (25c). The far right hand side of Equation (25c) uses flux conservation to write giving .
Thus the onset of DCLC changes the dynamics in the following sequence: DCLC acceleration sets in at ; the Larmor radius increases to equal , the initial flux width in the nose; and any further increase in energy could cause the ions to escape. But, in trying to escape, ions attached to field lines begin to spread out the flux, causing the flux width to expand so as to contain the most energetic ions. And as the flux expands, by Equation (25c). Thus, in fairly short order, system variables evolve to:
[TABLE]
where is the energy as ions emerge from the Central Column if we neglect kink mode acceleration in the nose, and is the energy of ions that enter the nose directly from the Central Column, as opposed to cold ions recycling from the ambient as discussed in Section 4.4.
At , field lines begin to untwist, giving so that the current carried dominantly by can no longer flow parallel to field lines. Yet this current must be maintained to satisfy our mean-field MHD equilibrium, since the displacement current can be ignored for this slowly evolving field, giving so that the poloidal current is continuous across the nose. Since the nose ion current must flow perpendicular to the dominant field , as field lines untwist, a hyper-resistive diffusive transport rate is required to carry the current, giving:
[TABLE]
That a net ion current results in Equation (27a) arises from the fact that by DCLC turbulence represents a random walk in momentum due to ion scattering analogous to collisional scattering. Because DCLC is resonant at the ion cyclotron frequency, scattering by DCLC affects ions only but does not affect electrons. That drift waves excited by DCLC do affect both ions and electrons plays a role in vertical transport in and out of the nose, in Section 4.4. Finally, Equation (27c) takes note of the transition from current density in the Central Column where current flows along field lines, to a reduction by a factor when finally diffusion must transport current perpendicular to , where we anticipate from Equation (29e) below.
It remains to verify that ion radiation can be neglected in the nose. We rewrite the momentum Equation (14a) to apply to the nose, giving:
[TABLE]
where we approximate and the radiation term, from Equation (18), represents both synchrotron radiation () and curvature radiation (). Before DCLC sets in, in Equation (23b) gives synchrotron radiation equal to curvature radiation. Then, even with , integrating Equation (28) gives for any before DCLC begins. When DCLC does set in, the dominant synchrotron radiation with is still negligible compared to DCLC acceleration approximated as , giving near , decreasing as increases.
3.4 DCLC Quasi-Linear Transport, Ion Energy Distribution
We now try to calculate the likely nonlinear outcome of DCLC instability in producing the accelerated ion distribution. The non-linear development of DCLC instability was shown to explain momentum transport in the 2XIIB mirror experiment (Berk & Stewart, 1977). Later Smith & Cohen (1983) applied Kaufman’s formal method (Kaufman, 1972) to derive an exact quasi-linear DCLC transport equation in action-angle space, a method that can be extended to give a Fokker-Planck equation including both diffusion and friction (Fowler & Gatto, 2007). For relativistic DCLC, the important action variables are the ion magnetic moment , and the radial canonical momentum with and vector potential giving . Here we will approximate the formal results in the two-dimensional and space to proceed directly to a 3D quasi-linear equation for transport in , and , as follows.
Here and hereafter, we drop the superscript so that, for example, will mean hyper-resistive diffusion in the radial component of momentum, and so on. We obtain:
[TABLE]
Here the source represents input from the Central Column (at energy ) plus recycling from the ambient discussed in Section 4.3. In Equation (29a), the first term on the right hand side represents ion acceleration in the nose. The other terms describe diffusion during acceleration. Equation (29b) applies the transform to variables , and . For DCLC, the formal theory shows that all transport processes arise from the electrostatic perturbation giving the relationships shown in Equations (29c), (29d) and (29e). Though and turn out to be of similar magnitude, physically derived from is analogous to a random walk by ion “collisions” due to cyclotron resonance, while includes also random walk by DCLC drift-wave motion. For strong instability with growth rate , both processes give a correlation time with the ion Larmor radius as the step size. That is the step size for motion follows from .
We approximate the solution of Equation (29b) by:
[TABLE]
We will return to the meaning of in Section 4.4.
We obtain Equation (30a) by setting the first term in Equation (29b) equal to zero, integrating and factoring out the ion density . We account for by taking the lower integration limit , to account for rapid acceleration of recycling ambient ions up to the energy of Central Column ions entering the nose. And we have set the spatial diffusion terms in Equation (29b) equal to zero. In Equation (30b), we first write for relativistic energies and then apply Equations (21b) and (29d), giving using also Equation (26a) with using in Equation (9g). To determine , we note that DCLC momentum diffusion only spreads energies downward (to fill the “loss cone”) so that ion energies at radius do not exceed , leading us to set to give . This downward spread of energies to achieve a state of marginal stability has been well documented in fusion mirror devices (Baldwin, 1977; Smith & Cohen, 1983); but inhibits downward energy spread in jets. Note that the final result for only depends on ratios of transport quantities, not their absolute value.
4 COSMIC RAY ENERGY SPECTRUM, INTENSITY ON EARTH
In this Section, we use the transport equation in Section 3.4 to compare our model with observations of UHE cosmic rays.
4.1 Cosmic Ray Energy Spectrum
In our model, cosmic rays are ejected via the term in Equation (29d). That this gives an energy spectrum similar to observations can be seen as follows.
Since UHE cosmic rays travel great distances, we treat AGN jet/radiolobes as point sources described by integrating the term in Equation (29b) over the entire volume of the nose, giving an energy spectrum . As several steps are required, we display them as follows:
[TABLE]
with given by Equation (26b).
In Equation (31b), we approximate together with a constant that we will use to guarantee energy conservation, as discussed below. The lower integration limit is defined by Equation (31f), chosen to eliminate values of not accessible by acceleration by . In Equation (31b), is given by Equation (31e) taken from Equation (30c), with in Equation (9f) with area for flux width in the nose, as in Section 3.3.1. Taking follows since DCLC only drives ion current to maintain current continuity around the loop. In Equation (31c), we first integrate in to cancel one factor , then substitute from Equation (26a) in remaining factors. In Equation (31d), we set using . While () can be order unity as discussed in Section 3.4, we choose to absorb () into giving here and hereafter the new fitting parameter in Equation (31g).
To see that Equation (31d) yields an approximate power law, we note that depends on and change variables from to , using:
[TABLE]
where volts is the full dynamo voltage from Equation (9g) with numbers in Section 2.8. We now apply this transformation to Equation (31d) with as the upper limit of integration and from Equation (30c). We obtain:
[TABLE]
The factor is given by:
[TABLE]
The dominant scaling factor arises naturally from diffusion across a flux width and the normalization of together with the change of variable whereby an integration in from becomes an integration in from .
4.1.1 The Fitting Parameter
We determine the fitting parameter by energy conservation. We delay discussion of the physical significance of to Section 4.4.
Neglecting the small power dissipated in the Central Column (1% by Equation (17)), the conservation of the jet power escaping as cosmic ray ions is given by:
[TABLE]
[TABLE]
where we have approximated
We estimate the conversion efficiency as follows. We approximate by an average value which guarantees that across the nose. This same approximation can be used to estimate the efficiency for converting electromagnetic energy into kinetic energy, from Equation (26b) with in Equation (24c). We obtain:
[TABLE]
4.2 Minimum Energy
As we noted at the end of Section 3, the fact that DCLC is driven by a hole in momentum space produces the peculiar result that DCLC momentum diffusion is not isotropic but mainly tries to fill this hole. At any radius , this limits the energy distribution to , where adjusts to maintain marginal stability and ( is the disk voltage) at .
Initially , the dominant energy of ions entering the nose from the Central Column, eV by Equation (19b). DCLC commences where the ion Larmor radius at energy satisfies the DCLC onset condition at marginal stability, given by the equality in Equation (23c). That is, DCLC onset does not determine . Rather, determines the radius of onset (giving in Equation (24c)). There is no DCLC turbulence at . On the other hand, as DCLC builds up the plasma pressure parameter , the field at is depressed, causing the Larmor radius to increase. In response, DCLC momentum diffusion at pushes proton energies downward just enough to re-establish marginal stability. As DCLC growth causes the pressure parameter to grow by Equation (25c), the field inside the nose decreases as , giving a lower DCLC threshold energy . We obtain:
[TABLE]
Here we have applied in Equation (25c) with in Equation (26a) at ; then introducing magnitudes from Equation (24c) and . More accurate determinations of DCLC marginal stability and the role of an additional Alfven Ion Cyclotron (AIC) instability occurring at high are discussed in Fowler (1981) using non-relativistic theory, giving with relativistic masses the same results as relativistic DCLC theory, as discussed in Appendix B2.
4.3 Calibrating the Model to Observations
We compare our with Figure 1 in Cronin (1999). This curve translates the count rate on Earth to an intensity proportional to our . As Cronin notes, like our , the basic scaling is , then reducing to at higher energies, perhaps for other reasons mentioned by Cronin. Cronin conveniently relates the plotted intensity to count rate above a specified energy giving for . This gave km2 year) for cosmic rays above about eV. The same scaling correctly predicted results giving km2 century) above eV for extra-galactic sources detected at the Pierre Auger Observatory (Pierre Auger Collaboration, 2007, 2014, 2017). Applying Cronin’s scaling to our model gives for a single AGN source:
[TABLE]
where is the ion flow rate for all energies , and Equation (38b) gives and from Equations (8b), (8c) and (9d). That implies both electron and ion recycling to maintain charge neutrality and a net disk current as required by . Recycling with the ambient is discussed in Section 4.4.
In terms of , data points on Cronin’s curve for AGNs at a distance correspond to a flux on Earth equal to To compare with the Pierre Auger result above, we take eV (Pierre Auger Collaboration (2007)). Summing over sources gives:
[TABLE]
The range in in Equation (31b) corresponds mainly to the range in , from the closest to the most distant AGNs accessible at the Pierre Auger Observatory (470 AGNs within a 240 Mlyr range, Pierre Auger Collaboration 2007). As examples, the lower limit in Equation (31b) would correspond to a single nearby source (, = 10 MPc) and the upper limit corresponds to several distant sources (, = 100 MPc).
How in our model can be consistent with Equation (40b) concerns density profile adjustments to maintain marginal stability to DCLC. We employed marginal stability to determine the onset of DCLC by taking the equality in Equation (23c) for ions of energy entering the nose from the Central Column. The onset of DCLC causes the plasma pressure parameter to grow.
As grows, the field inside the nose decreases as (, giving a lower DCLC threshold energy given by Equation (40c). Equation (40c) could accommodate Equation (31b), depending on the uncertain value of near the DCLC onset at . More accurate determinations of DCLC marginal stability and the role of an additional Alfven Ion Cyclotron (AIC) instability occurring at high are discussed in Fowler (1981) using non-relativistic theory, giving with relativistic masses the same results as relativistic theory, as in the case of DCLC (see Appendix B.2).
4.4 Recycling with the Ambient
The small value of in Equation (36) concerns recycling of protons as the jet nose pushes against the ambient. Recycling is given by:
[TABLE]
In Equation (41a), we use from Equation (36), giving a result independent of the power law exponent . Equation (41b) gives the distribution of cosmic ray power across the radial profile of the jet nose, with lowest energy ions nearest to the tip of the jet.
We note that the necessity for recycling is characteristic of a decreasing power law . An excess of ions beyond that needed to conduct the jet current results for any specified total power applied to the conservation of energy (to calibrate ), then applied to the conservation of particles (giving . Our model makes this quantitative by showing that the power law is characteristic of a single AGN source whose black hole mass , together with angular momentum conservation, determines all jet parameters as shown in Section 2.8. Then the total emission of extra-galactic cosmic rays includes some reaching Earth, some not.
We note also that the ambient can sustain recycling and recycling does not appear to change our model parameters. Physically, recycling is part of DCLC vertical transport, whereby adjustments in the Vlasov distribution functions for hot ions, ambient ions and electrons do in detail what adjusting does in our model. Thus recycling commences where DCLC commences, at in Equation (24c). Because of the low jet density in Equation (9f), snow-plowing an ambient density over an annular width at would be sufficient to accommodate the result in Equation (41a).The jet density in Equation (9f) already includes recycling, any excess being expelled vertically by DCLC diffusion to maintain just the ion density needed to conduct the current sustained by radial DCLC diffusion . Nor would recycling add much to ambient pressure, the ratio being using and eV, with .
5 ELECTRON RADIATION, PAIRS, FLARES, NEUTRINOS
This Section is mainly concerned with what we can say concretely about electron radiation in quasi-steady state. We will focus on electron synchrotron radiation extending to microwaves, which for decades has been and continues to be the main evidence for the jet/radiolobe structures that we claim also create the UHE cosmic rays reaching Earth.
We will neglect weak DCLC acceleration of electrons in the nose and focus on the Central Column, where electrons are accelerated to TeV’s, by Equation (19b). With steady state in mind, we omit discussion of transient flares, perhaps byproducts of kink mode turbulence. We neglect dilution of ions by positrons which, by two-stream interaction with their electron partners, cannot contribute to the jet current. Finally, we omit discussion of neutrinos, but note that the essentially mono-energetic proton beam produced by kink mode acceleration in the Central Column has associated with it lower energy protons following wandering field lines. As is mentioned in Section 5.4, these wandering protons accelerated mainly when they pass through the Central Column can spread proton energies downward, from eV to TeVs, while protons accelerated near where field line curvature is small would add a tail up to eV by Equation (17).
5.1 Electron Synchrotron Model
The two main features of our synchrotron model are: (1) angular momentum conservation in Equation (6c) which serves to project the magnetic field pattern at the disk all along the jet so that Central Column parameters at the disk determine synchrotron wavelengths; and (2) plasma turbulence which spreads this field projection over an ever-widening cone like the dashed path in Figure 1, finally expanding into giant radiolobes like the closed flux in that figure.
As noted earlier, what causes some field lines to wander away from the Central Column is MHD kink instability. As is explained in detail in Paper II, there are two distinct kinds of MHD kink modes driven by the jet current, confirmed in MHD simulations mentioned in Section 2.7. In Paper II, we derived the two kinds of kink modes from the MHD energy principle giving the free energy for magnetic perturbations of the form where subscript 1 denotes perturbations and is a displacement of magnetic field lines.
Internal kink modes with localized to the Central Column create the well-understood kink mode accelerator in Section 3. Internal kink modes are characterized by magnetic resonances in . External kink modes of interest here are those with small enough to avoid these resonances, requiring for shown in Figure 1, also dominated by describing a rigid body motion of the entire magnetic structure. The other point made in Paper II is that, because accretion disk power is concentrated near , should be thought of as arising at , , growing not in time but in along the jet length; and relativistic effects give “secular,” not exponential, growth, giving .
Thus we come to a picture of the synchrotron radiator localized in the Central Column together with an array of magnetic field lines emerging from the disk as the dotted cone in Figure 1, “wandering” in the sense of a twist of the 2D structure as a whole, first near the black hole but growing radially along the jet when viewed in 3D.
We will be able to predict the synchrotron power and dominant wavelengths without detailed knowledge of the wandering field line structure. The power is known from the calculation of the voltage drop in Equation (17). The magnetic field strength determining synchrotron wavelengths can be determined by angular momentum conservation, as follows.
Angular momentum projects the accretion disk field forward, giving synchrotron radiation at approximately the same magnetic field strength all along the conical field line structure. This can be seen by the following extension of Equation (6c) into the jet:
[TABLE]
Here by Equations (8d) and (9d). We extend the meaning of in Equation (42a) to include paths along wandering field lines, averaged over as indicated by giving correlated perturbations on the right hand side, omitting small mean field contribution at , giving the Diffuse Pinch in Figure 1. Equation (42b) is the MHD perturbation for wandering field lines at derived from Faraday’s Law, as in deriving above. Applying Equation (42b) to Equation (42a) gives Equation (42c). We set as is typically observed in simulations, still obeying at all r.
The path integral giving Equation (42a) is also the path of synchrotron-emitting electrons following field lines that leave the Central Column. We conclude that, although the mean field falls off as outside , the total field wherever accelerated electrons flow within the dashed cone in Figure 1, yielding the wavelength distribution in Section 5.3. Poloidal flux is properly conserved by Fourier expansion, giving the mean field that is uniquely the Grad-Shafranov solution in Figure 1 containing all of the flux (see Paper II). Also of order was confirmed in the tower-like MHD simulations in Carey, et al. (2011), giving in quasi-steady state. We have ignored the build up of flux-amplifying closed structures, as in certain fusion devices, negligible for astrophysical jets in which helicity injection by the accretion disk is mainly consumed in lengthening the jet (Fowler, et al., 2009).
5.2 Synchrotron Luminosity
The synchrotron radiation described in Section 5.1 dissipates a large part of the electron energy produced by kink mode acceleration in the Central Column. Electron acceleration is equal and opposite to ion acceleration in Section 3.2, giving as the synchrotron power, or luminosity, in good agreement with observations (Krolik, 1999). The synchrotron power is given by and in Equation (8), yielding:
[TABLE]
where the inequality accounts for non-synchrotron electron radiation (Li & Kusunose, 2000). This luminosity is consistent with observations of powerful AGN jets (Krolik, 1999).
5.3 Synchrotron Wavelengths
As noted above, a new feature of our model is that, both in the Central Column and in the surrounding structure (the cone and radiolobe regions in Figure 1), synchrotron radiation should always be calculated using the magnetic field Gauss at near the black hole. Then, wherever they occur, the frequency and wavelength for synchrotron emissions at the nth harmonic are given by:
[TABLE]
with numbers from Section 2.8.
For a given electron energy , radiation intensity is spread over a range of frequencies around the fundamental (Jackson 1998, Equation 14.31), the total radiation for all harmonics being that given by Equation (18) applied to electron cyclotron orbits in a magnetic field. Observed radiation can be found by integrating the radiation at a fixed . An example of this procedure for an assumed energy distribution is given in Li & Kusunose (2000).
In principle, the electron transport model in Appendix C.3 could yield a correct distribution, but this has proved more difficult than the corresponding ion calculation in Section 3.4; and our techniques do not yet treat flares and hot spots except to note that transients consume less power than steady state (Fowler, et al., 2009). That this does not matter for our main goal of explaining UHE cosmic rays follows from the limited power in all electron radiation, by Equation (43).
We do predict an upper limit on electron energies TeVs due to curvature radiation in the Central Column, giving in Equation (19b) both for electrons and ions in the Central Column. Our estimate of lower ion energies outside the Central Column, in Section 5 .4, suggests the same for electrons, indicating the same maximum synchrotron wavelengths everywhere in the jet/radiolobe system (the upper limit in Equation (44b)) but the shortest wavelengths and highest frequencies in and near the bright Central Column with radius by Equation (9a).
5.4 Synchrotron Opening Angle
As is mentioned in Section 2.7, wandering field lines due to MHD kink modes distribute synchrotron radiation power over the dashed cone in Figure 1, though the axi-symmetric Central Column remains magnetically collimated. In Paper II, it is shown that relativistic ion acceleration on kinking field lines yields secular growth of the field line displacement , giving then an opening angle of the dashed cone in Figure 1 given by:
[TABLE]
giving in radians. Here we combine Equations (39) and (53) in Paper II, giving as the average Lorentz factor of ions wandering in and out of the kink accelerator zone concentrated in the Central Column. An ion with for the most extreme field line excursions is to be compared with a maximum of order in Equation (19b).
6 SUMMARY, COMPARISON TO PREVIOUS WORK
This is the final paper in a series of four papers showing that the production of ultra high energy cosmic rays occurs naturally when processes exist to create an efficient dynamo by accretion around a supermassive black hole inside an AGN. Section 2 and Appendix A argue that MRI is the likely process for such a dynamo.
Given the ubiquitous presence of seed magnetic fields sufficient to excite a persisting MRI dynamo, as discussed in Appendix A.3, our model predicts that all AGNs should produce UHE cosmic rays, typically for a jet lifetime yrs that we equate to the time for accretion to exhaust ambient mass around a black hole (accretion rate ). Why actually only of AGNs are the radioloud variety we associate with UHE cosmic rays may concern the year jet duration in relation to the 100-times-longer life of the Universe. The fact that is independent of black hole mass may imply a connection to the Eddington limit on luminosity (Beskin, 2010). Why some AGNs eject one-sided jets (up or down) should, by our model, require a one-sided dissipation of magnetic helicity (a product of poloidal and toroidal fluxes) since helicity conservation requires two jets as discussed in Appendix B of Paper I.
In Sections 3–5, we have shown that the magnetically-collimated jet predicted in Section 2 automatically becomes a particle accelerator that explains many features of UHE cosmic rays and the synchrotron radiation by which AGN jets are observed. A common theme is the use of hyper-resistive diffusion in a two-dimensional, mean-field Ohm’s Law to evaluate the consequences of 3D turbulence, both MHD and kinetic, often avoiding the need to know turbulence magnitudes through Onsager-like relationships among observables (for example, in deriving the cosmic ray energy spectrum in Equation (33b)). Kinetic effects play key roles: to sustain a coronal short-circuit current for the duration of the initial conical jet (Equation (4a)); the two-stream instability, which allows a relativistic kink-mode accelerator to exist in the Central Column (Equation (13c)); and the DCLC instability that drives cosmic ray acceleration in the jet nose (Section 3.3.3).
Our approach differs from conventional wisdom in that, given the dominance of magnetic fields above disks if MRI is the creator of a dynamo, we regard accretion disks, jets and the jet nose highlighted here as a single magnetic system evolving according to the hyper-resistive Ohm’s Law given in Equation (2a). In this model, the bright jets ejected by AGNs serve mainly as a magnetically-guided conduit of power from the black hole dynamo to the nose end of the jet. That the jet is magnetically-collimated serves both to stabilize this linkage and also to limit energy loss by the synchrotron radiation by which jet/radolobe structures are observed, thus leaving most of the jet power to accelerate and eject cosmic ray ions in the nose.
We recognize that magnetic collimation of jets over Mpc dimensions is controversial. Three main schools of thought have contributed to understanding magnetized accretion disks and the jets they produce. One school, stimulated by the centrifugal ejection model of Blandford & Payne (1982), has focused on magnetic jets as accelerated “winds” (Beskin, 2010), either at the disk or by conversion of magnetic energy to kinetic energy early in the jet history. Another school has focused on simulations from first principles (GRMHD codes, McKinney, et al. (2012)). A third school, including our work, has focused on the analogy between AGN jets and magnetically-collimated jets produced in the laboratory (Colgate & Li, 2004; Lovelace, 1976; Zhai et al., 2014).
By using simple models to identify critical physics, we believe that, as summarized in Section 2, we have closed the loop among these different schools of thought, with the conclusion that well-understood laboratory jets are the correct model for AGN jets on astrophysical timescales, but only after an initial conical jet launched by rotation has evolved to the magnetically-collimated jets required for our model of cosmic ray acceleration. Interestingly, the same rotation that makes accretion jets different from laboratory experiments at launch eventually causes a jet to become magnetically-collimated, already occurring at 1 of its final length, as estimated in Section 2.7. Other differences from previous work are:
- (1)
In Paper I: The concept of jets generated by whirling magnetic flux frozen into the accretion disk (Frank, et al., 2002; Lynden-Bell, 1996) is replaced by non-ideal MHD in which MRI-driven hyper-resistivity is the mechanism of jet ejection. As in magnetized jets created in the laboratory, these jets are ultimately ejected vertically by an electrostatic sheath, with just enough ion mass to carry the current, as opposed to continuing centrifugal ejection (Blandford & Payne, 1982). The electrostatic sheath is discussed in Paper II, Appendix A.2. 2. (2)
In Paper II: We show that magnetic collimation produces MHD kink instability at short-wavelengths known to accelerate electrons and ions in the laboratory without destroying collimation, together with very long wavelength kinks producing wandering field lines carrying the power observed as synchrotron radiation by electrons in giant radiolobes and as a bright cone with opening angle 0.01 radians, consistent with our model. 3. (3)
In Paper III and this paper: We show that, in the nose-end of the jet, flaring field lines both suppress electromagnetic radiation and create the conditions for kinetic instability that accelerates ions to energies equal to the dynamo voltage. Unlike transient mechanisms along the lines of Fermi’s initial explanation of cosmic rays, we find that AGN jets become steady-state accelerators that can convert about half of the jet power into the kinetic energy of UHECRs. Our acceleration mechanisms are based on known plasma processes that can be further explored in the laboratory as a guide to astrophysical simulations and observations (see Paper III).
Our model predictions are consistent with nine observables, outlined in Table 1 with references to the text in this paper. Additional ways to test the model in future work are:
- (1)
In GRMHD and Relativistic MHD simulations of accretion disk jets, apply a diagnostic similar to Figure 5 in order to measure directly the short-circuit current predicted in Sections 2.2 and 2.6 as the mechanism of transition from a conical jet to a magnetic tower. 2. (2)
New laboratory experiments on ion acceleration outlined in Paper III. 3. (3)
Non-linear kinetic simulations of two-stream instability in a relativistic current-carrying column to verify the persistence of a kink-mode accelerator in the relativistic regime.
ACKNOWLEDGMENTS
We wish to express our deep appreciation for the intellectual foresight by our colleague Stirling Colgate that motivated this paper and Papers I and II coauthored with him, still in progress when Stirling passed away December 1, 2013. We thank Alexander Tchekhovskoy for suggesting the short-circuit current and Paul Bellan for its theoretical formulation in Equation (2b). We thank Joseph McClenaghan for the use of Figure 6. We also thank Mitch Negus for his help in preparing the manuscript. HL gratefully acknowledges the support of the U. S. Department of Energy Office of Science, the LANL/LDRD program and NASA/ATP program for this work. RA acknowledges the California Alliance Postdoctoral Fellowship Program, the Simons Foundation and the John Templeton Foundation.
Table 1. Model Predictions
Jet Properties
Jet length Equation (10c) 2. 2.
Radiolobe radius Equation (10d)
Synchrotron radiation 3. 3.
Synchrotron luminosity Equation (43) 4. 4.
Synchrotron wavelength cm Equation (44b) 5. 5.
Synchrotron opening angle Equation (45) 6. 6.
Maximum Electron Energy TeV Equation (19b)
UHR Cosmic Rays 7. 7.
Maximum UHECR energy Equation (38b) 8. 8.
UHECR intensity on Earth Equations (39a,b) 9. 9.
UHECR energy spectrum Equation (33b)
APPENDIX A: MRI-DRIVEN ACCRETION DISK DYNAMOS
The evolution of conical jets to magnetic towers in Section 2 depends on MHD hyper-resistivity inside the disk, given by (Fowler & Gatto, 2007). A finite mean-field emerges if correlations between and give a finite value when is averaged over fluctuations in the same way that gives the mean fields in Section 2, indicated here by . That conical jets evolve to the magnetic towers of our model follows from Equation (3d) in the main text, requiring together with vertical accretion velocity and justified by requiring a field line tilt , at the threshold of the vertical force balance of gravity against centrifugal ejection of mass (Paper I, Appendix A; Blandford & Payne (1982)).
The proof that is derived from:
[TABLE]
In Paper I, Appendix A, follows from assuming equipartition of and perturbations in the non-linear steady state . Then, given proper correlations, the numerator and denominator are equal without knowledge of and . As was stated in Paper I, with these assumptions Equation (A1b) is valid for the non-linear steady state for any kind of MHD-like fluctuations.
A.1: Hyper-Resistivity by MRI
We now evaluate Equation (A1a) using and from the linearized equations giving MRI instability in Balbus & Hawley (1998)—hereafter BH98. In applying linearized theory, we assume that MRI produces the non-linear correlations necessary to yield a finite , as kink modes are known to do as discussed in Section 3.2. Then linear theory provides an estimate of ratios of mode amplitudes. We obtain:
[TABLE]
Equation (A2b) is the MRI replacement for Equation (A1b), obtained as follows. We apply linearized MRI theory for the simplest case of plane waves in a uniform field, with growth rate (BH98, Equation (111)). As is common for ideal MHD, we write perturbations in terms of a vector field line displacement giving . Using this and obtained from the induction equation, Equation (A2c), we obtain Equation (A2d). Using this and Equation (A4d), we obtain Equations (A2a) and (A2b).
The approximation applies for the maximum MRI growth rate, , representing the linearized modes most likely to grow out of the noise. For these modes, with and (BH98, Equation (114)), we obtain exactly. We then assume that fluctuations giving during peak growth are also characteristic of the steady state (Kadomtsev, 1965).
A.2: Jet Evolution by MRI
Consider an initial condition — typical of GRMHD simulations — with a pre-existing black hole mass , a pre-existing finite poloidal magnetic field , but and an ambient environment that we take to be a hydrogen plasma at for an actual accretion disk. Then the dynamics inside the accretion disk is described by mass flow and the radial force equation and angular momentum equation, giving:
[TABLE]
where we have used Equation (A3a) to eliminate terms in , giving then Equation (A3b) on dividing by and Equation (A3c) on dividing by . We have omitted the viscosity in Eq. (A3c), which is small to order (Colgate, et al., 2014). Gravitational accretion begins for the few ions heading directly toward the black hole (the term in Equation (A3b)) and more importantly, for the much larger population near enough to the black hole so that gravity exceeds their centrifugal force .
In the limit that in Ohm’s Law in Equation (3a), accretion is advective, whereby accreting ions stuck to field lines compress the flux, possibly creating magnetic stagnation of accretion by the term on the right hand side of Equation (A3b) (true in some GRMHD simulations). Then, as arises by MRI, the terms dominate, giving a mainly Keplerian disk ejecting a conical jet that transitions to a magnetic tower as discussed in Section 2, where and Equation (6c) follow from integrating Equations (A3a) and (A3c) in over the disk half-height (Colgate, et al., 2014). The timescale for MRI evolution is typically . The coronal ratio characterizing magnetic towers is achieved as positive makes grow by Equation (3a) while negative halts advective compression of by Equation (3b).
A.3: Persistence of MRI Dynamos
For MRI as the dynamo process, the main issue is the persistence of as instability saturates. Given an adequate seed magnetic field, MRI can persist even as grows from zero to exceed (BH98, Section IV.F). To see this, let us revise the dispersion relation in BH98 (Section IV.B), by adding to the equilibrium and to the perturbation, giving with . We obtain:
[TABLE]
where is the growth rate.
Equation (A4a) is the linearized induction equation with . Equations (A4b) and (A4c) are momentum Equations (106) and (107) in BH98, obtained by inserting Equation (A4a) into . Including gives in Equation (A4b) with no effect in Equation (A4c).
Combining Equations (A4b) and (A4c) gives the dispersion relation:
[TABLE]
where for Keplerian rotation. For , Equation (A5a) reduces to BH98, Equation (111). Keeping gives instability with growth rate if the right hand side of Equation (A5a) is negative, requiring:
[TABLE]
with the limit on in BH98. Thus we see that MRI initiated by poloidal flux alone persists even as grows to exceed . Moreover, that adding also produces the coupling to perturbations in needed to give finite and in Equation (A2a). We get:
[TABLE]
The ratio of omitted electric force terms to magnetic terms can be shown to be of order at .
Concerned that MRI may not persist, Pariev & Colgate (2007) proposed an alternative mechanism whereby star-disk collisions provide an external means of sustaining dynamos, supported by simulations that maintain a dynamo even when MRI is suppressed (Pariev, et al., 2007). Star-disk collisions could be simulated in GRMHD codes by prescribed, random heat pulses.
Given any process that can sustain , any seed magnetic field around black hole can create an accretion-driven dynamo. Whatever the process, it is only necessary that embedded plasma be “magnetized,” meaning that the ion Larmor radius should be small on spatial coherence scales. For example, if an AGN accretion disk captures a single magnetized star, like our Sun, it might add a field Gauss, giving for a maximum gravitational energy of order GeV for hydrogen ions. Then a seed field coherent over dimensions would grow to astrophysical dimensions. Simple estimates show that even advection of primordial plasma can create a seed field, with density (Colgate & Li, 2004) attached by magnetization to the primordial field ( gauss, Ando & Kusenko (2016)), giving a dynamo if the field growth rate , already occurring at a radius about equal to in Figure 1b, as derived in Paper II, Equation (49). This suggests that, in actual disk conditions, MRI itself should be sufficient to sustain the dynamo.
A.4: Extending the Dynamo Inside the Central Column, Cowling’s Theorem
With the approximations above, strong MRI dynamo action exists across the Diffuse Pinch where is Keplerian and is finite, but the dynamo could die out toward the black hole where the breakdown of our Diffuse Pinch solution indicates a fall off of rotation below Keplerian, even though due to toroidal current outside remains large all the way to the black hole. That current inside the Central Column is essential to our model of UHECR acceleration is evident from Figure 4, showing that is of the asymptotic current re-entering the disk at a large radius. That the dynamo persists inside the Central Column is due to the fact that MRI is not local, but should be averaged in radius over the reciprocal not yet specified. A plausible value easily spans the gap from the Diffuse Pinch to the event horizon. Moreover, if rotation were to slow down too much for MRI to persist, a term omitted from Equation (6c) would build up rotation to restore MRI. Fortunately, details inside the Central Column are not important for kink stability analysis, as shown in Paper II, and our electric circuit representing the Central Column can include dynamo action at the event horizon (Blandford & Znajek, 1977; McDonald & Thorne, 1982).
We conclude that our estimates in Papers I and II are a reasonable guess, giving near the black hole (Equation (65b), Paper II):
[TABLE]
and, at the opposite extreme, growing poloidal flux near the O-point at , , where the poloidal field vanishes in Figure 1, given by (Equation (B1), Paper I):
[TABLE]
To defeat Cowling’s anti-dynamo theorem (Moffatt, 1978), it is shown in Papers I and II that negative by Equation (A2) becomes positive near , representing helicity flow from the Central Column to the O-point as in the creation of spheromaks in the laboratory.
APPENDIX B: KINETIC INSTABILITY IN AGN JETS; DCLC
In this Appendix, we review the relativistic derivation of Equation (22) giving, with all terms included, two-stream instability in the Central Column and also Drift Cyclotron Loss Cone (DCLC) instability in the nose that account for UHECR acceleration in our model.
B.1: Formalism
Both two-stream and DCLC dispersion relations used in this paper appear in the past literature. Here we show their common origin as ion excitation of electron waves, even at relativistic velocities. These are electrostatic modes with potential perturbation . We apply the relativistic Vlasov solution employing Fourier transforms and comment on Landau damping separately. The main features emerge in the “slab” (Cartesian) coordinate approximation, given by:
[TABLE]
where is the direction along , while lies in the flux surface, and is perpendicular to a flux surface (the direction of density gradients giving drift modes for DCLC).
The unperturbed distribution function is given by:
[TABLE]
In Equation (B2d), is the rest mass and , and also the total energy are relativistic constants of the motion in a uniform magnetic field, while is still constant and the magnetic moment is an adiabatic invariant in a non-uniform field. For electrons, dependence of on , giving the greatest effect of relativity, yields the two-stream instability, while dependence on adds spatial dependence giving the electron drift waves excited by DCLC. The mean-field density is defined by the normalization , and similarly the perturbation . We first ignore dependence, adequate for ions, giving . Then the perturbation is given by:
[TABLE]
We now integrate Equation (B3a) to obtain and from this the density perturbation . We follow Fowler (1981), Section VI. Several steps are involved. First, since is a Fourier transform, we relate at an earlier time to by a factor given by:
[TABLE]
Here sums over indices and from to . Introducing orbital information into Equation (B4a), following (Post & Rosenbluth, 1966), we obtain Equation (B4b), with Larmor spin angle and since and are constant along the orbit since and the energy are constants of the motion, and we take to be approximately uniform over the orbit. Introducing the identity Equation (B4e) into Equation (B4b) gives Equation (B4c), where anticipates the average in , giving in Equation (B4d). Changing variables to and and completing this integration gives the perturbed density, with normalization :
[TABLE]
Here ′ denotes quantities at time . The spin angle integration has already been performed in Equation (B5a) giving Equation (B5b), using Equation (B4d). Applying Equation (B5f) to replace gives this operation in the first line of Equation (B5c). Carrying out this operation and integrating on gives Equation (B5d), in which the leading term no longer involves an orbit integration (Post & Rosenbluth, 1966). Carrying out the remaining integrations on yields the familiar resonance denominators in Equation (B5e). Adding dependence on gives electron drift waves, yielding an additional term of the form:
[TABLE]
where . Integrating first over the orbit and then over momentum gives the drift contribution to the density perturbation. Integrating over momentum eliminates the first term. We obtain for the charge density:
[TABLE]
Equation (B7c) gives the drift term to be added to Equation (B5b). Integration on gives Equation (B7b) by analogy with Equation (B5c) (but missing the operation coming from a factor in that equation). Taking for electrons gives Equation (B7c), rewritten on the far right hand side as it usually appears in the literature (Post & Rosenbluth, 1966), though in fact this term does not contain mass so that there is no relativistic correction. Introducing (B7c) into Poisson’s Equation, , and dividing by gives the right hand side of Equation (22b) giving the electron drift wave contribution in the DCLC dispersion relation in Section 3.3.2. The electron term with in Equation (B5e) will give two-stream instability in Appendix C.
B.2: Relativistic DCLC Instability in the Nose
Growth of or a pre-existing adds ion cyclotron resonance terms to the left hand side of Equation (22a). Resonance gives large contributions from the ion Bessel functions in Equation (B5), but only if inside the plasma. That this cannot happen in the Central Column follows from:
[TABLE]
where the subscript [math] denotes rest mass quantities. We see that is not possible for any known black hole mass. In laboratory experiments, these ion cyclotron modes occurred when the non-relativistic (Post, 1981); or equivalently, when:
[TABLE]
Equation (B9) is an approximate condition for DCLC instability with pressure parameter in Equation (20a). This can be satisfied in the laboratory at low if , but in AGN jets it is satisfied only as ions enter the nose.
A DCLC dispersion relation like the non-relativistic form giving the DCLC instability threshold in Equation (23a) can be derived from Equation (B5e) applied to ions. The non-relativistic derivation, in Post & Rosenbluth (1966) reproduced in Fowler (1981), follows the same steps that led to Equation (B5e). Relativity would appear only in the mass defining the cyclotron frequency. The result, using also the drift term discussed above, is given by:
[TABLE]
Instability requires near representing the “hole” in perpendicular energy discussed in Section 3.3.2. It was found that, mathematically, a sufficient hole is present when the large Larmor orbits cause the -averaged distribution to be strongly positive while a Maxwellian in would give a negative value. Given a hole, Equation (B10a, b) is identical with the non-relativistic formulation, aside from swapping for in the integration, giving relativistic masses but otherwise identical with the DCLC instability in non-relativistic form with relativistic masses, as in Equation (23a).
APPENDIX C: ELECTRON SYNCHROTRON RADIATION, TWO-STREAM INSTABILITY
In this Appendix, we discuss synchrotron radiation, beginning with the dispersion relation for two-stream instability that converts MHD kink mode acceleration of electrons in the Central Column into orbital spin giving rise to synchrotron radiation, and concluding with the synchrotron energy spectrum.
C.1: Two-Stream Instability
To obtain the dispersion relation giving two-stream instability between counter-streaming ions and electrons in the Central Column, we set in Equation (B6) giving by Equation (B5e) with and we take the limit in all Bessel functions giving and for all . We assume ions are accelerated to velocity and electrons to velocity . Poisson’s equation becomes:
[TABLE]
In Equation (C1b), setting has eliminated terms with giving in Equation (B5e). Equation (C1c) follows on dividing by after applying the following result in Montgomery & Tidman (1964), Equations (10.51) and (10.52). Using and from Equation (B3d) with relativistic , we obtain:
[TABLE]
where Equation (C2b) is the surviving term in an integration by parts.
Combining Equations (C2b) and (C2c) gives the final two-stream dispersion relation, Equation (C1c), with the assumption that there is no acceleration perpendicular to B, giving the far right hand side of Equation (C2c).
We note that the ion term in Equation (C1c) is largest at resonance, suggesting that we take near the ion resonance (Goldston & Rutherford 1995, Section 23.4 ). Upon substituting and expanding, Equation (C1c) gives:
[TABLE]
where in Equation (C3e) is the unstable root. Here . In the Central Column, radiation-limited acceleration gives by Equation (18) including electron synchrotron radiation, giving then the ratio of rest masses in Equations (C3d) and (C3e) so that the only relativistic effect is the speed of light as the wave velocity.
C.2: Electron Scattering
Once initiated, the two-stream instability can scatter electrons to transfer parallel acceleration to synchrotron radiation, as in Section 5.1. To do so, the electron beam must excite waves with non-zero which adds terms to the dispersion relation in Equation (C1c). Keeping both the electron drift term from Equation (B7c) and potentially-resonant ion terms with finite , we obtain:
[TABLE]
Terms labeled “wave” generate the waves driven unstable by the terms labeled “drive,” with from Equation (19b). In the absence of scattering, becomes a delta-function giving the result in Equation (C4a). Ions do not scatter, as discussed in Section 3.2, giving in the drive term in Equation (C4b) the same ion term as that in Equation (C2a), together with a new term taking into account cyclotron resonance. This ion cyclotron term comes from Equation (B5e) but dropping the term of order compared to on the left hand side of Equation (C4a).
Lastly, we include in the drive term an imaginary contribution from the electron Landau pole, whereas the integral in Equation (C4a) is to be interpreted as the principal part. The factor in the Landau term arises from changing variables from to . The Landau pole gives damping or inverse damping depending on signs. The magnitude of the Landau damping term in Equation (C4d), comparable to other electron terms and much larger than ion terms away from resonance, indicates that Landau damping would cause electron momentum diffusion to stop short of the ion resonance. Resonance extends between positive , giving, by Equation (3e), the following allowed spread of electron velocities with no Landau damping:
[TABLE]
Thus, as noted in Section 3.2, a uniform spread in electron momentum gives negligible net velocity and negligible electron current while avoiding electron resonance giving Landau damping. That Landau damping exists but ceases if there are few resonant electrons has been demonstrated conclusively in laboratory experiments in the low density regime relevant to AGN jets (Post, 1981).
A necessary condition is that (surrogate for in cylindrical coordinates) fit inside the Central Column. Keeping only the wave terms in Equation (C4a), the dispersion relation in Equation (22) reduces to:
[TABLE]
The quantity in Equation (C6b) approximates in the principal-part integration in Equation (C6a) as a constant over a momentum spread limited by curvature radiation, in Equation (19b). The strong weighting at small gives the result shown with . Examining these results shows that waves fit with and .
C.3: Electron Transport in the Central Column
Synchrotron radiation mainly arises due to acceleration of electrons by MHD kink modes in the Central Column, as discussed in Section 3.1. Electron acceleration is equal and opposite to ion acceleration in Section 3.2. The counter-streaming electrons and ions excite the two-stream instabilities that scatter electrons to produce also electron synchrotron radiation. We describe these electron processes by a quasi-linear transport equation analogous to Equation (29b) for ions in the nose, now omitting spatial derivatives but keeping ( above) as an approximate constant of the motion and transforming into . We obtain:
[TABLE]
where determines the energy spread by two-stream instability limited by curvature radiation giving Equation (19b) and describes scattering from to that yields synchrotron radiation. In Equation (C7c), is electron synchrotron radiation (Equation (18), ) with electron Larmor radius and ) is constant.
Qualitatively, assuming factorability gives by setting the bracket and by setting the bracket , with coupling via . Strong coupling gives relatively flat in the domain in Equation (C5), with a tail terminating at TeV’s due to ; and perhaps best approximated by a power law as in Li & Kusunose (2000). By contrast, dominant with yields with determined by the transport coefficients.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1al-Karkhy, et al. (1993) al-Karkhy, A., Browning, P. K., Cunningham, G., Gee, S. I. & Rusbridge, M. G. 1993, Phys. Rev. Letts. 70, 1814
- 2Ando & Kusenko (2016) Ando, S. & Kusenko, A. 2016, Ap J Letters
- 3Balbus & Hawley (1998) Balbus, S. A. & Hawley, J. F. 1998, Rev. Mod. Phys. 70, 1
- 4Baldwin (1977) Baldwin, D. E. 1977, Rev. Mod. Phys. 49, 317
- 5Begelman, et al. (1984) Begelman, M. C., Blandford, R. D. & Rees, M. J. 1984, Rv MP, 56, 255
- 6Berk & Stewart (1977) Berk, H. L. & Stewart, J. J. 1977, Phys. Fluids 20, 1080
- 7Beskin (2010) Beskin, V. S. 2010, Phy U, 53, 1199
- 8Bierman (1997) Bierman, P. 1997, J. Phys. G 23, 1