# Critical self-organization of astrophysical shocks

## 1. Introduction

It is now generally recognized that most of the observed gamma radiation derives in one way or another from accelerated particles. Radio and x-ray spectral components from a variety of astrophysical objects are believed to have a similar origin. High energy neutrinos, whose detection is on the program for the future and existing water/ice detectors, must also be related to the ultra-high energy cosmic rays (UHECR). Their origin is, in turn, a mystery and the huge Auger detector complex is now being built to elucidate it (Blandford (1999); Cronin (1999)).

There has been essential progress in our understanding of how accelerated particles produce the radiation detected. The models concentrate on the synchrotron and inverse Compton emission for electrons and on the $`\gamma `$ ray and neutrino production in $`pp`$ and $`p\gamma `$ reactions for protons. However, the primary spectrum of accelerated particles remains a stumbling block making predictions of otherwise similar models so different (e.g., Mannheim et al.,  (1999); Waxman & Bahcall (1999)).

The “standard” mechanism of particle acceleration, capable of producing nonthermal power-law spectra extending over many decades in energy is the I-order Fermi or diffusive shock acceleration. It was originally suggested to explain the origin of galactic cosmic rays (CRs). For the purposes of the high-energy radiation and UHECR it is usually adopted as an axiom, and mostly only in its simplest, test particle (TP) or linear realization. In particular, it is assumed that any strong nonrelativistic shock routinely produces a $`E^2`$ spectrum of protons and/or electrons. In fact this spectrum arises from a simple formula $`FE^{(r+2)/(r1)}`$ where $`r`$ is the shock compression (for $`r=4`$). However, it is valid only if the shock thickness is much smaller than the particle mean free path. This, in turn, is true only if the energy content of accelerated particles is small compared to the shock energy (inefficient acceleration) so that the shock structure is maintained by the thermal, not by the high energy particles. Otherwise, the accelerated particles create the shock structure on their own and if so, then obviously on a scale that is larger or of the order of their mean free path, thus making the above formula invalid. Therefore, the TP regime requires a very low CR number density (the rate of injection $`\nu `$ into the acceleration process), which appears to be impossible in the parameter range of interest. It has been inferred from observations (e.g., Lee (1982)), simulations (e.g., Bennet & Ellison (1995)), and theory (e.g., Malkov (1998)) that the CR number density $`n_{\mathrm{CR}}`$ at a strong shock must be $`10^3`$ of the background density $`n_1`$ upstream. It is important to emphasize here that when the actual injection rate $`\nu `$ exceeds the critical value (denote it $`\nu _2`$), the test particle ($`E^2`$) solution simply does not exist. Simple measures, such as calculating corrections, are intrinsically inadequate. What happens is that two other solutions with considerably higher efficiencies branch off at a somewhat lower injection rate $`\nu _1<\nu _2`$, one of which disappears again at $`\nu =\nu _2`$, together with the test particle solution.

Thus, it seems to be difficult to put an accelerating shock into a regime in which the CR energy production rate (acceleration efficiency) could be gradually adjusted by changing parameters. It is either too low (TP regime) or it is close to unity. Note, that this situation is quite suggestive of that occurring in phase transitions or bifurcations.

Generally, neither of those extreme regimes provide an adequate description of particle spectra and related emission. Nevertheless, we argue in this Letter that despite this apparent lack of regulation ability, shocks must be still capable of self-regulation and self-organization. The transition region between the two acceleration regimes (critical region) is very narrow in control parameters like $`\nu `$. On the other hand, the self-regulation can work efficiently only when the parameters are within this region. This requirement determines them, and also resolves the question of the mechanism of self-regulation. The above consideration is similar to the concept of self-organized criticality (SOC) (e.g., Bak, Tang & Wiesenfeld (1987), Hwa & Kardar (1992), Diamond & Hahm (1995)).

## 2. Formulation of the problem

We use the diffusion-convection equation (e.g., Drury (1983)) for describing the distribution of high energy particles (CRs). We assume that the gaseous discontinuity (also called the subshock) is located at $`x=0`$ and the shock propagates in the positive $`x`$\- direction. Thus, the flow velocity in the shock frame can be represented as $`V(x)=u(x)`$ where the (positive) flow speed $`u(x)`$ jumps from $`u_2u(0)`$ downstream to $`u_0u(0+)>u_2`$ across the subshock and then gradually increases up to $`u_1u(+\mathrm{})u_0`$. In a steady state the equation reads

