Ultra High Energy Cosmic Ray Source Models: Successes, Challenges and General Predictions

In order to accelerate a proton to an energy $\sim 8\,{\rm EeV}$, it must follow a trajectory where $\int d{\bf r}\cdot{\bf E}$ exceeds $\sim 8$ EV. The simplest way to imagine this happening is in the context of a large scale magnetic field. Faraday’s law guarantees that if the magnetic field varies with time, then there will be an EMF and an opportunity for particle acceleration. However, even a stationary flow can produce large, accelerating, electric field when there are appropriate boundary conditions. One such configuration is the unipolar inductor (Weber, ). A spinning, conducting body is envisaged to be endowed with an axisymmetric magnetic field distribution. In the simplest case, the body rotates with angular velocity $\Omega$.

If sufficient plasma is continuously produced to sustain an electrical current, without having dynamical relevance, the electromagnetic field in the “magnetosphere” will be force-free. Imposing a boundary condition at the neutron star surface leads to an electric field ${\bf E}=-({\bf\Omega}\times{\bf r})\times{\bf B}$ which implies that ${\bf E}\cdot{\bf B}=0$. The magnetic flux surfaces will be equipotential and isorotational and the potential difference between two adjacent flux surfaces containing magnetic flux $d\Phi$ is $dV=\Omega d\Phi/2\pi$. There will also be an electrical current $I(\Phi)$ flowing within a flux surface labeled by $\Phi$ which is associated with a toroidal component of magnetic field and an outward Poynting flow of electromagnetic energy, $IdV$. This Poynting flux extracts rotational energy efficiently from the neutron star with very little dissipation in the neutron star. It can be radiated at great distance, for example in a Pulsar Wind Nebula such as the Crab Nebula.

The relationship between $I$ and $V$ depends upon the detailed boundary conditions but, if the electromagnetic field remains force-free, then $V\sim IZ_{0}/2\pi$, where $Z_{0}\equiv 377\,{\rm Ohm}$ is the impedance of free space. To order of magnitude, a potential difference of, say, $10\,{\rm EV}$ is then associated with an electromagnetic power $\sim 10^{36}\,{\rm W}$. Whether or not this potential difference is used for efficient UHECR particle acceleration depends upon circumstance but, if it is, then the associated power is inescapable.

Now consider the application of these general ideas to neutron stars. Spinning magnetized neutron stars are, manifestly, not axisymmetric but quite similar considerations apply. If the magnetic moment is that of a regular pulsar, then the spin period must be close to the maximal allowed value $\sim 1.5\,{\rm ms}$ and the source lifetime is only years (Fang:2012rx, ). Magnetars, with surface magnetic field $\sim 100\,{\rm GT}$ are more promising (Arons:2002yj, ). A spin period $\sim 30\,{\rm ms}$ suffices but the source lifetime would be less than a day. It is hard to see how UHECR could escape the environs of a neutron star so soon after its birth given the high density of matter and radiation which should cause catastrophic losses.

However, there is another possibility and this is that the power derives from the magnetic, not the rotational energy. This can be as much as $10^{41}\,{\rm J}$ and it could be released intermittently, following magnetic flares, where $\Omega$ would be replaced by the reciprocal of the time it would take light to cross a fraction of the magnetic surface, a few $\mu s$. Furthermore, the magnetic energy stored below the surface might be hundreds of times larger than that within the magnetosphere and accessible for much longer. This is important because the magnetar birth rate should be less than $\sim 10^{-6}\,{\rm Mpc}^{-3}\,{\rm yr}^{-1}$. requiring an energy per magnetar of $>10^{43}\,{\rm J}$ to account for the UHECR luminosity density.

Despite this energetic discouragement, it is instructive to consider the general principles of particle acceleration exhibited by an axisymmetric pulsar magnetosphere. In this case, the magnetic field is stationary and the EMF formally vanishes. The magnetic surfaces are electrostatically equipotential. However, a current must flow which may require the continuous formation of electron-positron pairs in the magnetosphere. It has long been supposed that this can happen at a “gap”, a possibly transient potential drop along the magnetic field where ${\bf E}\cdot{\bf B}\neq 0$ and single charged particles can be accelerated to sufficient energy to produce $\gamma$-rays which can form pairs and breakdown the vacuum. There has been much modeling of gaps and generally it has been found that the associated potential drops across a gap are order of magnitude smaller than what would be needed to create UHECR directly at gaps. They could, of course be associated with the emission of much less energetic $\gamma$-rays and cosmic rays.

If there is minimal potential drop along the flux surfaces, then the power must associated with electrical current crossing between the flux surfaces. There are many different ways that this has been modeled. For example, there could be enough plasma present that there is a Lorentz force density, ${\bf j}\times{\bf B}$, that gradually converts the electromagnetic energy into kinetic energy flux with minimal dissipation. This outflow might then be an efficient particle accelerator in a jet or wind (see below). Alternatively, the current may complete dissipatively as particle acceleration. At one extreme, all of the power may heat the thermal plasma and cause it to radiate. An intermediate possibility is that much of this power creates a suprathermal distribution of intermediate energy cosmic rays.

