A Diffusive Shock Acceleration Model for Protons in Weak Quasi-parallel Intracluster Shocks

TL;DR

This paper develops an improved diffusive shock acceleration model for protons in weak quasi-parallel shocks within the intracluster medium, predicting acceleration efficiencies and emphasizing the importance of shock strength.

Contribution

It introduces a self-consistent DSA model that better estimates injection momentum and accounts for CR feedback effects in weak shocks, advancing understanding of cosmic-ray acceleration in galaxy clusters.

Findings

CR acceleration efficiency ranges from 10^-3 to 10^-2 for shocks with Mach numbers 2.25 to 5.

Proton acceleration is negligible in shocks with Mach numbers below 2.25.

The model predicts a decrease in postshock gas temperature due to CR energy transfer.

Abstract

Low sonic Mach number shocks form in the intracluster medium (ICM) during the formation of the large-scale structure of the universe. Nonthermal cosmic-ray (CR) protons are expected to be accelerated via diffusive shock acceleration (DSA) in those ICM shocks, although observational evidence for the $\gamma$-ray emission of hadronic origin from galaxy clusters has yet to be established. Considering the results obtained from recent plasma simulations, we improve the analytic test-particle DSA model for weak quasi-parallel ($Q_\parallel$) shocks, previously suggested by \citet{kang2010}. In the model CR spectrum, the transition from the postshock thermal to CR populations occurs at the injection momentum, $p_{\rm inj}$, above which protons can undergo the full DSA process. As the shock energy is transferred to CR protons, the postshock gas temperature should decrease accordingly and the subshock strength weakens due to the dynamical feed of the CR pressure to the shock structure. This results in the reduction of the injection fraction, although the postshock CR pressure approaches an asymptotic value when the CR spectrum extends to the relativistic regime. Our new DSA model self-consistently accounts for such behaviors and adopts better estimations for $p_{\rm inj}$. With our model DSA spectrum, the CR acceleration efficiency ranges $\eta\sim10^{-3}-0.01$ for supercritical, $Q_\parallel$-shocks with sonic Mach number $2.25\lesssim M_{\rm s}\lesssim5$ in the ICM. Based on \citet{ha2018b}, on the other hand, we argue that proton acceleration would be negligible in subcritical shocks with $M_{\rm s}<2.25$.

Equations

$$p_{\mathrm{inj}}=Q_{\mathrm{i}}\cdot p_{\mathrm{th},\mathrm{p}},$$

$$f_{\mathrm{CR}}(p)\approx \psi \cdot f_{\mathrm{N}}{\left(\frac{p}{p_{\mathrm{inj}}}\right)}^{-q}\exp \left[-{\left(\frac{p}{p_{\max }}\right)}^2\right].$$

$$f_{\mathrm{N}}=\frac{n_2}{{\pi }^{1.5}}{p}_{\mathrm{th},\mathrm{p}}^{-3}\exp (-Q_{\mathrm{i}}^2),$$

$$\xi \equiv \frac{4\pi }{n_2}{\int }_{p_{\min }}^{p_{\max }}{f}_{\mathrm{CR}}(p)p^{2}dp,$$

$$E_{\mathrm{CR},2}=4\pi c{\int }_{p_{\min }}^{p_{\max }}(\sqrt{p^2+(m_{p}c{)}^2}-m_{p}c)f_{\mathrm{CR}}(p)p^{2}dp.$$

$$R_{\mathrm{T}}=\frac{E_{\mathrm{th}}(T_2,0)-E_{\mathrm{CR},2}}{E_{\mathrm{th}}(T_{2,0})}.$$

Full text

Considering the results obtained from recent plasma simulations, we improve the analytic test-particle DSA model for weak quasi-parallel ($Q_{\parallel}$) shocks, previously suggested by Kang & Ryu (2010). In the model CR spectrum, the transition from the postshock thermal to CR populations occurs at the injection momentum, $p_{\rm inj}$, above which protons can undergo the full DSA process. As the shock energy is transferred to CR protons, the postshock gas temperature should decrease accordingly and the subshock strength weakens due to the dynamical feed of the CR pressure to the shock structure. This results in the reduction of the injection fraction, although the postshock CR pressure approaches an asymptotic value when the CR spectrum extends to the relativistic regime. Our new DSA model self-consistently accounts for such behaviors and adopts better estimations for $p_{\rm inj}$. With our model DSA spectrum, the CR acceleration efficiency ranges $\eta\sim 10^{-3}-0.01$ for supercritical, $Q_{\parallel}$-shocks with sonic Mach number $2.25\lesssim M_{\rm s}\lesssim 5$ in the ICM. Based on Ha et al. (2018b), on the other hand, we argue that proton acceleration would be negligible in subcritical shocks with $M_{\rm s}<2.25$.

acceleration of particles – cosmic rays – galaxies: clusters: general – shock waves

††journal: The Astrophysical Journal

1 Introduction