$$u\frac{f}{x}+\kappa (p)\frac{^2f}{x^2}=\frac{1}{3}\frac{du}{dx}p\frac{f}{p},$$
(1)   
where $`f(x,p)`$ is the isotropic (in the local fluid frame) part of the particle distribution. This is assumed to vanish far upstream ($`f0,x\mathrm{}`$), while the only bounded solution downstream is obviously $`f(x,p)=f_0(p)f(0,p)`$. The most plausible assumption about the cosmic ray diffusivity $`\kappa (p)`$ is that of the Bohm type, i.e., $`\kappa (p)=Kp^2/\sqrt{1+p^2}`$ (the particle momentum $`p`$ is normalized to $`mc`$). In other words $`\kappa `$ scales as the gyroradius, $`\kappa r_\mathrm{g}(p)`$. The reference diffusivity $`K`$ depends on the $`\delta B/B`$ level of the MHD turbulence that scatters the particles in pitch angle. The minimum value for $`K`$ would be $`Kmc^3/eB`$ if $`\delta BB`$. Note that this plain parameterization of this important quantity is perhaps the most serious incompleteness of the theory which will be discussed later.

To include the backreaction of accelerated particles on the plasma flow three further equations are needed. First, one simply writes the conservation of the momentum flux in the smooth part of the shock transition ($`x>0`$, the so-called CR-precursor)

$$P_\mathrm{c}+\rho u^2=\rho _1u_1^2,x>0$$
(2)   
where $`P_\mathrm{c}`$ is the pressure of the high energy particles

$$P_\mathrm{c}(x)=\frac{4\pi }{3}mc^2_{p_0}^{p_1}\frac{p^4dp}{\sqrt{p^2+1}}f(p,x)$$
(3)   
It is assumed here that there are no particles with momenta $`p>p_1`$ (they leave the shock vicinity because there are no MHD waves with sufficiently long wave length $`\lambda `$, since the cyclotron resonance requires $`p\lambda `$). The momentum region $`0<p<p_0`$ cannot be described by equation (1) and the behavior of $`f(p)`$ at $`pp_0`$ is described by the injection parameters $`p_0`$ and $`f(p_0)`$ (Malkov (1997), \[M97\]). The plasma density $`\rho (x)`$ can be eliminated from equation (2) by using the continuity equation $`\rho u=\rho _1u_1`$. Finally, the subshock strength $`r_\mathrm{s}`$ can be expressed through the Mach number $`M`$ at $`x=\mathrm{}`$ (e.g., Landau & Lifshitz )

$$r_\mathrm{s}\frac{u_0}{u_2}=\frac{\gamma +1}{\gamma 1+2R^{\gamma +1}M^2}$$
(4)   
where the precursor compression $`Ru_1/u_0`$ and $`\gamma `$ is the adiabatic index of the plasma.

The system of equations (1,2,4) describes in a self-consistent manner the particle spectrum and the flow structure. An efficient way to solve it is to reduce this system to one integral equation (M97). A key dependent variable is an integral transform of the flow profile $`u(x)`$ with a kernel suggested by an asymptotic solution of the system (1)-(2) which has the form

$$f(x,p)=f_0(p)\mathrm{exp}\left[\frac{q}{3\kappa }\mathrm{\Psi }\right]$$

where

$$\mathrm{\Psi }=_0^xu(x^{})𝑑x^{}$$

and the spectral index downstream $`q(p)=d\mathrm{ln}f_0/d\mathrm{ln}p`$. The integral transform is as follows

$$U(p)=\frac{1}{u_1}_0^{\mathrm{}}\mathrm{exp}\left[\frac{q(p)}{3\kappa (p)}\mathrm{\Psi }\right]𝑑u(\mathrm{\Psi })$$
(5)   
and it is related to $`q(p)`$ through the following formula

$$q(p)=\frac{d\mathrm{ln}U}{d\mathrm{ln}p}+\frac{3}{r_\mathrm{s}RU(p)}+3$$
(6)   
Thus, once $`U(p)`$ is found both the flow profile and the particle distribution can be determined by inverting transform (5) and integrating equation (6). Now, using the linearity of equation (2) ($`\rho u=const`$), we derive the integral equation for $`U`$ by applying the transformation (5) to the $`x`$ derivative of equation (2) (M97). The result reads