A third possibility is that it is the highest energy particles that carry most of the electrical current across magnetic surfaces. The manner in which this may happen is interesting. A modest energy particle will gyrate around and move along the magnetic field in the local frame in which the electric field vanishes. However, for larger particle energies the gyro radius, $r_{L}$, increases and the collisionless particle motions also include drift velocities of their guiding centers with speeds $\sim cr_{L}/L$, where $L$ is the scale length. It is possible that despite the rarity of the highest energy particles, their greater mobility can allow them to carry most of the electrical current and, consequently, be accelerated to energies as large as those for which $r_{L}$ approaches $L$, specifically, those with energy $\sim ZecBL$. This basic mechanism has been invoked to account for intermediate energy particle acceleration at a pulsar wind termination shock. For this to be possible, in principle, requires that the electrical current flows along and not across magnetic surfaces all the way from the pulsar magnetosphere to the surrounding nebula.

Many of these basic physics ideas are translatable to black hole magnetospheres with some important differences. We suppose that magnetic flux is trapped by the inertia and pressure of orbiting gas and some of this flux passes in and out of the event horizon. The black hole is not a perfect conductor. In fact, a formal surface conductivity can be associated with the horizon associating a resistance of $\sim Z_{0}/2\pi$ with the black hole. This has two consequences. The magnetic field lines rotate with an angular velocity roughly half that of the black hole. There is as much power dissipated behind the horizon as there is extracted outside the black hole.

For the case of massive black holes, the power associated with UHECR acceleration is then, at least, comparable with that radiated by a moderate AGN such as a FRI radio source or a Seyfert galaxy. There is the same need to supply electrical charge continuously and gaps have been invoked as sources of $\gamma$-rays (e.g. (Levinson:2010fc, )). As with pulsars, the potential differences that need to be applied to keep the current flowing are generally tiny compared with those needed to account for UHECR. To drive home this point, if the current carrying pair plasma had the minimum density and energy $\sim 1\,{\rm MeV}$ required, then its energy density would be only a fraction $\sim$ the ratio of this energy to the total potential difference — $\sim 1\,{\rm MeV}$ to $\sim 10\,{\rm EV}$ or $10^{-13}$. (This is also the ratio of the minimum to the maximum gyro radii.) In fact gaps will not even be required if either the magnetic field lines are able to undergo interchange instability to keep the magnetosphere supplied with electron-ion plasma or externally produced $\gamma$-rays can produce pairs directly. For all of these reasons, gaps are not seen as UHECR sources.

As with pulsars, it also possible to imagine particle acceleration occurring when current crosses the toroidal magnetic field lines associated with the electromagnetic outflow produced by a magnetized, spinning black hole. This takes place in the “Espresso” model of UHECR acceleration (Caprioli_2015, ) where lower energy cosmic rays are injected into the outflow and emerge with rigidity comparable with the total potential difference available — an expression of the famous Hillas criterion (1984ARA&A..22..425H, ) (see section 3.2.1 below). Assuming $\Gamma\sim 30$ (as in powerful blazars), a one-shot boost of a factor of $\sim\Gamma^{2}$ in energy, can transform the highest-energy galactic CRs at 100 PeV in the highest-energy UHECRs at 100 EeV. If this mechamism can account for the whole UHECR data is still an open question.

There is a second, quite general constraint on locating UHECR acceleration close to an AGN black hole. A $\sim 100\,\mu$ far infrared photon will be able to create pions in the Coulomb field of a $\sim 10\,{\rm EeV}$ cosmic ray. The UHECR acceleration site must be sufficiently far away from the AGN that the photon density does not prevent energy loss or photodisintegration. Similar remarks apply to the pair production threshold which involves $\sim 30\,{\rm GHz}$ radio photons. As discussed below, even more stringent constraints must be satisfied to avoid photodisintegration of the heaviest nuclei. The end result is that the only way that processes directly connected to an AGN magnetosphere could account for UHECR is if they involve low luminosity AGN and operate at some remove from the actual black hole magnetosphere.

Similar arguments can be used to rule out UHECR production around or directly associated with the magnetospheres of stellar mass black holes either close to formation in GRB or later on in binary X-ray sources. We now turn to the viability of UHECR acceleration further away from massive or stellar black holes, within relativistic jets.

Let us recall the jet power requirement for UHECR production. In the following we consider a relativistic plasma jet with half opening angle is $\theta$, $r$ is the jet height above the launching region. We use masses in units of solar mass $\bar{m}={M}/{M_{\odot}}$, luminosities in units of $L_{Edd}=1.26\,10^{38}\bar{m}$ erg/s and radii in units of gravitational radius $\bar{r}={r}/{r_{g}}$ with $r_{g}=3\cdot 10^{5}\bar{m}$ cm. Quantities in the jet comoving frame are denoted with a prime symbol. The jet expansion time in the comoving frame is $t^{\prime}_{\rm exp}=r/(\Gamma\beta c)$, where $\Gamma$ is the bulk Lorentz factor of the jet $\Gamma=(1-\beta^{2})^{-1/2}$. The condition for jet stability implies that the comoving signal crossing timescale is smaller than $t_{\rm exp}$, which gives the condition $\Gamma\theta\lesssim\beta_{s}/\beta$ where $\beta_{s}c$ is the Alfvén velocity. In relativistic magnetized jets $\beta\approx\beta_{s}\approx 1$. The acceleration region cannot extend further than the transverse distance $r\theta$. The escape time sideways is thus $t^{\prime}_{\rm esc}=(r\theta)/(\Gamma c)$.