Hierarchical clustering of the large-scale structure of the universe induces supersonic flow motions of baryonic matter, which result in the formation of weak shocks with sonic Mach numbers $M_{\rm s}\lesssim 4$ in the hot intracluster medium (ICM) (e.g., Ryu et al., 2003; Vazza et al., 2009; Ha et al., 2018a). In particular, shocks associated with mergers of subcluster clumps have been observed in X-ray and radio (e.g., Brunetti & Jones, 2014; van Weeren et al., 2019). These ICM shocks are thought to accelerate cosmic ray (CR) protons and electrons via diffusive shock acceleration (DSA) (Bell, 1978; Drury, 1983). Although the acceleration of relativistic electrons can be inferred from the so-called giant radio relics (e.g., van Weeren et al., 2019), the presence of the CR protons produced by ICM shocks has yet to be established (e.g., Pfrommer & Enßlin, 2004; Pinzke & Pfrommer, 2010; Zandanel & Ando, 2014; Vazza et al., 2016; Kang & Ryu, 2018). The inelastic collisions of CR protons with thermal protons followed by the decay of $\pi^{0}$ produce diffuse $\gamma$-ray emission, which has not been detected so far (Ackermann et al., 2016). Previous studies using cosmological hydrodynamic simulations with some prescriptions for CR proton acceleration suggested that the non-detection of $\gamma$-ray emission from galaxy clusters would constrain the acceleration efficiency $\eta\lesssim 10^{-3}$ for ICM shocks with $2\lesssim M_{\rm s}\lesssim 5$ (e.g., Vazza et al., 2016); the acceleration efficiency is defined in terms of the shock kinetic energy flux, as $\eta\equiv E_{\rm CR,2}u_{2}/(0.5\rho_{1}u_{\rm sh}^{3})$ (Ryu et al., 2003). Hereafter, the subscripts $1$ and $2$ denote the preshock and postshock states, respectively. And $\rho$ is the density, $u$ is the flow speed in the shock-rest frame, $u_{\rm sh}$ is the shock speed, and $E_{\rm CR,2}$ is the postshock CR proton energy density.

Proton injection is one of the key processes that govern the DSA acceleration efficiency. In the so-called thermal leakage model, suprathermal particles in the tail of the postshock thermal distribution were thought to re-cross the shock from downstream to upstream and participate in the Fermi I process (e.g. Malkov, 1997; Kang et al., 2002). Through hybrid simulations, however, Caprioli & Spitkovsky (2014a, CS14a, hereafter) showed that in quasi-parallel ($Q_{\parallel}$, hereafter, with $\theta_{\rm Bn}\lesssim 45\degree$) shocks, protons are injected through specular reflection off the shock potential barrier, gaining energy via shock drift acceleration (SDA), and that the self-excitation of upstream turbulent waves is essential for multiple cycles of reflection and SDA. Here, $\theta_{\rm Bn}$ is the obliquity angle between the shock normal and the background magnetic field direction. They considered relatively strong ($M_{\rm s}\gtrsim 6.5$) $Q_{\parallel}$-shocks in plasmas with $\beta\sim 1$, where $\beta=P_{\rm gas}/P_{\rm B}$ is the ratio of the gas to magnetic pressures. As CRs are accelerated to higher energies, the CR energy density increases in time before the acceleration saturates at $E_{\rm CR,2}/(E_{\rm CR,2}+E_{\rm th,2})\approx 0.06-0.13$ for $M_{\rm s}\approx 6.3-63$ (see Figure 3 of CS14a) with the injection fraction, $\xi\sim 10^{-4}-10^{-3}$ (see Equation [4] below). Here, $E_{\rm th,2}$ is the energy density of postshock thermal protons. As a result, the postshock thermal distribution gradually shifts to lower temperatures as the CR power-law tail increases its extent (see Figure 1 of CS14a).

Moreover, CS14a found that in the immediate postshock region, the proton momentum distribution can be represented by three components: the Maxwellian distribution of thermal particles, $f_{\rm th}(p)$, the CR power-law spectrum, $f_{\rm CR}(p)$, and the suprathermal ‘bridge’ connecting smoothly $f_{\rm th}$ and $f_{\rm CR}$ (see Figure 2 of CS14a). This suprathermal bridge gradually disappears as the plasma moves further downstream away from the shock, because the electromagnetic turbulence and ensuing kinetic processes responsible for the generation of suprathermal particles decrease in the downstream region. Far downstream from the shock, the transition from the Maxwellian to CR distributions occurs rather sharply at the so-called injection momentum, which can be parameterized as $p_{\rm inj}\approx Q_{\rm i}\ p_{\rm th,p}$, where $p_{\rm th,p}=\sqrt{2m_{p}k_{B}T_{2}}$ is the postshock thermal proton momentum and $Q_{\rm i}\sim 3-3.5$ is the injection parameter. Here, $T_{2}$ is the temperature of postshock thermal ions, $m_{p}$ is the proton mass, and $k_{B}$ is the Boltzmann constant. They suggested that the CR energy spectrum can be modeled by the DSA power-law attached to the postshock Maxwellian at $p_{\rm inj}$, although their hybrid simulations revealed a picture that is quite different from the thermal leakage injection model. Later, Caprioli et al. (2015, CPS15, hereafter) presented a minimal model for proton injection that accounts for quasi-periodic shock reformation and multicycles of reflection/SDA energization, and predicted the CR spectrum consistent with the hybrid simulations of CS14a.