$`U(t)`$ $`=`$ $`{\displaystyle \frac{r_\mathrm{s}1}{Rr_\mathrm{s}}}+{\displaystyle \frac{\nu }{Kp_0}}{\displaystyle _{t_0}^{t_1}}𝑑t^{}\left[{\displaystyle \frac{1}{\kappa (t^{})}}+{\displaystyle \frac{q(t^{})}{\kappa (t)q(t)}}\right]^1`$ (7)   
$`\times `$ $`{\displaystyle \frac{U(t_0)}{U(t^{})}}\mathrm{exp}\left[{\displaystyle \frac{3}{Rr_\mathrm{s}}}{\displaystyle _{t_0}^t^{}}{\displaystyle \frac{dt^{\prime \prime }}{U(t^{\prime \prime })}}\right]`$   
where $`t=\mathrm{ln}p`$, $`t_{0,1}=\mathrm{ln}p_{0,1}`$. Here the injection parameter

$$\nu =\frac{4\pi }{3}\frac{mc^2}{\rho _1u_1^2}p_0^4f_0(p_0)$$
(8)   
is related to $`R`$ by means of the following equation

$`\nu `$ $`=`$ $`Kp_0\left(1R^1\right)`$ (9)   
$`\times `$ $`\left\{{\displaystyle _{t_0}^{t_1}}\kappa (t)𝑑t{\displaystyle \frac{U(t_0)}{U(t)}}\mathrm{exp}\left[{\displaystyle \frac{3}{Rr_\mathrm{s}}}{\displaystyle _{t_0}^t}{\displaystyle \frac{dt^{}}{U(t^{})}}\right]\right\}^1`$   
The equations (4,7,9) form a closed system that can be easily solved numerically. We analyze the results in the next section.

## 3. Mechanisms of critical self-organization

The critical nature of this acceleration process is best seen in variables $`R,\nu `$. The quantity $`R1`$ is a measure of shock modification produced by CRs, in fact $`(R1)/R=P_\mathrm{c}(0)/\rho _1u_1^2`$ (eq.) and may be regarded as an order parameter. The injection rate $`\nu `$ characterizes the CR density at the shock front and can be tentatively treated as a control parameter. It is convenient to plot the function $`\nu (R)`$ instead of $`R(\nu )`$ (using equation ), since $`R(\nu )`$ is not always a single-valued function, Fig. Critical self-organization of astrophysical shocks.

The injection rate $`\nu `$ at the subshock should be calculated given $`r_\mathrm{s}(R)`$ (M97) with the self-consistent determination of the flow compression $`R`$ on the basis of the $`R(\nu )`$ dependence obtained. However, in view of its critical character, this solution can be physically meaningful only in regimes far from criticality, i.e., when $`R1`$ (test particle regime) or $`R1`$ (efficient acceleration). But, it is difficult to see how this system could stably evolve remaining in one of these two regimes. Indeed, if $`\nu `$ is subcritical it will inevitably become supercritical when $`p_1`$ is sufficiently high. Once it happened, however, the strong subshock reduction (equation ) will reduce $`\nu `$ and drive the system back to the critical regime, Fig. Critical self-organization of astrophysical shocks.

The maximum momentum $`p_1`$ is subject to self-regulation as well. Indeed, when $`R1`$, the generation and propagation of Alfven waves is characterized by strong inclination of the characteristics of wave transport equation towards larger wavenumbers $`k`$ on the $`kx`$ plane due to wave compression. Thus, considering particles with $`pp_1`$ inside the precursor, one sees that they are in resonance with waves that must have been excited by particles with $`p>p_1`$ further upstream but, there are no particles with $`p>p_1`$. Therefore, the required waves can be excited only locally by the same particles with $`pp_1`$ which substantially diminishes the amplitude of waves that are in resonance with particles from the interval $`p_1/R<p<p_1`$. (The left inequality arises from the resonance condition $`kcpeB/mc`$ and the frequency conservation along the characteristics $`ku(x)=const`$). This will worsen the confinement of these particles to the shock front. The quantitative study of this process is the subject of current research. What can be inferred from Fig. Critical self-organization of astrophysical shocks now, is that the decrease of $`p_1`$ straightens and rise the curve $`\nu (R)`$, so that it returns to the monotonic behaviour. However once the actual injection becomes subcritical (and thus $`R1`$) then $`p_1`$ will grow again restoring the two extrema on the curve $`\nu (R)`$.

The above dilemma is quite typical for dynamical systems that are close to criticality or marginal stability. A natural way to resolve it consists in collapsing the extrema into an inflection point so that a self-organized critical (SOC) acceleration regime is established being determined by the conditions $`\nu ^{}(R^{})=\nu ^{\prime \prime }(R^{})=0`$. These conditions not only determine unique critical values $`R^{}`$ and $`\nu ^{}\nu (R^{})`$ but also yield the maximum momentum $`p_1`$ as a function of $`M`$, which is shown in Fig. Critical self-organization of astrophysical shocks. A few particle spectra that develop in the SOC states for different Mach numbers are shown in Fig. Critical self-organization of astrophysical shocks, along with the asymptotic $`M=\mathrm{}`$ non SOC spectrum. The latter can be calculated in a closed form (M97). Note that the hardening of the spectra in about the last decade below the cut-off is entirely due to the abrupt cut-off itself.