The jet luminosity $L_{j}=\pi(\theta r)^{2}\beta\Gamma^{2}cu^{\prime}_{j}$ where $u^{\prime}_{j}$ is the comoving energy density. The magnetic field of a relativistic jet with magnetic luminosity $L_{B}\equiv\xi_{B}L_{j}$, bulk Lorentz factor $\Gamma$ and half opening angle $\theta$ is

	$\displaystyle B^{\prime}$	$\displaystyle=$	$\displaystyle\frac{2\left(\xi_{B}L_{j}/c\right)^{1/2}}{\theta\Gamma r}$		(1)
		$\displaystyle\approx$	$\displaystyle 4.3\cdot 10^{8}(\xi_{B}\bar{L}_{j}/\bar{m})^{1/2}\bar{r}^{-1}(\theta\Gamma)^{-1}{\rm G}\,.$

The magnetic field $B^{\prime}$ is related to the total energy density $u^{\prime}_{j}$ by: $B^{\prime}=\sqrt{4\pi\xi_{B}u^{\prime}_{j}}$.

The acceleration time $t^{\prime}_{acc}=t^{\prime}_{L}/\xi_{acc}=E^{\prime}/(ZeB^{\prime}c\xi_{acc})$ with $\xi_{acc}\lesssim 1$ (in numerical experiences $\xi_{acc}\lesssim 0.1$ making the usually assumed Bohm diffusion overoptimistic). The condition $t^{\prime}_{acc}=t^{\prime}_{esc}$ (following (1984ARA&A..22..425H, )) gives:

	$\displaystyle E^{\prime}_{\rm max}$	$\displaystyle=$	$\displaystyle\frac{\xi_{acc}r\theta ZeB^{\prime}}{\Gamma}\approx 300\,\xi_{acc}Z[B^{\prime}_{\rm G}r_{\rm cm}](\Gamma\theta)\Gamma^{-2}\,{\rm eV}$		(2)
		$\displaystyle\approx$	$\displaystyle 3.\cdot 10^{16}\xi_{acc}Z\left(\xi_{B}\bar{L}_{j}\bar{m}\right)^{1/2}\Gamma^{-2}{\rm eV}\,,$

which gives the condition on the observed jet luminosity:

$$L_{j}>8.4\cdot 10^{44}\xi_{B}^{-1}\Gamma^{2}\xi_{acc}^{-2}\left(\frac{E_{\rm max}}{Z\,10^{20}{\rm eV}}\right)^{2}{\rm erg\,s^{-1}}\ ,$$

(3)

where $E_{\rm max}/Z$ is the maximum rigidity reachable in the observer’s frame.

It is important to note that here only the mean intensity of the magnetic field $B^{\prime}$ is considered. However in the usual model of diffusive shock acceleration (DSA), the largest turbulence scale of the magnetic field plays an important role, as the relevant escape time is the escape time upstream of the shock front. As the shocks microphysics is poorly understood (especially in the mildly-relativistic regime), the usual prescription for the particles confinement is that a cosmic ray remains accelerated at the shock if the Larmor radius of the particle is smaller than the largest turbulence scale of the magnetic field in the shocked region, i.e. $r_{L}\lesssim\lambda_{\rm max}$. It is important to keep in mind that $\lambda_{\rm max}$ can be a small fraction of the size of the system. Therefore, in the case of DSA, $t^{\prime}_{\rm esc}=\xi_{esc}(r\theta)/(\Gamma c)$ with $\xi_{esc}\lesssim 1$.
Also, the above estimate of the maximum rigidity reachable is based on a geometrical argument. However there are several energy losses that will limit the acceleration process. The maximum comoving energy $E^{\prime}_{\rm max}$ of the accelerated particles is estimated by equating their acceleration time $t^{\prime}_{acc}$ with the relevant loss time $t^{\prime}_{loss}$, and the related "opacity" is defined as $\tau_{loss}\equiv{t^{\prime}_{esc}}/{t^{\prime}_{loss}}$.

The synchrotron opacity is given by

$$\tau_{syn}=\frac{(r\theta)Z^{4}\sigma_{T}B^{2}E}{6\pi m_{e}^{2}c^{4}\Gamma}\left(\frac{m_{e}}{m_{A}}\right)^{4}\,.$$

(4)

Concerning photodissociation, the lowest energy and highest cross-section process is the giant dipole resonance. The energy loss rate due to GDR interaction with the photons of energy $\epsilon^{\prime}$ is $(t^{\prime}_{\rm GDR})^{-1}\simeq\sigma_{GDR}n^{\prime}(\epsilon^{\prime})c{\Delta\epsilon_{\rm GDR}}/{\epsilon_{\rm GDR}}\,\,{\rm s^{-1}}$ with $\sigma_{GDR}=1.45\,10^{-27}A$ cm², $\Delta\epsilon_{\rm GDR}=8$ MeV and $\epsilon_{\rm GDR}\approx 42.65A^{-0.21}$ MeV. The photon energy $\epsilon^{\prime}$ is given by the threshold condition for GDR interaction: $E^{\prime}_{\rm N}(A)\epsilon^{\prime}\approx 2.5\,10^{17}(A/56)$ eV². The GDR opacity is $\tau_{GDR}\equiv{t^{\prime}_{esc}}/{t^{\prime}_{GDR}}$. With an opacity close to 1, we are in the optimum working condition for the energy loss process considered. Only when the shock becomes transparent the maximum rigidity is limited by the escape from the magnetized region. Particle escape acts as a high pass filter: only particles close to the maximum rigidity $\simeq E^{\prime}_{\rm max}/Z$ can escape from the magnetized region upstream of the shock wave. We expect that the spectrum for the escaping particles is thus very hard, almost a delta function around $\simeq E^{\prime}_{\rm max}/Z$.