Recently, Ha et al. (2018b, HRKM18, hereafter) studied, through Particle-in-Cell (PIC) simulations, the early acceleration of CR protons in weak ($M_{\rm s}\approx 2-4$) $Q_{\parallel}$-shocks in hot ICM plasmas where $\beta\sim 100$ (e.g., Ryu et al., 2008). In the paper, they argued that only supercritical $Q_{\parallel}$-shocks with $M_{\rm s}\gtrsim 2.25$ develop overshoot/undershoot oscillations in their structures, resulting in a significant amount of incoming protons being reflected at the shock and injected into the DSA process. Subcritical $Q_{\parallel}$-shocks with $M_{\rm s}\lesssim 2.25$, on the other hand, have relatively smooth structures, and hence the preacceleration and injection of protons into DSA are negligible. Thus, it was suggested that ICM $Q_{\parallel}$-shocks may accelerate CR protons only if $M_{\rm s}\gtrsim 2.25$.

Although the simulations followed only to the very early stage of DSA where the maximum ion momentum reaches up to $p_{\rm max}/m_{i}c\sim 0.5$ ($m_{i}$ is the reduced ion mass111Throughout the paper, we differentiate $m_{i}$ from $m_{p}$, because the effects of the reduced mass ratio, $m_{i}/m_{e}$ ($m_{e}$ is the electron mass), in PIC simulations remain to be fully understood. In the simulations of HRKM18, for example, $m_{i}=100-800\ m_{e}$ was used.), HRKM18 attempted to quantify proton acceleration at ICM $Q_{\parallel}$-shocks. The simulated CR spectrum indicated the injection parameter of $Q_{\rm i}\approx 2.7$, which led to a rather high injection fraction, $\xi\approx 2\times 10^{-3}-10^{-2}$, for shocks with $M_{\rm s}=2.25-4$. If we simply extrapolate this injection fraction to the relativistic regime of $p_{\rm max}/m_{i}c\gg 1$, the ensuing DSA efficiency would be rather high, $\eta>0.01$, which is in strong disagreement with the existing observations of $\gamma$-rays from galaxies clusters.

In a ‘fluid-version’ of numerical studies of DSA, on the other hand, the time-dependent diffusion-convection equation for the isotropic part of the momentum distribution function, $f_{\rm CR}(p)$, is solved, adopting a Bohm-type spatial diffusion coefficient ($\kappa\propto p$) and a ‘macroscopic’ prescription for thermal leakage injection ($\tau_{\rm esp}$) (e.g., Kang et al., 2002). Previous studies using this approach managed to follow the evolution of CR proton spectrum into the relativistic energies of up to $p_{\rm max}/m_{p}c\sim 50$ for shocks with a wide range of sonic Mach numbers (e.g., Kang & Jones, 2005). They showed that, as the CR pressure increases in time, the subshock weakens and $T_{2}$ decreases accordingly, resulting in the gradual reduction of the injection rate and $f_{\rm CR}(p_{\rm inj})$ [see Figure 5 of Kang & Jones (2005)]. This leads to the decrease of the injection fraction $\xi(t)$ with time, although the postshock CR pressure reaches an approximate time-asymptotic value [see Figure 6 of Kang et al. (2002)]. These results are consistent with those of the hybrid simulations described above.

Previously, Kang & Ryu (2010) considered an analytic model for $f_{\rm CR}(p)$ in the test-particle regime of DSA for weak ICM shocks. They suggested that the test-particle solution of $f_{\rm CR}(p)$ could be valid only if $Q_{i}\gtrsim 3.8$, which results in the injection fraction $\xi\lesssim 10^{-3}$ and the CR pressure $P_{\rm CR,2}/\rho_{1}u_{\rm sh}^{2}<0.1$. In that study, however, the changes of $T_{2}(t)$ and $\xi(t)$ with the increase of $p_{\rm max}$ were not included self-consistently, because $Q_{\rm i}$, although a free parameter, has a fixed value, and $T_{2}$ was estimated simply from the Rankine-Hugoniot relation, relying on the test-particle assumption. Hence, the model failed to incorporate the full aspect of DSA observed in the previous simulations.

Based on the earlier studies of DSA using hybrid, PIC, and fluid simulations, we here propose an improved analytic model that is designed to approximately emulate the CR proton spectrum of DSA for given shock parameters. The basic formulation is still based on the test-particle solution with a thermal leakage injection recipe with a free parameter, $Q_{\rm i}$, as in Kang & Ryu (2010). The main improvement is, however, the inclusion of the reduction of the postshock thermal energy density due to the transfer of the shock energy to the CR population in a self-consistent manner; also the model considers a more realistic range of $Q_{\rm i}\approx 3.0-3.5$ that reflects the results of the hybrid simulations of CS14a and CPS15.

