Re-acceleration of Cosmic Ray Electrons by Multiple ICM Shocks

Since the thickness of the shock transition zone is of the order of the gyroradius of the postshock thermal protons, both protons and electrons need to be pre-accelerated to suprathermal momenta greater than the so-called injection momentum,

in order to diffuse across the shock transition layer and fully participate in the DSA process (e.g., Caprioli et al., 2015; Ryu et al., 2019; Kang, 2020). Here $p_{\rm thp}=\sqrt{2m_{\rm p}k_{\rm B}T_{2}}$ and $Q_{\rm p}\approx 3.5-3.8$ is the injection parameter (e.g., Kang & Ryu, 2010; Caprioli et al., 2015; Ha et al., 2018b). Throughout the paper, common symbols in physics are used: e.g., $m_{\rm p}$ for the proton mass, $m_{\rm e}$ for the electron mass, and $c$ for the speed of light. Adopting the traditional thermal leakage injection model (Kang et al., 2002), the distribution function of the injected/accelerated CR protons can be approximated for $p\geq p_{\rm inj}$ as

where $p_{\rm max}$ is the maximum DSA-accelerated momentum. Then the CRp injection fraction is determined by $Q_{\rm p}$ and $q$ as

Since $\xi$ is expected to increase with the shock Mach number (Ha et al., 2018b), $Q_{\rm p}$ should be smaller for higher $M$. However, the quantitative behavior of $Q_{\rm p}(M)$ based on plasma kinetic simulations has not been fully explored yet. Thus we adopt a constant value, $Q_{\rm p}=3.8$ for simplicity, since the main focus of this study is to examine the qualitatively effects of multiple shocks.

The electron injection at $Q_{\perp}$-shocks are known to involve the following key processes (e.g., Guo et al., 2014; Kang et al., 2019; Trotta & Burgess, 2019; Ha et al., 2021; Kobzar et al., 2021): (1) the reflection of some of incoming electrons at the shock ramp due to magnetic deflection, leading to the excitation of upstream waves through the EFI, (2) the generation of ion-scale waves via the AIC due to the dynamics of the reflected protons in the shock transition zone, (3) the energy gain due to the gradient-drift along the motional electric field in the shock transition zone. Thus the electron pre-acceleration occurs mainly through the so-called shock drift acceleration (SDA), rather than DSA.

Several numerical studies of electron pre-acceleration have indicated that a suprathermal tail develops with the power-law form, $p^{-q}$, for $p\gtrsim p_{\rm ref}$, which extends above the DSA injection momentum $p>p_{\rm inj}$ (Guo et al., 2014; Park et al., 2015; Trotta & Burgess, 2019; Kobzar et al., 2021). Here, $p_{\rm ref}$ represents the lowest momentum of the reflected electrons (see Figure 1 of Kang, 2020). This is again parameterized as

where $p_{\rm the}=\sqrt{2m_{\rm e}k_{\rm B}T_{2}}$ is the postshock electron thermal momentum, and so $p_{\rm inj}/p_{\rm ref}=\sqrt{m_{\rm p}/m_{\rm e}}\approx 43$ (see Figure 2). Then the spectrum of injected electrons is assumed to follow the DSA power-law for $p\geq p_{\rm ref}$:

Note that the electrons with $p_{\rm ref}\lesssim p\lesssim p_{\rm inj}$ are referred as ‘suprathermal’ electrons, whereas those with $p\gtrsim p_{\rm inj}$ are defined as CR electrons. Again the quantitative estimation for $Q_{\rm e}(M)$ has yet to come, so we adopt the same injection parameter $Q_{\rm e}=3.8$ as $Q_{\rm p}$, which results in the CRe to CRp number ratio, $K_{e/p}=(m_{\rm e}/m_{\rm p})^{(q-3)/2}$, for $p\geq p_{\rm inj}$. Hereafter, we focus on CR electrons and omit the character ‘e’ from the subscript, i.e, $f_{\rm inj,e}(p)=f_{\rm inj}(p)$. So the spectrum of CR electrons injected/accelerated at the $k$-th shock will be represented by $f_{\rm inj,k}(p)$.

From the equilibrium condition that the DSA momentum gains per cycle are equal to the synchrotron/iC losses per cycle, a maximum momentum can be estimated as follows (Kang, 2011):