## 4. Discussion

The detailed microphysics behind the SOC is extremely complex and must include the self-consistent turbulence evolution and particle acceleration with their strong backreaction on the shock structure. We have simplified it and argued that the most important dynamical components, the bulk plasma flow and the high-energy particles, must be in a balance that constitutes a certain equipartition of the shock energy between the two. This was done by identifying the factors that prevent either of them from prevailing alone.

The above situation is similar to that in e.g., a simple sandpile paradigm of the SOC. It is impossible (in fact unnecessary) to describe the individual grain dynamics, but it is clear that when the critical macroscopic characteristic of the system (the slope of the sandpile) becomes too steep due to the action of external factors, like tilting of the entire system or addition of sand at the top, the sandpile relaxes bringing the slope to its critical magnitude.

## 5. Possible feedback from observations

Perhaps the most significant observational aspect is the particle spectrum. Although the conversion of detected radiation spectra into the primary particle spectra is ambiguous, in some cases it may be compared with the theory. The most striking prediction is that in shocks with very high maximum particle energy, the spectra must be harder than $`q=4`$ (or 2 in the normalization $`f(E)dE`$). This is because of a very low injection requirement for efficient acceleration in such shocks (Fig.Critical self-organization of astrophysical shocks). If we (conservatively) set $`n_{\mathrm{CR}}/n_110^3`$, then $`\nu 10^3cp_0/mu_1^2`$ which may easily exceed $`\nu _2`$ (local maximum on the $`\nu (R)`$ curve) already for $`p_110^410^5`$, putting the acceleration into a strongly nonlinear regime. This should have important observational consequences.

First, the particle energy is concentrated at the upper cut-off instead of being evenly distributed over the logarithmic energy bands as in the test particle $`E^2`$ solution. This makes the upper bounds on CR generated neutrino fluxes (see e.g., Bahcall & Waxman (1999) and Mannheim et al.,  (1999) and references therein) rather ambiguous. Indeed, the UHECR spectrum is normalized to the observed one at $`E_{\mathrm{norm}}10^{19}`$ eV while being obscured by the galactic background at the energies $`E10^{18}`$ eV. Thus, CR spectra harder than $`E^2`$ imply that the upper limit on the neutrino fluxes e.g., derived by Waxman & Bahcall (1999) should be even lower than the $`E^2`$ spectrum implies for the energies $`E_\nu <510^{17}`$ eV ($`E_\nu /E_{\mathrm{CR}}0.03`$ is a typical energy relation). According to the same logic, it should be increased for higher energies. Note that the upper cut-offs in individual shocks contributing to the UHECR must be still much higher than $`E_{\mathrm{norm}}`$ to validate our simplified handling of particle losses. On the other hand if the sources with $`E_{\mathrm{max}}<E_{\mathrm{norm}}`$ contribute significantly, the measurements at $`E_{\mathrm{norm}}`$ tell us nothing about their normalizations and the upper bound on neutrino fluxes may be increased up to the level dictated by the CR observations in the lower energy range, as suggested by Mannheim et al.,  (1999). This scenario is supported by Fig. Critical self-organization of astrophysical shocks, provided that there are many strong shocks in the ensemble. The observed CR steep power-law spectrum is then essentially a superposition of flatter or even non-power-law spectra from individual sources properly distributed in $`E_{\mathrm{max}}`$.

To summarize our conclusions, the main factor that should determine the particle primary spectra, and thus the neutrino flux, is how the accelerating shocks are distributed in cut-off momenta, which in a SOC state means in Mach numbers. There is no universal spectral form for individual shocks at the current state of the theory (except a not quite representative case of $`M\mathrm{}`$). Therefore, one should understand particle losses mechanisms, since they determine the shock structure and thus the spectra, directly and through $`E_{\mathrm{max}}`$. These mechanisms are inseparable from the dynamics of strong compressible MHD turbulence generated by those same particles. The further progress in its study will improve our understanding of the acceleration process and related radiation.

This work was supported by U.S. DOE under Grant No. FG03-88ER53275. We also acknowledge helpful discussions with F. A. Aharonian and J. P. Rachen.