In the next section, we first review what has been learned about proton injection and acceleration at $Q_{\parallel}$-shocks from recent plasma simulations. In Section 3, we describe our analytic DSA model for the CR proton spectrum produced at weak $Q_{\parallel}$-shocks, along with the injection fraction and acceleration efficiency that characterize the DSA of CR protons. A brief summary follows in Section 4.

2 Implications from Plasma Simulations

Although the structure and time variation of collisionless shocks are primarily governed by the dynamics of reflected protons and the waves excited by them in the foreshock region, the roles of electron kinetic processes in proton injection to DSA has not yet been fully explored (e.g., Balogh & Truemann, 2013). Only PIC simulations can follow from first principles various microinstabilities and wave-particle interactions due to ion and electron kinetic processes. Owing to greatly disparate time and length scales of ion and electron processes, however, the runs of PIC simulations are limited to only several $\times 10^{2}\ \Omega_{\rm ci}^{-1}$, depending on $m_{i}/m_{e}$, $\beta$, and the dimension of simulations. Here, $\Omega_{\rm ci}^{-1}={m_{i}c}/{eB_{0}}$, is the ion cyclotron period where $c$ is the speed of light, $e$ is the electron charge, and $B_{0}$ is the background magnetic field strength. Typically, the injection and early acceleration of protons can be followed up to the maximum momentum of $p_{\rm max}/p_{\rm th,i}\sim 30$ ($p_{\rm th,i}=\sqrt{2m_{i}k_{B}T_{2}}$) in PIC simulations (e.g., Park et al., 2015, HRKM18).

Hybrid simulations, in which electrons are modeled as a charge-neutralizing fluid, can be run to several $\times 10^{2}-10^{3}\ \Omega_{\rm cp}^{-1}$ (where $\Omega_{\rm cp}^{-1}={m_{p}c}/{eB_{0}}$), neglecting details of electron kinetic processes. Yet they can follow proton acceleration only up to $p_{\rm max}/m_{p}u_{\rm sh}\sim 30$ or so (e.g., CP14a).

With currently available computational resources, both PIC and hybrid simulations can only study the early development of suprathermal and nonthermal protons. Thus, it would be a rather challenging task to extrapolate what we have learned about DSA from existing plasma simulations to the relativistic regime of $p_{\rm max}/m_{p}c\gg 1$.

2.1 Hybrid Simulations

As discussed in the introduction, the injection and acceleration of protons at $\beta\approx 1$ $Q_{\parallel}$-shocks with $M_{\rm s}\gtrsim 6.3$ were studied extensively through 2D hybrid simulations (CS14a and CPS15). A small fraction of incoming protons can be injected to DSA after undergoing two to three cycles of SDA, followed by reflection off the shock potential drop. In addition, at low-$\beta$ ($\beta\lesssim 1$) shocks, the proton reflection can be facilitated by the magnetic mirror force due to the compression of locally perpendicular magnetic fields in upstream MHD turbulence, which are self-excited by back-streaming protons (e.g., Sundberg et al., 2016).

The efficiency of proton injection could be quantitatively different at weak ICM shocks with $\beta\sim 100$, because the shock potential drop is smaller at lower $M_{\rm s}$ shocks and the magnetic mirror force is weaker in higher $\beta$ plasmas. Caprioli & Spitkovsky (2014b, CS14b, hereafter), on the other hand, showed that the magnetic field amplification due to resonant and non-resonant streaming instabilities increases with the Alfvén Mach number, $M_{\rm A}\approx\beta^{1/2}M_{\rm s}$. Hence, the level of upstream turbulence is expected to be higher for higher $\beta$ shocks at a given $M_{\rm s}$. Therefore, higher $\beta$ could have two opposite effects on the efficiency of proton injection, i.e., weaker magnetic mirror but stronger turbulence in the foreshock. Unfortunately, so far hybrid simulations for high-$\beta$ ($\beta\gg 1$) shocks have not been published in the literature yet.

CPS15 suggested that the proton injection at weak shocks may be different from their findings for strong shocks in the following senses: (1) the overshoot in the shock potential is smaller at weaker shocks, leading to a smaller reflection fraction at each confrontation with the shock, (2) the fractional energy gain at each SDA cycle is smaller, so more SDA cycles are required for injection, (3) the levels of turbulence and magnetic field amplification are weaker. As a result, the proton injection and acceleration efficiencies should be smaller at weaker shocks. According to Figure 3 of CS14a, for the $M_{\rm s}\approx 6.3$ shock ($M=5$ in their definition), the DSA efficiency is $\eta\approx 0.036$, so a smaller $\eta$ is expected for ICM shocks with $M_{\rm s}\lesssim 4$.