Here $B_{\rm e}^{2}=B^{2}+B_{\rm rad}^{2}$ is the “effective” magnetic field strength, and $B_{\rm rad}=3.24~{}{\mu\rm G}(1+z)^{2}$ takes into account the iC cooling due to the cosmic background radiation at redshift $z$, and $B$ is given in units of $~{}{\mu\rm G}$. We set $z=0.2$ as a reference epoch, and so $B_{\rm rad}=4.7~{}{\mu\rm G}$. For typical ICM shock parameters,

The upstream spectrum of the $k$-th shock, $f_{\rm u,k}(p)$, contains the electrons injected and re-accelerated by all previous shocks, which are decompressed and cooled in the postshock region behind the $(k-1)$-th shock. Then, the downstream spectrum, $f_{\rm d,k}(p)$, re-accelerated at the $k$-th shock, can be calculated by the following re-acceleration integration (Drury, 1983):

Again we assume for simplicity that suprathermal electrons with $p_{\rm ref,k}\lesssim p\lesssim p_{\rm inj,k}$ can be re-accelerated via DSA in the same way as for $p\gtrsim p_{\rm inj,k}$, although the re-acceleration of these suprathermal electrons has not been fully explored through plasma simulations. An alternative choice for the lower bound of the integral is $p_{\rm inj,k}$, since only particles above the injection momentum could diffuse back and forth across the shock transition and fully participates in DSA. The result of the re-acceleration integral, $f_{\rm d,k}(p)$, depends somewhat weakly on the the lower bound of the integral, as illustrated in Figure 2. The exponential cutoff, $\exp(-p^{2}/p_{\rm eq}^{2})$, should be applied to Equation (8) as well.

The synchrotron emission from mono-energetic electrons with $\gamma_{\rm e}$ peaks around the characteristic frequency,

The synchrotron radiation spectrum emitted by the power-law population in Equation (5) has the power-law form, $j_{\nu}\propto\nu^{-\alpha_{\rm sh}}$, where $\alpha_{\rm sh}=(q-3)/2=0.5(M^{2}+3)/(M^{2}-1)$ (Kang, 2015).

The volume-integrated radio spectrum, $J_{\nu}$, calculated with $F(p)$ is expected to steepen to $\alpha_{\rm int}=\alpha_{\rm sh}+0.5$ for $\nu\gtrsim\nu_{\rm br}$, that corresponds roughly to the characteristic frequency of the electrons with $\gamma_{\rm e,br}$. Inserting Equation (11) into Equation (12) results in

Note that the transition of the spectral index of $J_{\nu}$ from $\alpha_{\rm sh}$ to $\alpha_{\rm sh}+0.5$ occurs rather gradually over the broad frequency range, $\nu\sim 0.01-1~{}{\rm GHz}$, because more abundant lower energy electrons also contribute to the emission at $\nu_{\rm br}$ (see Figure 1 of Carilli et al. (1991)).

As in Paper I, we consider the ICM plasma that consists of fully ionized hydrogen atoms and free electrons with $T_{1}=5.8\times 10^{7}$ K (5 keV) and $n_{1}=10^{-4}~{}{\rm cm^{-3}}$. The normalization of $f(p)$ presented in the figures below scales with $n_{1}$.

We consider the models with three shocks specified with Mach numbers, $M_{1}$, $M_{2}$, and $M_{3}$, to explore how the strength of preceding shocks affect the CR spectrum at the third shock. Effects of multiple shocks are expected to depend on the time between two consecutive shock passages and the magnetic field strength. Time intervals, $t_{\rm p1,2}=t_{\rm p1}=t_{\rm p2}\approx 5-100$ Myr are considered. The third shock is assume to last for $t_{\rm p3}\approx 100$ Myr to produce the postshock region of $\sim 100$ kpc, as observed in typical giant radio relics (e.g. van Weeren et al., 2010, 2016). The fiducial value of the preshock magnetic field is set as $B_{1}=1~{}{\mu\rm G}$. Table 1 summarizes the model parameters considered in this study.

For MHD shocks, the postshock magnetic field strength depends on $M$ and $\theta_{\rm Bn}$. In $Q_{\perp}$ shocks, $B_{2}\approx rB_{1}$, so, for instance, $B_{2}\approx 3~{}{\mu\rm G}$ for M=3 ($r=3$). Since the effective field strength in the postshock region, $B_{\rm e,2}=(B_{2}^{2}+B_{\rm rad}^{2})^{1/2}$ with $B_{\rm rad}=4.7~{}{\mu\rm G}$ ($z=0.2$), the electron cooling remains significant even for weak magnetic fields with $B_{1}\ll 1~{}~{}{\mu\rm G}$. On the other hand, the break frequency in Equation (13) scales with $B_{2}$, as long as $B_{\rm e,2}\sim 5~{}~{}{\mu\rm G}$, and so it decreases to $\nu_{\rm br}\sim 1$ MHz for $B_{1}\sim 0.01~{}~{}{\mu\rm G}$.