Gamma-ray bursts (GRBs) are associated with relativistic jets from newly born stellar-mass black holes. There are three different types of GRBs, which are associated with different progenitors. Long GRBs (lGRBs) are associated with the core-collapse of massive stars (more precisely, to SN Ib/Ic). They are the most powerful type of GRBs, with a jet luminosity able to reach $\sim 10^{52}$ erg/s. They are also the rarest, with a rate of $\sim 1$ Gpc^-3yr^-1 and a $\gamma$-ray luminosity density $\sim 6\,10^{42}$ erg Mpc^-3yr^-1. Short GRBs (sGRBs) are the aftermath of binary neutron star and black hole neutron star mergers and have jet luminosities of typically $\lesssim 10^{49}$ erg s^-1. Finally, there is another class of GRBs, the low luminosity GRBs (LLGRBs) which have lower luminosities (typically $L_{\rm jet}\lesssim 10^{50}$ erg s^-1, $\sim 100$ times the rate of lGRBs (2007ApJ…662.1111L, ).

In the case of GRBs jets, there are two types of shocks that have been widely considered for diffusive shock acceleration: external shocks and internal shocks (1995PhRvL..75..386W, ; 1995ApJ…453..883V, ). External shocks are caused by the jet deceleration in the ambient medium (these are responsible for the GRB afterglow) and are superluminal. Many authors pointed out the inefficiency of Fermi acceleration at ultra-relativistic shocks, unless a strong small-scale turbulence is present in the downstream medium (2006ApJ…641..984N, ; 2006ApJ…645L.129L, ). Acceleration in ultra-relativistic external shocks can produce cosmic rays with energies above a few $Z\times 10^{15}$ eV (PeVatrons) see e.g., (1999MNRAS.305L…6G, ; 2021ApJ…915L…4G, )- but not the highest energy cosmic rays.

Internal shocks are expected to form inside the jet once the fast layers of the plasma jet catch up with the slower parts. In that case, the shock in the comoving frame is mildly relativistic (with Lorentz factor of $\sim$1.1 to $\sim$2) so efficient diffusive shock acceleration can take place. It has been shown that the spectrum of UHECR emitted by GRB internal shocks can reach maximum rigidities of 8 EV (iron at $\sim$100 EeV) for the highest luminosities ($L_{\rm jet}\gtrsim 10^{52}$ erg/s), making diffusive shocks acceleration at internal shocks a promising candidate for UHECR production (10.1093/mnras/stv893, ). However, the caveat for the long GRB model is that they fall short to account for the required UHECR luminosity density (typically by a factor $\sim 100$) unless they have a poor $\gamma$-ray efficiency (2010ApJ…722..543E, ; 10.1093/mnras/stv893, ); or, it has been proposed that one single Galactic event can be responsible for most of the UHECR flux if the local extragalactic magnetic field is large and can trap efficiently $\lesssim$8 EV rigidity particles (Eichler:2016mut, ).

The jets associated with massive black holes in active galactic nuclei, where the outflow is magnetically independent of a central black hole/accretion disk source, have also been proposed as UHECR accelerators. (By magnetically independent, we mean that the large-scale electrical current that flows along the inner jet has mostly left the jet, to be replaced by small scale current supporting small scale magnetic field.) Here two main mechanisms have been proposed, relativistic shock fonts and magnetic reconnection in a relativistic boundary layer. As just pointed out, relativistic shocks are qualitatively different from their nonrelativistic counterparts and they may be much less efficient in accelerating the highest energy particles. Partly for this reason, it has been suggested that it is the nonrelativistic backflow after a relativistic jet passes through a relativistic termination shock that is where more efficient UHECR acceleration can occur (10.1093/mnras/sty2936, ). As the number of suitable sources within the UHECR horizon is small, there is the future possibility of making tentative associations with the highest energy particles may be possible (Bell:2021pkk, ).

There is good observational evidence that relativistic jets decelerate, dissipate, accelerate and radiate through their surfaces, especially close to their massive black hole sources. Presumably gas is entrained into the electromagnetic jet from a much more slowly moving external medium. These can be characterized as boundary layers from a gas dynamical perspective and as a cylindrical, relativistic current sheet from an electromagnetic perspective. It is presumably the second description that is most promising as an accelerator. Either way, this is a site where magnetic reconnection will take place. Non-relativistic reconnection is a process which has been extensively studies and there is much space physics data to compare with simulations. Overall it is not very efficent in creating very high energy particles. However, relativistic acceleration has recently been studied quite extensively and it involves novel features, especially the launching of high speed plasmoids (Sironi:2020mzd, ). There does not seem to be any natural way to accelerate particles to the very highest energy but this is certainly worth further consideration.

There has been much recent observational attention paid to Tidal Disruption Events (TDE), where individual stars orbit massive ($M\lesssim 10^{8}\,{\rm M}_{\odot}$) black holes in otherwise dormant galaxies passing by so close that they are they are ripped apart by tidal forces with the debris falling back onto the black hole on timescales of months to years and radiating in the optical and X-ray bands especially. A minority of TDEs are found to be jetted and there is evidence that the jets have Lorentz factors $\gtrsim 10$, like blazars (10.1093/mnras/stad344, ). They have also been proposed as sources of UHECR (Farrar:2014yla, ).

The estimated rate of jetted TDEs beamed towards us is $2\times 10^{-11}\,{\rm Mpc}^{-3}\,{\rm yr}^{-1}$ (Andreoni:2022afu, ). A conventional estimate of the total rate of jetted TDE is perhaps a hundred times larger. On this basis, we can estimate the energy required per jetted TDE to account for ${\cal L}_{\rm UHECR}$ (Sec. 2) as $\sim 0.03\,{\rm M}_{\odot}c^{2}$. This is not impossible but it seems implausible given the very much lower radiative efficiencies typically inferred.

An alternative class of UHECR acceleration models invokes Diffusive Shock Acceleration (DSA, see (BE87, ) for a review) at giant intergalactic shock fronts, most notably those surrounding rich clusters of galaxies but also at shock fronts surrounding the filaments which connect these clusters. There are many possibilities. Here, we will strict attention to the more limited but more prescriptive and, consequently, more refutable idea that the entire cosmic ray spectrum observed at Earth is the result of successive DSA from distinct source classes.

Specifically, we consider a hierarchical model in which the observed $\sim 1\,{\rm GeV}-\sim 3\,{\rm PeV}$ cosmic rays originate in nearby supernova remnants. The pressure from these cosmic rays drives a magnetocentrifugal Galactic wind which passes through a Galactic Wind Termination Shock (GWTS) where re-acceleration takes place. Some of the highest energy particles from the GWTS escape upstream and contributes to the shin part of the spectrum; the remainder are transmitted into the circumgalactic and intergalactic media. GWTS acceleration is even more powerful around galaxies associated with AGN and starbursts.

These intermediate energy cosmic rays permeate the intergalactic medium and are input to DSA at the large intergalactic shocks. This “holisticââ interpretation features propagation as well as acceleration. Of course, there are many alternative acceleration sites that could dominate different parts of the spectrum in a hierarchical scheme, as discussed above, or, alternatively, the intergalactic acceleration might operate mainly on mildly relativistic particles injected directly at the shock front.

DSA is a specific mechanism for converting a significant fraction of the kinetic energy flux, measured in the frame of the shock, of the gas upstream of a strong shock front into high energy particles that are transmitted downstream. It relies upon scattering by hydromagnetic disturbances/waves. Individual particles of rigidity (momentum per unit charge) $R$, with mean free paths $\ell(R)$, will diffuse against the flow with an upstream scale height $L\sim\ell c/u$, where $u$ is the speed of the shock relative to the upstream gas. In the simplest, test particle version of this mechanism at a planar shock front, relativistic particles will increase their rigidity by a fractional amount $\Delta R/R\sim 4(1-1/r)u/3c$, where $r$ is the shock compression ratio, every time they make a double crossing of the shock front. On kinematic grounds, the probability that an individual cosmic ray not return upstream is $4u/rc$. The probability that a relativistic particle incident upon the shock front with initial rigidity $R$ be transmitted downstream with rigidity greater than $xR$ is

$$P(>x)=x^{\frac{3}{r-1}};\quad x>\geq 1.$$

(5)

In this limit the shock behaves as a linear system, convolving the input upstream momentum space particle distribution function with a power law, Greensâ function to give the transmitted downstream distribution function. The slope of the power law does not depend upon the details of the scattering and no particles escape upstream. The maximum rigidity, $R_{\rm max}$ that can be accelerated by a given shock is usually determined by the condition that $L$ must be less than the size of the shock, typically its radius of curvature. Note that this condition is, at least, $c/u$ times more stringent that the Hillas criterion. Alternatively, it takes $\sim(c/u)^{2}$ scatterings, not one orbit, to double a particle’s energy. DSA is inherently slower than the more violent mechanisms discussed above. The question is “Can it still account for the highest energy particles we actually observe?”.

Realistically, the test particle approximation fails in many ways. The acceleration is influenced by the cosmic-ray pressure gradient decelerating the gas ahead of the shock and changing the kinematics. In addition the cosmic rays can account for a significant fraction of the total energy of the flow and, on account of their different effective specific heat ratio, change the determination of $r$ through the conservation laws. The accelerated cosmic rays also create and sustain the wave turbulence that scatters them (10.1093/mnras/172.3.557, ). The most direct, and, arguably, the most important, form of scattering involves waves with wavelengths resonant with the particle gyro radii $r_{\rm L}(R)$. If the rms magnetic field amplitude of these waves is $\delta B(r_{\rm L})=\delta B(R)$, then $\ell\sim(B/\delta B)^{2}r_{\rm L}$, assuming random phases and Gaussian statistics (both questionable). If we assume that the waves are nonlinear over some range of $R$, then $\ell\sim r_{\rm L}$ and we call this Bohm diffusion. The efficiency of cosmic ray acceleration, $\epsilon_{\rm CR}$, which can be measured by the ratio of the downstream cosmic ray pressure to the upstream momentum flux, is then determined by the processes that determine the rms magnetic field strength.

There are several possible sources of this wave turbulence. If cosmic rays stream with mean speed in excess of the Alfvén speed in a uniform magnetic field, then they will excite the growth of resonant Alfvén waves propagating along the direction of their mean velocity. These waves will then scatter the cosmic rays in pitch angle and reduce their mean motion so that it is not much more than the Alfvén speed. The linear growth rate is roughly $(n_{\rm CR}/n_{i})c/r_{\rm L}$, where $n_{\rm CR}$ measures the number of resonant particles and $n_{i}$ is the background ion density. An alternative possibility is that waves are created hydromagnetically with wavelength much larger than $r_{\rm L}$ which initiates a turbulence spectrum where the waves with shorter wavelengths will have linear amplitudes. A third possibility is that the waves are generated with wavelengths shorter than $r_{\rm L}$ by the motion of electrons carrying the return current that must balance the current associated with the accelerated cosmic rays. Extensive “PIC” and related simulations have vindicated the generation of wave turbulence and elucidated how these and other processes interplay in DSA.

Although we are mostly concerned with questions of UHECR acceleration, it is instructive to contrast intergalactic shocks with precursor acceleration of lower energy particles at smaller shock fronts. It is a reasonable conjecture that analogous physical processes operate at all high Mach number shocks, scaled by $L$ and the momentum flux. While there has long been strong evidence associating supernova remnants with the majority of Galactic cosmic rays up to PeV energy, the argument for the reminder of the spectrum being due to DSA, is no better than circumstantial. What is clear is that if DSA accounts for most of the total spectrum, the acceleration must occur with near maximal efficiency with respect to $R_{\rm max}$ and $\epsilon_{\rm CR}$ in all three sites. For this reason, it is interesting to look at this problem inductively — to ask what would be required from the complex nonlinear plasma processes in the vicinity of the shock to account for the observations— rather than deductively, to attempt to calculate these processes from first principles.

Despite this commonality, there are important distinctions to be made. The first involves the kinematics. Most supernova remnants exhibit a quasi-spherical shock, convex on the upstream side, expanding into a near uniform interstellar medium, though stellar winds and molecular clouds can certainly introduce observable differences in many examples. This geometry facilitates the escape of the highest energy particle upstream as opposed to their transmission downstream, as assumed under the test particle approximation. Over half of the very highest energy particles can leave the remnant in this fashion. The original model of DSA essentially ignored these particles, supposing, instead, that the cosmic rays that we observed were transmitted downstream and were then released, following adiabatic decompression at the end of the remnant’s lifetime. By contrast, it could be that most of the observed cosmic rays escape upstream with $R\sim R_{\rm max}$ which decreases as the remnant ages.

To be quantitative, the cosmic ray luminosity, per unit area, of the Galactic disk is generally estimated to be at least ten percent of the local supernova remnant luminosity per unit disk area. ($E_{\rm SNR}\sim 10^{51}\,{\rm erg}$ released every $\sim 30\,{\rm yr}$ throughout the Galaxy.) Now, in order to accelerate protons to $\sim 3\,{\rm PeV}$ energy at a SNR expanding with speed $u_{\rm SNR}$ and radius $R_{\rm SNR}$, requires that the gyro radius be $\lesssim u_{\rm SNR}R_{\rm SNR}/c$ or $B\gtrsim 100(u_{\rm SNR}/10,000\,{\rm km\,s}^{-1})^{-1}(R_{\rm SNR}/1\,{\rm pc})^{-1}\,\mu{\rm G}$, perhaps two orders of magnitude large than the rms Galactic field. Furthermore, the magnetic field strength is strongly bounded above by requiring that its contribution to the momentum flux, $B^{2}/4\pi$ be less than that of the gas, which is $\sim 0.3E_{\rm SNR}R_{\rm SNR}^{-3}$, adopting the Sedov solution, or $B\lesssim 10(R_{\rm SNR}/1{\rm pc})^{-3/2}\,{\rm mG}$. This interval allows efficient cosmic ray acceleration by young SNR up to the knee in the spectrum.

The next level in this proposed hierarchy involves galactic winds. A large motivation for proposing a wind in our Galaxy comes from cosmic-ray observations. The measured ratio of light to medium nuclides at $\sim\,{\rm GeV}$ energy allows us to infer that these cosmic rays escape the disk after traversing a grammage $\sim 6\,{\rm g\,cm}^{-2}$. This implies a residence time $\sim 30-100\,{\rm Myr}$, much shorter than a disk rotation period. Cosmic rays leave the disk continuously and it is reasonable to suppose that they carry gas with them. The wind provides a chimney or exhaust for the cosmic rays (as well as the hot gas which is unable to cool radiatively) created by the SNR. The existence of an outflowing wind, occupying most of the volume of the Galactic halo, is not inconsistent with the presence of denser gas clouds falling inward, as observed.

Various models of this outflow have been entertained. One recent model posits the presence of a vertical magnetic field at the disk and that $\sim\,{\rm GeV}$ cosmic rays stream along it just faster than the Alfvén speed so that Alfvén waves are excited and couple the cosmic rays to the gas, allowing it to achieve the escape speed from the solar radius. The cosmic rays alone drive the wind (Mukhopadhyay:2023kel, ).

A more powerful wind is possible if the vertical magnetic field corotates with the disk. In this case the cosmic rays just elevate the gas to modest altitude $\sim$ a few kpc, at which point the field develops a large radial component and the gas is flung out magnetocentrifugally and it is the orbiting interstellar medium that provides most of the power for the wind. Under this scenario, the interstellar medium should be seen as an open system losing a significant fraction of its mass, and angular momentum as well as energy over a Galactic lifetime. This possibility is not generally entertained in models of the interstellar medium.

The outflow will pass through three magnetohydrodynamic (Alfvénic and magnetosonic) critical points before crossing a termination shock, typically at a radius $R_{\rm GWTS}\sim 200\,{\rm kpc}$, with speed $u_{\rm GWTS}\sim 1000\,{\rm km\,s}^{-1}$, before joining the circumgalactic medium. The flow expands with density $\rho$ decreasing by roughly three orders of magnitude. The cosmic ray rigidities decrease $\propto\rho^{1/3}$, roughly one order of magnitude. The poloidal component of the magnetic field decreases from $\sim 2\,\mu{\rm G}$ at the disk to $B_{p}\sim 3\,{\rm nG}$ at the shock, while the toroidal field is roughly ten times larger, $B_{\phi}\sim 30\,{\rm nG}$.

The termination shock can reaccelerate $\sim{\rm GeV-PeV}$ cosmic rays to fill out the $\sim\,{\rm PeV-EeV}$ shin in the spectrum. The shock surface is likely stationary and concave on the upstream side with radius. This changes DSA in significant and observable ways. For modest cosmic-ray rigidity, say below the reduced peak of the cosmic-ray spectrum $R\lesssim 3\,{\rm PV}$, the shock will be effectively planar and a reaccelerated spectrum will be transmitted downstream. However, when $\ell\gtrsim u_{\rm GWTS}R_{\rm GWTS}/c$, the high energy cosmic rays are able to cross the Galactic halo and continue their acceleration at a quite different part of the shock. This acceleration can continue up to a rigidity where $\ell\sim R_{\rm GWTS}$, perhaps $\sim 100\,{\rm PV}$. When, as is usually the case, cosmic ray energy rather than rigidity spectra are displayed, the composition of the shin particles becomes progressively heavier, as is observed.

Again, to be quantitative, the wind may account for as much as $\sim 1\,{\rm M}_{\odot}{\rm yr}^{-1}$mass loss from the galaxy and a corresponding power $\sim 5\times 10^{41}\,{\rm erg\,s}^{-1}$, more than the cosmic ray power. A turbulent magnetic field strength $B_{\rm RMS}(100\,{\rm PV})\sim 0.3B_{\phi}(R_{\rm GWTS})$ suffices to sustain acceleration within a closed shock surface throughout the shin observed at Earth.

We now turn to the final and most relevant level in this simple hierarchy when the $\sim\,{\rm PeV-EeV}$ cosmic rays accelerated by galaxies, especially those that are more active than our Galaxy, are incident upon intergalactic shock fronts (e.g. (1995ApJ…454…60N, ; 1997MNRAS.286..257K, ), Simeon et al. in preparation) . The most powerful shocks, prominent in cosmological simulations, surround rich clusters of galaxies which show Mach numbers that can exceed $\sim 100$. The assembly of clusters, over cosmological time, is complex but we will idealize the shock as a stationary spherical surface surrounding most of the cluster. These shocks have not been observed directly as yet. Deep searches at low radio frequencies seems well-motivated and would be very prescriptive, if successful. Simulations also show quasi-cylindrical filaments that may connect clusters and which may be surrounded by quasi-cylindrical infall shocks. The Milky Way may lie within such a filament.

The kinematics of a stationary cluster accretion shock, CAS, is different from that at SNR and GWTS. The convergence of gas ahead of the shock results in adiabatic heating of lower energy particles. This acceleration, which supplements DSA, will be less efficient for higher energy particles, which diffuse outward, relative to the gas flow. We continue to view the problem inductively and ask what scattering would have to be present in order for nearby cluster shocks to account for the observed UHECR spectrum. If we take an accretion shock surrounding the Virgo cluster then we observation combined with cosmological simulations suggest a shock radius $R_{\rm CAS}\sim 2\,{\rm Mpc}$ an infall speed and density ahead of the shock $u_{-}\sim 1000\,{\rm km\,s}^{-1}$ and $\rho_{-}\sim 10^{-29}\,{\rm g\,cm}^{-3}$. The shock is strong and the rate at which gas kinetic energy crosses it is $\sim 10^{45}\,{\rm erg\,s}^{-1}$. If we adopt a rich cluster density of $\sim 10^{-5}\,{\rm Mpc}^{-3}$, then the cluster gas kinetic energy luminosity density is $\sim 30{\cal L}_{\rm UHECR}$. There seems to be sufficient power to account for the observed UHECR flux. The larger challenge is to accelerate heavy nuclei, optimally, to $R_{\rm max}\sim 10\,{\rm EV}$.

In this model, we detect the UHECR that escape upstream with a spectrum peaking at a rigidity $\sim 1-8\,{\rm EV}$ which will vary from cluster to cluster. The cosmic rays which are transmitted downstream will never be directly observable, though they might be seen at much lower energy through their $\gamma$-ray emission. In order to accelerate cosmic rays to $R_{\rm max}\sim 10\,{\rm EV}$, we require that $R_{\rm CAS}>(c/3u_{-})\ell(R_{\rm max}$. If we adopt Bohm diffusion, then the rms, resonant field ahead of the shock on these scales is $B_{\rm rms}\sim 1\,\mu{\rm G}$. This contrasts with a general field in the IGM that is generally expected to be $\lesssim 1\,{\rm nG}$. The isotropic magnetic pressure associated with this field is roughly an order of magnitude smaller than the nominal momentum flux carried by the background gas $\sim 10^{-13}\,{\rm dyne\,cm}^{-2}$. Provided that the turbulence does not dissipate and pre-heat the gas, thereby reducing the shock Mach number and the efficacy of DSA, a maximum rigidity $R_{\rm max}\sim 10\,{\rm EV}$ seems attainable. However, if, for example, protons with significantly larger rigidity than this are identified, it seems hard to see how accretion shocks around relatively well-observed and simulated, local clusters of galaxies could account for their acceleration.

Now turn to the input. Of course, we do not know what is the intergalactic spectrum produced by all the varied types of galaxy that surround the clusters. There are many possible contributors that we we may not be observing at Earth, for example starbursts and the winds they drive, pulsar wind nebulae, jets associated with spinning black holes and a network of weaker intergalactic shock fronts. (It is possible to associate the observed light ankle with filament shock acceleration in which case this comprises the spectrum transmitted downstream.) If we just confine our attention to the SNR spectrum produced by our local Galactic disk, then we expect a source spectrum $S(R)\propto R^{-2.2}$, as produced by all interstellar shocks over their lifetimes. This should extend up to a rigidity of order a few PV. At this point the rigidity spectrum should steepen as these additional sources fill out the shin part of the spectrum. (An energy spectrum will exhibit progressively heavier nuclei as observed.) Incident cosmic rays with rigidity well below the UHECR range will be convected adiabatically into the shock and see it as an essentially planar shock accelerator. Provided that the compression ratio of the subshock is neither too large, specifically $2.5\lesssim r\lesssim 3.5$, then it will be the cosmic rays near the knee in the spectrum at few $PV$ that will dominate the UHECR spectrum that is reflected upstream back into the IGM. This will be reflected in the UHECR composition. In this model, the GeV particles that dominate the GCR spectrum and freshly injected particles are minor contributors to the output UHECR spectrum.

How reasonable is it for a CAS to behave in this extreme fashion? And, by extension, how reasonable is it for SNR and GWTS to accelerate cosmic rays to energies $\sim 3$ and $\sim 300\,{\rm PeV}$ respectively? Here, a rather different approach from that followed in most of the impressive simulations that have been completed may be helpful. Instead of starting from a quiescent initial state and monitoring the progressive acceleration of cosmic rays, and their associated local turbulence, to ever greater rigidity, limited by dynamic range and runtime, it might be interesting to start with the near-maximally turbulent state described and see if it can be self-sustaining. The dynamics is dominated by the highest energy particles and the background gas. PeV cosmic rays can be injected at the shock front. A key feature of this approach is that, as the rigidity increases, the particles diffuse further ahead of the shock front until they either escape upstream or are transmitted downstream with probabilities $\sim$ a half, ending their acceleration. It will also be necessary to include the curvature of the shock which facilitates this escape. These particles, with $R\sim R_{\rm max}$, will provide the first encounter that the inflowing intergalactic gas has with the shock. The gas itself is, presumably, still quite weakly magnetised with a pressure much smaller than the (anisotropic) pressure of the UHECR. Indeed, the pressure of the cosmic rays is likely to dominate the thermal pressure of the IGM. There is the possibility of hydromagnetic instability, similar to the firehose and mirror instabilities but also including the mean velocity of the UHECR fluid through the gas. This turbulence should grow as the gas falls further into the cluster with wavelengths longer than the local UHECR gyro radius ($\sim 1-10B_{\mu{\rm G}}^{-1}\,{\rm kpc}$), but shorter than the distance from the shock. The turbulence, which should be quasi-stationary in space, will have decreasing outer scale as the shock is approached, and should evolve to shorter wavelengths, progressively scattering lower rigidity particles through strongly nonlinear, resonant interaction as the shock is approached.

This “bootstrap” mechanism (adopting the physics as opposed to the statistics usage (Blandford:2007zza, )) provides a different approach to exploring the complex physics of DSA. In the context of cluster shocks, it needs to be augmented by the inclusion of photopion and photodisintegration loss which are important because, unlike with the SNS, GWTS cases, the acceleration timescale is not short compared with the propagation timescale. Indeed, intergalactic and Galactic diffusion needs to be included in the problem alongside the background velocity field which transitions from inflow to cosmological outflow about 7 Mpc away from the cluster for Virgo (Shaya:2017wnh, ). The complex thermal and multiphase nature of the IGM, which is likely to be influenced by the hypothesized magnetic turbulence, is also a large part of the problem. Finally, the downstream boundary condition on particle acceleration, which is normally idealized as strong scattering in a uniform flow could turn out to be important.