Moreover, CS14b showed in their Figure 9 that the normalization (amplitude) of postshock $f_{\rm CR}$ decreases as $p_{\rm max}(t)$ increases with time. We interpret that this trend is caused by the increase in the number of SDA cycles required for injection to DSA, because the subshock weakens gradually due to the CR feedback, and so the energy gain per SDA cycle is reduced. Considering that the ratio of $p_{\rm max}/p_{\rm th,p}$ reaches only to $\sim 30$ in these hybrid simulations, the normalization of $f_{\rm CR}$ may continue to decrease as the CR spectrum extends to the relativistic region with $p_{\rm max}/m_{p}c\gg 1$.

2.2 Particle-in-cell Simulations

HRKM18 explored for the first time the criticality of high-$\beta$ $Q_{\parallel}$ shocks and showed that protons can be injected to DSA and accelerated to become CRs only at supercritical shocks with $M_{\rm s}\gtrsim 2.25$. Figure 7 of HRKM18 showed that the shock criticality does not sensitively depend on $m_{i}/m_{e}$ and numerical resolution, but the acceleration rate depends slightly on $\beta$. As mentioned before, turbulence is excited more strongly for higher $\beta$ cases due to higher $M_{\rm A}$. But the reflection fraction is smaller for higher $\beta$ due to weaker magnetic mirror forces, leading to lower reflection fraction and lower amplitude of $f_{\rm CR}$ near $p_{\rm inj}$.

In order to get a glimpse of the long-term evolution of the CR proton spectrum, we extend the 1D PIC simulation reported in HRKM18 from $\Omega_{\rm ci}t_{\rm end}=90$ to 270 for the model of $M_{\rm s}=3.2$, $\theta_{\rm Bn}=13^{\circ}$, $m_{i}/m_{e}=100$, $\beta=100$, and $T_{1}=8.6$ keV ($10^{8}$ K). Details of numerical and model setups can be found in HRKM18 (see their Table 1). The main change is that a different computation domain, $[L_{x},L_{y}]=[3\times 10^{4},1]~{}(c/w_{\rm pe})^{2}$, is adopted here in order to accommodate the longer simulation time. Because of severe computational requirements, in practice, it is difficult to extend this kind of PIC simulations to a much larger box for a much longer duration. In this simulation the average velocity of ions is $\sqrt{18.36}$ times higher than that of real protons for the given temperature.

Figure 1 shows the time evolution of the postshock energy spectra of ions, $dN/d\gamma$ (where $\gamma$ is the Lorentz factor), and the injection fraction, $\xi(t)$ [see Eq. (11) of HRKM18]. We adopt the simulation data of HRKM18 for $\Omega_{\rm ci}t\approx 94$ (red), while the data from the new extended simulation is used for $\Omega_{\rm ci}t\approx 240$ (blue). The region of $(1.5-2.5)r_{L,i}$ behind the shock is included, where $r_{L,i}$ is the ion Larmor radius defined with the incoming flow speed. Note that the spectrum near the energy cutoff might not be correctly reproduced due to the limited size of the simulation domain.

We notice the following features in Figure 1(a) : (1) the postshock temperature decreases slightly with time, (2) the injection parameter, $Q_{\rm i}=p_{\rm inj}/p_{\rm th,i}$, increases from $\sim 2.7$ to $\sim 3.0$ as the time increases from $\Omega_{\rm ci}t\approx 90$ to 240, and (3) the amplitude of $dN/d\gamma(\gamma_{\rm inj})$ decreases gradually.

Figure 1(b) shows the resulting gradual decrease of $\xi(t)$, which may continue further in time. As in HRKM18, a somewhat arbitrary value of $p_{\rm min}=\sqrt{2}p_{\rm inj}$ is adopted (see the next section for a further discussion). We interpret the bump in the evolution of $\xi(t)$ near $\Omega_{\rm ci}t\approx 210$ as a consequence of shock reformation.

3 Analytic Model for CR Proton Spectrum

The analytic model presented here inherits the test-particle DSA model with a thermal leakage injection recipe, which was suggested by Kang & Ryu (2010). It describes the downstream CR proton spectrum, $f_{\rm CR}(p)$, for weak shocks. For the preshock gas with the density, $n_{1}$, and the temperature, $T_{1}$, the postshock vales, $n_{2}$, and $T_{2,0}$222Here, $T_{2,0}$ denotes the temperature of the thermal gas when the postshock CR energy density, $E_{\rm CR,2}$, is negligible, reserving $T_{2}$ for the cases of non-negligible $E_{\rm CR,2}$., can be calculated from the Rankine–Hugoniot jump condition. For example, the shock compression ratio is given as $r=n_{2}/n_{1}=(\gamma_{\rm g}+1)/(\gamma_{\rm g}-1+2/M^{2}_{s})$, where $\gamma_{\rm g}=5/3$ is the gas adiabatic index.

Following Kang & Ryu (2010), we parameterize the model as follows: (1) The CR proton spectrum follows the test-particle DSA power-law, as $f_{\rm CR}(p)\propto p^{-q}$, where $q=3r/(r-1)$. (2) The transition from the postshock thermal to CR spectra occurs at the injection momentum

[TABLE]