Figure 3 shows the CR spectrum in four models: model E ($M=3$) without cooling, with $t_{\rm p1,2}=5$ Myr, and $20$ Myr, and model A with $t_{\rm p1,2}=5$ Myr. Panels (a.1)-(d.1) show the downstream spectrum, $f_{\rm sh}=f_{\rm sum,k}$ (solid lines), and the far-downstream spectrum, $f_{\rm sum,k}^{\prime}$ (dotted lines). Panels (a.2)-(d.2) show the power-law slopes, $q_{\rm sh}=-\partial(\ln f_{\rm sh})/\partial(\ln p)$ (solid lines) and $q_{\rm sum}^{\prime}=-\partial(\ln f_{\rm sum}^{\prime})/\partial(\ln p)$ (dotted lines). Panels (a.1)-(a.2) demonstrate the effects of DSA by multiple shocks, i.e., amplification and flattening of the CR spectrum, when energy losses are ignored.

Panels (a.1)-(c.1) show that, for $\gamma_{e}\approx 10^{4}$, the amplitude of $f_{\rm sh}(p)$ increases by a factor of about $5-10$ at each passage of a $M\sim 3$ shock, while panel (d.1) shows that the amplification factor is about 2 for a $M=2.3$ shock. However, panels (b.2)-(d.2) indicate that flattening of $f_{\rm sh}(p)$ almost disappear for higher energy electrons with $\gamma_{\rm e}>10^{4}$, when postshock cooling is included. The shock slope at the third shock, $q_{\rm sh,3}$ (blue solid), is smaller (flatter) than the DSA slope, $q_{3}=3r_{3}/(r_{3}-1)$, for low-energy electrons with $\gamma_{\rm e}\lesssim 10^{4}$, retaining the flattening effect of multiple shock passages. In the case of the far-downstream spectrum, $f_{\rm sum,3}^{\prime}(p)$, behind the third shock (blue dotted lines), the flattening effects disappear for $\gamma_{\rm e}\gtrsim 10^{3}$. Thus, Figure 3 demonstrates that spectral flattening due to multiple shocks could remain significant for $t_{\rm p1,2}\lesssim 20$ Myr, while the amplification factor of CR spectrum at $\gamma_{e}\approx 10^{4}$ can range $2-10$, depending on the shock Mach number ($M\lesssim 3$).

In Figure 4, model A with $t_{\rm p1,2}=20$ Myr (left column), and $t_{\rm p1,2}=5$ Myr (middle column), and model B with $t_{\rm p1,2}=5$ Myr (right column) are presented. Panels (a.1)-(c.1) show the volume-integrated energy spectrum, $F(p)p^{4}$, at each shock and the postshock spectrum, $f_{\rm sh}(p)p^{4}$, at the third shock. Panels (a.2)-(c.2) show $q_{\rm int}=-\partial(\ln F)/\partial(\ln p)$ and $q_{\rm sh,3}=-\partial(\ln f_{\rm sh})/\partial(\ln p)$. The green solid lines display the DSA slope for the third shock, $q_{3}=3r_{3}/(r_{3}-1)$ (green solid), while the green dot-dashed lines display $q_{3}+1$ (green dot-dashed).

In the next subsection we focus on the radio synchrotron spectrum in the frequency range of $\nu_{1}\lesssim\nu\lesssim\nu_{2}$, whose emission comes mainly from electrons in the energy range of $\gamma_{1}\lesssim\gamma_{e}\lesssim\gamma_{2}$ (see Equations (12)):

where $\nu_{1}=0.1$ GHz, $\nu_{2}=10$ GHz, and $\langle B_{2}\rangle$ is the mean postshock magnetic field strength. In the upper two rows of Figure 4, the vertical black dotted lines indicate $\gamma_{1}$ and $\gamma_{2}$ for the third shock with $M_{3}$ and $B_{1}=1~{}{\mu\rm G}$.