where $Q_{\rm i}$ is the injection parameter. The main improvement here is that the postshock temperature, $T_{2}$, decreases slightly from $T_{2,0}$, and hence $p_{\rm th,p}$ does too, as the fraction of the shock energy transferred to CRs increases.

Our model leads to the following form of the CR proton spectrum,

[TABLE]

Here, the maximum momentum of CR protons, $p_{\rm max}$, increases with the shock age (e.g., Kang & Ryu, 2010). The normalization factor can be approximated as

[TABLE]

assuming the CR power-law spectrum is hinged to the postshock Maxwell distribution at $p_{\rm inj}$. Therefore, in our model, $Q_{\rm i}$ is the key parameter that controls $f_{\rm N}$. In addition, we introduce an additional parameter, $\psi\sim 1$, to accommodate any uncertainties in determining the value of $Q_{\rm i}$ and the resulting amplitude, $f_{\rm N}$. Throughout this paper, however, $\psi=1$ is used. Figure 2(a) shows the model spectrum, $f_{\rm CR}(p)$, calculated with Equations (2)-(3), which illustrates the transition from the thermal to nonthermal CR spectra at $p_{\rm inj}$.

The PIC simulation described in section 2.2 indicates $Q_{\rm i}\approx 3$ when $p_{\rm max}/m_{i}c\approx 0.5$, but $Q_{\rm i}$ may further increase for $p_{\rm max}/m_{i}c\gg 1$, as noted above. On the other hand, hybrid simulations for strong shocks of $\beta\approx 1$ (CS14a, CPS15) showed that it is expected to range as $Q_{\rm i}\approx 3.0-3.5$. As discussed in section 2.1, higher $\beta$ could have two opposite effects on proton reflection, weaker magnetic mirror but stronger upstream turbulence. So it is difficult to make quantitative predictions on the long-term evolution of $Q_{i}$ in high-$\beta$ shocks without performing plasma simulations of very long duration. Here, we will consider the range of $Q_{\rm i}=3.3-3.5$ as an educated guess from the previous plasma simulation. Moreover, in our analytic model, $Q_{\rm i}<3.3$ would give the DSA efficiency of $\eta\gtrsim 0.01$ for $3\lesssim M_{\rm s}\lesssim 5$ (see Figure 4 below), which would be incompatible with the non-detection of $\gamma$-ray emission from galaxy clusters (Vazza et al., 2016).

From the model $f_{\rm CR}(p)$ in Equation (2), we calculate the injection fraction of CR protons by

[TABLE]

as in HRKM18. The postshock CR energy density is estimated by

[TABLE]

In the case of very weak shocks, where the CR spectrum is dominated by low energy particles, both $\xi$ and $E_{\rm CR,2}$ depend sensitively on the lower bound of the integrals, $p_{\rm min}$ (e.g., Pfrommer & Enßlin, 2004). We here adopt $p_{\rm min}\approx p_{\rm inj}$ for fiducial models, while $p_{\rm min}=780\ {\rm MeV}/c$, the threshold energy of $\pi$-production reaction, will be considered as well for comparison.

As mentioned above, $E_{\rm CR,2}$ may increase, as $f_{\rm CR}(p)$ extends to higher $p_{\rm max}$, resulting in the decrease of the postshock gas temperature from $T_{2,0}$ to $T_{2}$ (see Figure 5 of Kang & Jones (2005) and Figure 1 of CS14a). Thus, we introduce the temperature reduction factor,

[TABLE]

Then, $T_{2}=R_{\rm T}T_{2,0}$ is the reduced postshock temperature.333The fraction of thermal particles that becomes CR protons is assumed to be small, i.e., $\xi\ll 1$.

CPS15 suggested that when the postshock CR energy density approaches to $E_{\rm CR,2}\approx 0.1E_{\rm sh}=0.1(\rho_{1}u_{\rm sh}^{2}/2)$, the subshock weakens substantially, which suppresses the proton reflection and injection. Hence, the normalization of $f_{\rm CR}$ is expected to decrease as $p_{\rm max}$ increases. Our model is designed to mimic such a behavior by finding the self-consistent postshock thermal distribution with a lower temperature, while $p_{\rm inj}$ is assumed to be fixed. Then, the injection parameter increases as $Q_{\rm i}=Q_{\rm i,0}/\sqrt{R_{\rm T}}$, where $Q_{\rm i,0}$ is the initial value, leading to smaller values of $f_{\rm N}$. Note that $p_{\rm inj}$ at shocks with different parameters ($M_{s}$, $\theta_{\rm Bn}$, and $\beta$) is controlled by a number of complex kinetic process, and hence should be studied through long-term plasma simulations, beyond the current computational capacity. Considering that the proton injection into DSA is yet to be fully understood, fixing $p_{\rm inj}$ while slightly increasing $Q_{\rm i}$ in our model should be regarded as a reasonable assumption.

Figure 2(a) shows the model spectrum, including that of the self-consistent thermal distribution, in a $M_{\rm s}=3.2$ shock for $Q_{\rm i,0}=3.0-3.5$; the spectrum depends on the adopted value of $Q_{\rm i,0}$. Panels (b)-(d) illustrate the change of the model spectrum as $p_{\rm max}$ increases in shocks with $M_{\rm s}=2.5,$ 3.2, and 4.0, respectively. As $p_{\rm max}$ and also $E_{\rm CR,2}$ increase, the Maxwellian part shifts to slightly lower $T_{2}$, and $R_{\rm T}$ decreases accordingly. Because $p_{\rm inj}$ is assumed to be fixed, $Q_{\rm i}$ increases and thus the normalization factor $f_{\rm N}$ decreases in our model.

Figure 3 shows the change of $\xi$, $R_{\rm T}$, $E_{\rm CR,2}/E_{\rm sh}$, and $\eta$, calculated with Equations (2)-(6), as $p_{\rm max}$ increases for $Q_{\rm i,0}=3.3$ (dashed lines) and 3.5 (solid lines) in shocks with $M_{\rm s}=2.25-4.0$. The CR acceleration efficiency is related to the postshock CR energy density, as $\eta=E_{\rm CR,2}/rE_{\rm sh}$. Figure 3(b) plots how $R_{T}$ decreases, as $p_{\rm max}$ increases.

The injection fraction, $\xi$, increases with increasing $p_{\rm max}$ during the early acceleration phase, but decreases for $p_{\rm max}/m_{p}c\gg 1$. The latter behavior results from the gradual reduction of $f_{\rm CR}(p_{\rm inj})$, which is caused by the self-adjustment of the shock structure, that is, the cooling of the postshock thermal protons, the growing of the precursor, and the weakening of the subshock due to the dynamical feedback of the CR pressure. $E_{\rm CR,2}/E_{\rm sh}$ and $\eta$, on the other hand, monotonically increase and approach to asymptotic values for $p_{\rm max}/m_{p}c\gtrsim 10^{2}$.

Figure 4 shows the asymptotic values of those quantities as a function of $M_{\rm s}$ (filled circles and lines) for $Q_{\rm i,0}=3.3$ (black) and 3.5 (red), which would cover the most realistic range for ICM shocks (CS14a). As mentioned in the introduction, HRKM18 showed that ICM $Q_{\parallel}$-shocks with $M_{\rm s}<2.25$ may not inject protons into the DSA process, resulting in inefficient CR proton acceleration. We here include the $M_{\rm s}=1.5$ and $2.0$ cases (connected with dotted lines) for illustrative purposes, showing the values estimated with our model.

Note that the asymptotic value of $\xi(M_{\rm s})$ decreases with increasing $M_{\rm s}$ for supercritical shocks with $M_{\rm s}\geq 2.25$. This behavior is opposite to the relation, $\xi\propto M_{\rm s}^{1.5}$, during the very early acceleration stage of the PIC simulations reported in HRKM18. In those PIC simulations, $p_{\rm max}/m_{p}c\lesssim 0.5$, so the CR feedback effect is not very significant. However, our analytic model is designed to take account for the dynamic feedback of the CR pressure to the shock structure when $p_{\rm max}/m_{p}c\gg 1$, so $\xi$ could be smaller at higher $M_{\rm s}$.

With the adopted value of $Q_{\rm i,0}=3.3-3.5$, $E_{\rm CR,2}/E_{\rm sh}<0.1$, so the test-particle assumption should be valid. The acceleration efficiency increases with $M_{\rm s}$ and is close to $\eta\approx 0.01-0.02$ in the range of $M_{\rm s}=3-5$. Obviously, if the injection parameter is larger than what the hybrid simulations of CS14a indicated, that is, $Q_{\rm i,0}>3.5$, then DSA would be even less efficient.

In the studies of $\gamma$-ray emission from simulated galaxy clusters, the lower bound of $f_{\rm CR}$ is often taken as $p_{\rm min}=780~{}{\rm MeV}/c$, as noted above. The open triangles in Figure 4 show $E_{\rm CR,2}/E_{\rm sh}$ and $\eta$ calculated with this $p_{\rm min}$, otherwise adopting the same analytic spectrum given in Equations (2)-(3). For $M_{\rm s}=2.25$, the acceleration efficiency with $p_{\rm min}=780~{}{\rm MeV}/c$ is smaller by a factor of 3.3 than that with $p_{\rm min}=p_{\rm inj}$. But the two estimations are similar for $M_{\rm s}\gtrsim 4$. The efficiency with $p_{\rm min}=780~{}{\rm MeV}/c$ is $\eta\sim 0.01$ in the range of $M_{\rm s}=3-5$, while $\eta\sim 10^{-3}$ for $M_{\rm s}=2.25$. If this result is extended to the case of $M_{\rm s}\sim 6$, $\eta$ would be still close to 0.01, which is about three times smaller than the efficiency reported by CS14a (i.e., $\eta\approx 0.036$ at $M_{\rm s}\sim 6.3$). Note that this estimate is somewhat larger than the upper limit of $\eta\lesssim 10^{-3}$, quoted to be consistent with the non-detection of $\gamma$-ray emission from galaxy clusters by Vazza et al. (2016). On the other hand, as HRKM18 shown, $\eta$ may be very small and negligible for shocks with $M_{\rm s}<2.25$, for which the fraction of the total shock dissipation in the ICM was shown to be substantial (e.g., Ryu et al., 2003). Hence, the consistency of our model for proton acceleration with the non-detection of cluster $\gamma$-rays should be further examined by considering the details of the characteristics of shocks in simulated galaxy clusters.