Comparing $f_{\rm sh}(p)$ of radio emitting electrons among the three models shown in Figure 4, we find the following. (1) Model A with $t_{\rm p1,2}=5$ Myr suffers less cooling, and so retains more substantial effects of multiple shocks in the range of $\gamma_{1}\lesssim\gamma_{e}\lesssim\gamma_{2}$, compared to the model with $t_{\rm p1,2}=20$ Myr. (2) Model B with $M_{3}=3$ has a flatter spectrum than model A with $M_{3}=2.3$, but $q_{\rm sh,3}$ exhibits very little multiple shock effects in the same range of $\gamma_{e}$ (see the magenta line in panel (c.2)).

The amplitude of $F(p)$ at the first and second shocks, especially at low energies, is larger in the model with $t_{\rm p1,2}=20$ Myr than that with $t_{\rm p1,2}=5$ Myr, because the postshock advection length increases with time as $d=u_{2}t_{\rm p1,2}$. For both the models, however, $t_{\rm p3}=100$ Myr is the same, so the amplitude of $F(p)$ at the third shock (blue solid lines) is similar. As expected, $F(p)$ steepens gradually by one power of $p$ above the break momentum at $\gamma_{\rm e,br}\sim 10^{4}$ behind the third shock. So $q_{\rm int,3}$ (blue solid) approach to $q_{3}+1$, while $q_{\rm sh,3}$ (magenta dashed) approaches to $q_{3}$ at high energies. On the other hand, both the slopes are smaller (flatter) than the expected DSA slopes for $\gamma_{e}\lesssim 10^{3}-10^{5}$ due to the re-acceleration by multiple shocks. In particular, the effects of multiple shocks could persist for radio-emitting electrons with $\gamma_{1}\lesssim\gamma_{e}\lesssim\gamma_{2}$ (between the two vertical black dotted lines), if $t_{\rm p1,2}\lesssim 20$ Myr and the Mach numbers of the preceding shocks are higher than that of the last shock (i.e., model A).

The left panels of Figure 5 compare $F(p)$ for the four cases with different $B_{1}=0.01$, 0.1, 1, and 2 $~{}{\mu\rm G}$ for model A. Panels (a.1)-(a.2) show that $F(p)$ depends rather weakly on $B_{1}$, because the effective field strength, $B_{\rm e}\sim 5~{}{\mu\rm G}$, varies a little for the range of $B_{1}$ considered here.

The middle and right panels of Figure 5 show models C and D with $t_{\rm p1,2}=5$ My and $B_{1}=1~{}{\mu\rm G}$ (see Table 1). Again they demonstrate that the effects of multiple shocks are important only for the case in which the preceding shocks are stronger than the last shock (model C). In other words, the CR electrons, accelerated by weaker preceding shocks and then cooled in the postshock region, provide only low-energy seed electrons. So they could increase somewhat the amplitude of $f_{\rm sh}$, but do not affect substantially the power-law slope (model D).

The lower two rows of Figures 4 and 5 compare the radio sychrotron emission spectra for the respective models shown in the upper two rows. Panels (a.3)-(c.3) show the volume-integrated radio spectrum, $\nu J_{\rm int}(\nu)$, at each shock and the postshock spectrum, $\nu j_{\rm sh}(\nu)$, at the third shock, while panels (a.4)-(c.4) show $\alpha_{\rm int}=-d\ln J_{\rm int}/d\ln\nu$ (solid lines) and $\alpha_{\rm sh}=-d\ln j_{\rm sh}/d\ln\nu$ (magenta dashed lines). The green solid lines display the DSA slope for the third shock, $\alpha_{3}=(q_{3}-3)/2$, while the green dot-dashed lines display $\alpha_{3}+0.5$. Exceptions are panels (a.3) and (a.4) of Figure 5, where only $J_{\rm int}$ and $\alpha_{\rm int}$ are shown for model A with different values of $B_{1}$.

As noted above, the slope of $J_{\rm int}(\nu)$ increases gradually from $\alpha_{\rm sh}$ to $\alpha_{3}+0.5$ over a very broad range of the frequency (i.e., so-called CI case), where the break frequency, $\nu_{\rm br}\sim 0.3$ GHz. At high frequencies, $\alpha_{\rm int,3}$ (blue solid) approaches to $\alpha_{3}+0.5$, while $\alpha_{\rm sh,3}$ (magenta dashed) approaches to $\alpha_{3}$. On the other hand, panel (b.4) of Figure 4 shows that $\alpha_{\rm sh,3}<\alpha_{3}$ and $\alpha_{\rm int,3}<\alpha_{3}+0.5$ for $\nu\lesssim 10$ GHz, reflecting the flattening effects of multiple shocks in model A. In model C, panel (b.4) of Figure 5 exhibit the same behaviors for even higher frequencies.

Unlike $F(p)$, $J_{\rm int}(\nu)$ depends sensitively on $B_{1}$, as can be seen in panels (a.3)-(a.4) of Figure 5. This is because cooling is dominated by iC scattering off background photons for the range of $B_{1}$ considered here, while the radio synchrotron spectrum scales roughly with $B_{1}$. For model A with $B_{1}=0.01~{}~{}{\mu\rm G}$ (green solid), $J_{\rm int}(\nu)$ is almost a single power-law and $\alpha_{\rm int,3}\approx\alpha_{3}+0.5$ for $\nu>1$ GHz.

In short, the flattening of the radio spectrum at the shock, $j_{\rm sh}(\nu)$, or the volume-integrated radio spectrum, $J_{\rm int}(\nu)$, due to multiple re-acceleration is significant, only if the preceding shocks are stronger than the last shock as in models A and C. The volume-integrated spectrum $J_{\rm int}(\nu)$ steepens gradually over a very broad frequency range.

If a radio relic is generated by three shock passages as in model A, $J_{\rm int}(\nu)$ would not be a single power-law, but exhibit a spectral curvature at high frequencies, as shown in Figure 4. Then, inferring the Mach number of the radio relic shock from the relation $M_{\rm rad}=[(\alpha_{\rm int}+1)/(\alpha_{\rm int}-1)]^{1/2}$ may result in incorrect results, when the slope $\alpha_{\rm int}$ is estimated between two observation frequencies, for instance, $0.1{\rm GHz}\lesssim\nu_{\rm obs,1},\nu_{\rm obs,2}\lesssim 10{\rm GHz}$.

To examine this problem, we plot the spectral index $\alpha_{0.61}$ of $j_{\rm sh,3}(\nu)$ at 0.61 GHz versus the predicted DSA index, $\alpha_{3}=0.5(M_{3}^{2}+3)/(M_{3}^{2}-1)$, of the third shock in the left panel of Figure 6. The spectral index $\alpha_{0.15}^{1.5}$ of the volume-integrated spectrum, $J_{\rm int}$, between 0.15 and 1.5 GHz is shown against $\alpha_{3}+0.5=(M_{3}^{2}+1)/(M_{3}^{2}-1)$ of the third shock in the right panel of Figure 6. In models A and C , $\alpha_{0.61}$ is smaller than $\alpha_{3}$ due to the multiple re-acceleration effects, and the difference between the two indices is greater for smaller $t_{\rm p1,2}$. In the case of $t_{\rm p1,2}=100$ Myr (asterisks), $\alpha_{0.61}\approx\alpha_{3}$, because the multiple shock effects disappear due to postshock cooling. On the other hand, $\alpha_{0.61}\approx\alpha_{3}$ in models B and D, regardless of $t_{\rm p1,2}$, so the three symbols almost overlap each other.

By contrast, $\alpha_{\rm int}$ is smaller than $\alpha_{3}+0.5$ for all the cases shown except model A with $B_{1}=0.01~{}{\mu\rm G}$ (magenta triangle). This is because $J_{\rm int}$ steepens and exhibits spectral curvature over a very broad range of frequencies. The difference between the two indices is the greatest in model C (red symbols). For model A with $B_{1}=0.01~{}{\mu\rm G}$, the break frequency is low enough, $\nu_{\rm br}\sim 1$ MHz, so $J_{\rm int}$ becomes almost a single power-law between 150 MHz and 1.5 GHz.

This exercise illustrates that the estimation of the shock Mach number from radio spectral indices, $\alpha_{\rm sh}$ or $\alpha_{\rm int}$, at certain observation frequencies should be made with cautions, since the emission spectrum could be affected by any preceding shocks. In model A with $M_{3}=2.3$ and $t_{\rm p1,2}=5$ Myr, for instance, the Mach number estimated from $\alpha_{0.61}\approx 0.78$ would be $M_{\rm rad}\approx 2.9$, while that estimated from $\alpha_{0.15}^{1.5}\approx 1.08$ would be $M_{\rm rad}\approx 5.0$. Hence, both the estimates would be higher than the X-ray inferred value, $M_{\rm X}\approx M_{3}$.

In observations of real radio relics, however, $J_{\rm int}$ depends on the three-dimensional shape of the postshock region and the viewing angle relative to the shock surface. Modeling of more realistic configuration is beyond the scope of this study.