4 Summary

The DSA efficiency for CR protons at low $M_{\rm s}$ $Q_{\parallel}$-shocks in the high-$\beta$ plasmas of the ICM has yet to be investigated through kinetic plasma simulations. HRKM18 studied the injection and the early acceleration of protons up to $p_{\rm max}/m_{i}c\approx 0.5$ at such shocks through 1D PIC simulations, adopting reduced mass ratios of $m_{i}/m_{e}$. On the other hand, CS14a, CS14b, and CPS15 carried out hybrid simulations to study the DSA of protons, but considered only high $M_{\rm s}$ shocks in $\beta\approx 1$ plasmas. Here, we revisited the test-particle DSA model for low $M_{\rm s}$ shocks with a thermal leakage injection recipe that was previously presented in Kang & Ryu (2010). Reflecting new findings of recent plasma simulations, we improved the analytic DSA model by accounting for the transfer of the postshock thermal energy to the CR energy and the weakening of the subshock due to the dynamical feedback of the CR pressure to the shock structure.

We first set up an approximate analytic solution, $f_{\rm CR}(p)$, for CR protons in weak $Q_{\parallel}$-shocks. We then calculated the injection fraction, $\xi$, the postshock CR energy fraction, $E_{\rm CR,2}/E_{\rm sh}$, and the acceleration efficiency, $\eta$, of CR protons. The main aspects of our model and the main results are summarized as follows.

1. In weak shocks with $M_{\rm s}\lesssim 5$, above the injection momentum, $p_{\rm inj}=Q_{i}\ p_{\rm th,p}$, $f_{\rm CR}(p)$ follows the test-particle DSA power-law, whose slope is determined by the shock compression ratio.

2. According to plasma simulations such as CS14a, CPS15, and HRKM18, as CR protons are accelerated to higher energies, the postshock gas temperature $T_{2}$ and the normalization of $f_{\rm CR}$ decreases (see Figure 2). Thus, in our model, while the injection momentum, $p_{\rm inj}$, is assumed to be fixed, the injection parameter increases as $Q_{i}=Q_{i,0}/\sqrt{R_{\rm T}}$, where $R_{\rm T}$ is the reduction factor of the postshock temperature. Then $Q_{i}$ determines the CR spectrum according to Equations (2)-(6). We adopt $Q_{i,0}\approx 3.3-3.5$, extrapolating the results of previous hybrid simulations.

3. In our model, as $f_{\rm CR}(p)$ extends to higher $p_{\rm max}/m_{p}c\gg 1$, $\xi$ first increases and then decreases due to the reduction of $T_{2}$ and the increase of $Q_{\rm i}$, although $\eta$ monotonically increases and approaches a time-asymptotic value. Such a behavior was previously seen in fluid DSA simulations (e.g., Kang et al., 2002).

4. Both $\xi$ and $E_{\rm CR,2}/E_{\rm sh}$ depend on $Q_{i,0}$ and also the lower bound of the integrals, $p_{\rm min}$, especially in the case of very weak shocks (see Figure 4). For $p_{\rm min}\approx p_{\rm inj}$ and $Q_{i,0}=3.5$, the CR acceleration efficiency ranges as $\eta\approx 3.5\times 10^{-3}-0.01$ for $2.25\lesssim M_{\rm s}\lesssim 5.0$. If $p_{\rm min}\approx 780~{}{\rm MeV}/c$ is adopted, it decreases to $\eta\approx 1.1\times 10^{-3}-0.01$ for the same Mach number range. If $Q_{i,0}=3.3$ is adopted, $\eta$ becomes larger by a factor of $1.5-2$, compared to the case with $Q_{i,0}=3.5$.

5. In subcritical shocks with $M_{\rm s}<2.25$, protons may not be efficiently injected into DSA, so we expect that $\eta$ would be negligible at these very weak shocks (HRKM18).

In a parallel paper (Ha et al., 2019), we will investigate the $\gamma$-ray emission as well as the neutrino emission from simulated galaxy clusters due to the inelastic collisions of CR protons and ICM thermal protons, based on the analytic CR proton spectrum proposed in this paper. In particular, we will check whether the prediction for $\gamma$-ray emission complies with the upper limits imposed by Fermi LAT observations.

We thank the anonymous referee for critical comments that help us improve this paper from its initial form. D.R. and J.-H. H. were supported by the National Research Foundation of Korea (NRF) through grants 2016R1A5A1013277 and 2017R1A2A1A05071429. H.K. was supported by the Basic Science Research Program of the NRF through grant 2017R1D1A1A09000567.