Shock Acceleration Model for the Toothbrush Radio Relic

TL;DR

This paper investigates the origin of the Toothbrush radio relic by comparing weak and strong shock models, incorporating turbulent acceleration, and finds strong-shock models with low-energy seed electrons are more plausible.

Contribution

It introduces a comprehensive analysis of DSA models including turbulent acceleration to explain the relic's properties, favoring strong-shock scenarios with low-energy seed electrons.

Findings

Both models fit the integrated radio spectrum well.

Turbulent acceleration explains the broad transverse profile.

Strong-shock models with low-energy seed electrons are preferred.

Abstract

Although many of the observed properties of giant radio relics detected in the outskirts of galaxy clusters can be explained by relativistic electrons accelerated at merger-driven shocks, significant puzzles remain. In the case of the so-called Toothbrush relic, the shock Mach number estimated from X-ray observations ($M_{\rm X}\approx1.2-1.5$) is substantially weaker than that inferred from the radio spectral index ($M_{\rm rad}\approx2.8$).Toward understanding such a discrepancy, we here consider the following diffusive shock acceleration (DSA) models:(1) weak-shock models with $M_{\rm s}\lesssim 2$ and a preexisting population of cosmic-ray electrons (CRe) with a flat energy spectrum,and (2) strong-shock models with $M_{\rm s}\approx3$ and either shock-generated suprathermal electrons or preexisting fossil CRe. We calculate the synchrotron emission from the accelerated CRe, following the time evolution of the electron DSA, and subsequent radiative cooling and postshock turbulent acceleration (TA). We find that both models could reproduce reasonably well the observed integrated radio spectrum of the Toothbrush relic, but the observed broad transverse profile requires the stochastic acceleration by downstream turbulence, which we label "turbulent acceleration" or TA to distinguish it from DSA. Moreover, to account for the almost uniform radio spectral index profile along the length of the relic, the weak-shock models require a preshock region over 400~kpc with a uniform population of preexisting CRe with a high cutoff energy ($\gtrsim 40$ GeV). Due to the short cooling time, it is challenging to explain the origin of such energetic electrons. Therefore, we suggest the strong-shock models with low-energy seed CRe ($\lesssim 150$~MeV) are preferred for the radio observations of this relic.

Tables

Table 1. Table 1: Model Parameters

Model	$k{T}_1$	$M_{s,i}$	${\left[M_{s,o}\right]}^a$	${\left[k{T}_{2,o}\right]}^b$	$B_1$	${\left[B_{2,o}\right]}^c$	$s$	${\gamma }_{e,c}$	${\tau }_{\mathrm{acc},0}$	Remarks
Name	(keV)			(keV)	($\mu G$)	($\mu G$)			(Myr)
W1.7a	$5.2$	$1.7$	$1.64$	$8.56$	$1.5$	$2.7$	$4.4$	${10}^5$	$100$	no injection
W1.7b	$5.2$	$1.7$	$1.64$	$8.56$	$1.5$	$2.7$	$4.4$	$4\times {10}^4$	$100$	no injection
W1.7c	$5.2$	$1.7$	$1.64$	$8.56$	$1.5$	$2.7$	$4.4$	${10}^4$	$100$	no injection
W1.7aN	$5.2$	$1.7$	$1.64$	$8.56$	$1.5$	$2.7$	$4.4$	$8\times {10}^4$	-	no injection
W2.0a	$4.3$	$2.0$	$1.87$	$8.23$	$1.5$	$2.5$	$4.4$	$8\times {10}^4$	$100$	no injection
W2.0b	$4.3$	$2.0$	$1.87$	$8.23$	$1.5$	$2.5$	$4.4$	$4\times {10}^4$	$100$	no injection
W2.0c	$4.3$	$2.0$	$1.87$	$8.23$	$1.5$	$2.5$	$4.4$	${10}^4$	$100$	no injection
W2.0d	$4.3$	$2.0$	$1.87$	$8.23$	$1.5$	$2.5$	$4.4$	$8\times {10}^4$	$50$	no injection
W2.0aN	$4.3$	$2.0$	$1.87$	$8.23$	$1.5$	$2.5$	$4.4$	$8\times {10}^4$	-	no injection
S3.6a	$3.0$	$3.6$	$3.03$	$11.2$	$1$	$2.5$	$4.6$	$3\times {10}^2$	$100$	$\kappa =1.6$
S3.6b	$3.0$	$3.6$	$3.03$	$11.2$	$1$	$2.5$	$4.6$	$3\times {10}^2$	$100$	seed CRe
S3.6c	$3.0$	$3.6$	$3.03$	$11.2$	$1$	$2.5$	$4.6$	$3\times {10}^2$	$100$	no decay ($r_{dec}\to \infty$)
S3.6aN	$3.0$	$3.6$	$3.03$	$11.2$	$1$	$2.5$	$4.6$	$3\times {10}^2$	-	$\kappa =1.6$
S3.6bN	$3.0$	$3.6$	$3.03$	$11.2$	$1$	$2.5$	$4.6$	$3\times {10}^2$	-	seed CRe

Equations

$$\Delta l_{\nu }\approx 120\mathrm{kpc}\left(\frac{u_2}{10^3\mathrm{km}{\mathrm{s}}^{-1}}\right)\cdot Q\cdot {\left[\frac{{\nu }_{\mathrm{obs}}(1+z)}{0.61\mathrm{GHz}}\right]}^{-1/2},$$

$$Q(B_2,z)\equiv \left[\frac{(5\mu \mathrm{G}{)}^2}{B_2^2+B_{\mathrm{rad}}(z{)}^2}\right]{\left(\frac{B_2}{5\mu \mathrm{G}}\right)}^{1/2},$$

$$\frac{\partial g_{\mathrm{e}}}{\partial t}+u\frac{\partial g_{\mathrm{e}}}{\partial r}$$

$$=\frac{1}{3r^2}\frac{\partial (r^{2}u)}{\partial r}\left(\frac{\partial g_{\mathrm{e}}}{\partial y}-4g_{\mathrm{e}}\right)+\frac{1}{r^2}\frac{\partial }{\partial r}\left[r^2\kappa (r,p)\frac{\partial g_{\mathrm{e}}}{\partial r}\right]$$

$$+p\frac{\partial }{\partial y}\left[\frac{D_{pp}}{p^3}\left(\frac{\partial g_{\mathrm{e}}}{\partial y}-4g_{\mathrm{e}}\right)\right]+p\frac{\partial }{\partial y}\left(\frac{b}{p^2}{g}_{\mathrm{e}}\right),$$

$$\kappa (r,p)={\kappa }^{\ast }\left(\frac{p}{m_{\mathrm{e}}c}\right),$$

$$t_{\mathrm{rad}}({\gamma }_{\mathrm{e}})=9.8\times 10^7\mathrm{yr}{\left[\frac{(5\mu \mathrm{G}{)}^2}{B^2+B_{\mathrm{rad}}(z{)}^2}\right]}^2{\left(\frac{{\gamma }_{\mathrm{e}}}{10^4}\right)}^{-1}.$$

$$D_{pp}=\frac{p^2}{4{\tau }_{\mathrm{acc}}},$$

$${\tau }_{\mathrm{acc}}\sim {\left(\frac{c}{a}\right)}^2\frac{1}{⟨k⟩c}\frac{P}{W_f}.$$

$${\tau }_{\mathrm{acc}}={\tau }_{\mathrm{acc},0}\cdot \exp \left[\frac{(r_{\mathrm{s}}-r)}{r_{\mathrm{dec}}}\right]$$

$$f_{\mathrm{inj}}(r_{\mathrm{s}},p)=f_N{\left(\frac{p}{p_{\mathrm{inj}}}\right)}^{-q}\exp \left[-{\left(\frac{p}{p_{\mathrm{eq}}}\right)}^2\right],$$

$$p_{\mathrm{eq}}={\gamma }_{eq}{m}_{\mathrm{e}}c=\frac{m_{\mathrm{e}}^2{c}^2{u}_{\mathrm{s}}}{\sqrt{4e^{3}q/27}}{\left(\frac{B_1}{B_{e,1}^2+B_{e,2}^2}\right)}^{1/2}{k}_{\mathrm{Bohm}}^{-1},$$

$$f_{\mathrm{reacc}}(r_{\mathrm{s}},p)=q\cdot p^{-q}{\int }_{p_{\mathrm{inj}}}^p{p}^{\prime q-1}{f}_{\mathrm{pre}}(p^{\prime })d{p}^{\prime }$$

$$f_{\mathrm{pre}}(p)=f_o\cdot p^{-s}\exp \left[-{\left(\frac{p}{p_{e,c}}\right)}^2\right],$$

$$I_{\nu }(R)={\int }_0^{h_{1,\max }}{j}_{\nu }(r)d{h}_1+{\int }_0^{h_{2,\max }}{j}_{\nu }(r)d{h}_2,$$

Full text

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

††journal: The Astrophysical Journal††software: CRASH (Kang & Jones, 2006)

1 INTRODUCTION

Some galaxy clusters contain diffuse, peripheral radio sources on scales as large as $\sim$2 Mpc in length, called ‘giant radio relics’ (see, e.g., Feretti et al., 2012; Brüggen et al., 2012; Brunetti & Jones, 2014, for reviews). Typically they show highly elongated morphologies, radio spectra relatively constant along the length of the relic, but steepening across the width, and high linear polarization (Enßlin et al., 1998; van Weeren et al., 2010). Moreover, they have integrated radio spectra with a power-law form at low frequencies, but that steepen above gigaherts frequencies (Stroe et al., 2016). Previous studies have demonstrated that such observational features can often be understood as synchrotron emission from $\gtrsim 10$ GeV electrons in $\sim\mu$G magnetic fields, accelerated via diffusive shock acceleration (DSA) at merger-driven shock waves in the cluster periphery (e.g., Kang et al., 2012).

Yet, significant questions remain in the merger-shock DSA model of radio relics. Three of the troublesome issues are (1) low DSA efficiencies predicted for electrons injected in situ and accelerated at weak, $M_{\rm s}\lesssim 3$, shocks that are expected to form in merging clusters (e.g., Ryu et al., 2003; Kang et al., 2012); (2) inconsistencies of the X-ray based shock strengths with radio synchrotron-based shock strengths with the X-ray measures typically indicating weaker shocks (e.g., Akamatsu & Kawahara, 2013; Ogrean et al., 2014); and (3) a low fraction ($\lesssim 10$ %) of observed merging clusters with detected radio relics (e.g., Enßlin & Gopal-Krishna, 2001; Kang, 2016a). According to structure formation simulations, the mean separation between shock surfaces is $\sim 1$ Mpc, and the mean lifetime of intracluster medium (ICM) shocks is $t_{\rm dyn}\sim 1$ Gyr (e.g., Ryu et al., 2003; Pfrommer et al., 2006; Skillman et al., 2008; Vazza et al., 2009). So, actively merging clusters are expected to contain several shocks, and we might actually expect multiple radio relics in typical systems. Some of these difficulties could be accounted for in a scenario in which a shock may light up as a radio relic only when it encounters a preexisting cloud of fossil relativistic electrons in the ICM (e.g., Enßlin, 1999; Kang & Ryu, 2015).

Here, we focus on issue (2) above. To set up what follows, we note that in the test-particle DSA model for a steady, planar shock, nonthermal electrons that are injected in situ and accelerated at a shock of sonic Mach number $M_{\rm s}$ form a power-law momentum distribution, $f_{\rm e}(p,r_{\rm s})\propto p^{-q}$ with $q=4M_{\rm s}^{2}/(M_{\rm s}^{2}-1)$ (Drury, 1983). Based on this result and the relation $\alpha_{\rm sh}=(q-3)/2$ between $q$ and the synchrotron spectral index at the shock, $\alpha_{\rm sh}$ (with $j_{\nu}\propto\nu^{-\alpha}$), the Mach number of the hypothesized relic-generating shock is then commonly inferred from its radio spectral index using the relation $\alpha_{\rm sh}=(M_{\rm rad}^{2}+3)/2(M_{\rm rad}^{2}-1)$. On the other hand, the shock Mach number can be also estimated from the temperature discontinuity obtained from X-ray observations, using the shock jump condition, $T_{2}/T_{1}=(M_{\rm X}^{2}+3)(5M_{\rm X}^{2}-1)/16/M_{\rm X}^{2}$, where the subscripts, $1$ and $2$, identify the upstream and downstream states, respectively. Sometimes, if the temperature jump is poorly constrained, $M_{\rm X}$ is assessed from an estimate of the density jump, $\sigma=\rho_{2}/\rho_{1}=4M_{\rm X}^{2}/(M_{\rm X}^{2}+3)$. Although the radio and X-ray shock measures can agree, sometimes the synchrotron index, $\alpha_{\rm sh}$, implies a significantly higher Mach number, $M_{\rm rad}$, than $M_{\rm X}$.

Without subsequent, downstream acceleration, the effects of synchrotron and inverse Compton (iC) “cooling” ($\dot{p}\propto-p^{2}$) will truncate the postshock electron spectrum above energies that drop with increasing distance from the shock, since cooling times for $\gtrsim 10$ GeV electrons under cluster conditions are generally $<100$ Myr (e.g., Brunetti & Jones, 2014). That is the standard explanation for observed spectral steepening across the width of the relic (downstream of the shock). For reference, we note that if the shock is steady and planar, and the postshock magnetic field is uniform, this energy loss prescription translates into an integrated synchrotron spectral index, $\alpha_{\rm int}=\alpha_{\rm sh}+1/2$ (Heavens & Meisenheimer, 1987).

The spectral index along the northern, “leading” edge of the head portion (B1) of the so-called Toothbrush relic in the merging cluster 1RXS J060303.3 at $z=0.225$ is estimated to be $\alpha_{\rm sh}\approx 0.8$ ($q\approx 4.6$) with the corresponding radio Mach number $M_{\rm rad}\approx 2.8$ (van Weeren et al., 2016). But the gas density jump along the same edge in B1 inferred from X-ray observations implies a much weaker shock with $M_{\rm X}\approx 1.2-1.5$ (van Weeren et al., 2016). The associated radio index, $\alpha_{\rm sh}\sim 2-5$ ($q\sim 7-13$), is much too steep to account for the observed radio spectrum.

Toward understanding this discrepancy between $M_{\rm rad}$ and $M_{\rm X}$ for the Toothbrush relic, we here consider the following two scenarios for modeling the radio observations of this relic: (1) weak-shock models ($M_{\rm s}\lesssim 2$) with flat-spectrum, preexisting cosmic-ray elections (CRe), and (2) strong-shock models ($M_{\rm s}\approx 3$) with low-energy seed CRe. In the weak-shock models, we adopt a preshock, preexisting population of CRe with the “right” power-law slope, for example, $f_{\rm pre}(p)\propto p^{-s}\exp[-(p/p_{e,c})^{2}]$ with $s=2\alpha_{\rm sh}+3\sim 4.4$, where $p_{e,c}/m_{\rm e}c>10^{4}$ is an effective cutoff to the spectrum. In the strong-shock models, on the other hand, $M_{\rm s}\approx M_{\rm rad}$ is chosen to match the observed radio spectral index, and low-energy seed CRe ($p/m_{\rm e}c\sim 30$) are assumed to come from either the suprathermal tail population generated at the shock or a preexisting fossil population.

In the weak-shock models, the value for $p_{e,c}$ is critical, since the observed emissions at frequencies $\gtrsim 100$ MHz generally come from electrons with 10 GeV or higher energies ($p/m_{\rm e}c\gtrsim 10^{4}$). So, if $p_{e,c}/m_{\rm e}c\ll 10^{4}$, the preexisting electron population provides just “low-energy seed electrons” to the DSA process, which for $M_{\rm s}\lesssim 2$ would still lead to $\alpha_{\rm sh}\gtrsim 1.2$, and cannot produce the observed radio spectrum with $\alpha_{\rm sh}\approx 0.8$. Hence, the weak-shock models with $M_{\rm s}\lesssim 2$ can reproduce the observed spectral index profile of the Toothbrush relic only with $p_{e,c}/m_{\rm e}c\sim 7-8\times 10^{4}$ (Kang, 2016a). Consequently, in order to explain the fact that $M_{\rm rad}>M_{\rm X}$ for the Toothbrush or similar relics by the weak-shock models, one should adopt a potentially radio-luminous, preexisting electron population with $p_{e,c}/m_{\rm e}c\gg 10^{4}$. Preshock radio emission might be observable in this case, of course, unless the shock has already swept through the region containing fossil electrons, so no preshock electrons remain.

Such a requirement for large $p_{e,c}$, however, poses a question about the origin of these preexisting CRe, since the electrons with $p/m_{\rm e}c\sim 8\times 10^{4}$ in a $~{}{\mu\rm G}$ field cool on a brief time scale of $t_{\rm rad}\sim 10$ Myr. In the B1 portion of the Toothbrush relic, the spectral index is observed to be uniform over 400 kpc along its leading edge. This means that in the weak-shock models the length scale of the preshock region containing preexisting CRe with a uniformly flat spectrum with large $p_{e,c}$ should be as long as 400 kpc, tangential to the shock surface. Moreover, since the shock compression ratio is $\sigma\gtrsim 2$ for $M_{\rm s}\gtrsim 1.5$, while the observed radio width of the head is at least $150$ kpc, the width of this uniform preshock region should be $\gtrsim 300$ kpc along the shock normal direction. Considering the short cooling times for high-energy electrons, it should be difficult to explain the origin of such a uniform cloud of preexisting CRe by fossil CRe that were deposited in the past by an active galactic nucleus (AGN) jet, for instance, unless the effective electron dispersion speed across the preshock structure was $\gtrsim 0.1$c, or there was a uniformly effective turbulent acceleration (TA) in effect across that volume.

In the strong-shock models with $M_{\rm s}\approx 3$, on the other hand, the challenge to account for the uniform spectrum at the relic edge becomes less severe, since the models require only low-energy seed CRe ($p/m_{\rm e}c\sim 30$) that could be provided by either the shock-generated suprathermal electrons or preexisting fossil CRe ($p/m_{\rm e}c\lesssim 300$) with long cooling times ($t_{\rm rad}\gtrsim 3.5$ Gyr). In the latter case, the additional requirement for low-energy preexisting CRe to enhance the radio emission may explain why only a fraction of merger shocks can produce radio relics. Such low-energy CRe may originate from previous episodes of shock or turbulence acceleration or AGN jets in the ICM (e.g., Enßlin, 1999; Pinzke et al., 2013).

According to simulations for the large-scale structure formation of the universe, the surfaces of merger-driven shocks responsible for radio relics are expected to consist of multiple shocks with different $M_{\rm s}$ (see, e.g., Skillman et al., 2008; Vazza et al., 2009). From mock X-ray and radio observations of relic shocks in numerically simulated clusters, Hong et al. (2015) showed that the shock Mach numbers inferred from an X-ray temperature discontinuity tend to be lower than those from radio spectral indices. This is because X-ray observations pick up the part of shocks with higher shock energy flux but lower $M_{\rm s}$, while radio emissions come preferentially from the part with higher $M_{\rm s}$ and so higher electron acceleration. In the strong-shock models, we assume that the B1 portion of the Toothbrush relic represents a portion of the shock surface with $M_{\rm s}\approx 3$, extending over 400 kpc along the length of the relic.

It is important to note that the transverse width across the B1 component of the Toothbrush relic is about two times larger than that of another well-studied radio relic, the so-callled Sausage relic in CIZA J2242.8+5301. The FWHM at 610 MHz, for example, is about 110 kpc for the Toothbrush relic (van Weeren et al., 2016), while it is about 55 kpc for the Sausage relic (van Weeren et al., 2010). For the high-frequency radio emission from electrons with $p/m_{\rm e}c\sim 10^{4}$ radiatively cooled downstream from the shock, the characteristic width of the relic behind a spherical shock is

[TABLE]

where $u_{2}$ is the flow speed immediately downstream of the shock and $z$ is the redshift of host clusters (Kang, 2016b). We will argue below that relic-producing merger shocks are more spherical in geometry than planar. Note, then, that since the downstream flow speed in the shock rest frame increases behind a spherical shock, the advection length in a given time scale is somewhat longer than that estimated for a planar shock (Donnert et al, 2016). The factor $Q$ depends on the postshock magnetic field strength, $B_{2}$, as

[TABLE]

where $B_{2}$ and $B_{\rm rad}=3.24~{}{\mu\rm G}(1+z)^{2}$ are expressed in units of $\mu$G. The factor $Q$ evaluated, for instance, for $z=0.225$ peaks with $Q_{\rm max}\approx 0.6$ with $B_{2}\approx 2.8~{}{\mu\rm G}$. Then, with $u_{2}\approx 10^{3}~{}{\rm km~{}s^{-1}}$, the maximum width at 610 MHz becomes $\Delta l_{\nu}\approx 65$ kpc. Being only about half the observed width of the B1 region of the Toothbrush relic at this frequency, it seems too small to allow the width to be set by radiative cooling alone following acceleration at the shock surface.

To overcome such a mismatch, we here consider and include the process in which electrons are additionally accelerated stochastically by MHD/plasma turbulence behind the shock, that is, TA. Along somewhat similar lines, Fujita et al. (2015) recently suggested that radio spectra harder than predicted by the DSA in a weak shock could be explained if relativistic electrons are reaccelerated through resonant interactions with strong Alfvénic turbulence developed downstream of the relic shock. However, on small scales, Alfvénic MHD turbulence is known to become highly anisotropic, so resonant scattering is weak and ineffective at particle acceleration (e.g., Brunetti & Lazarian, 2007). On the other hand, fast-mode compressive turbulence remains isotropic down to dissipation scales, so it has become favored in treatments of stochastic reacceleration of electrons producing radio halos during cluster mergers (Brunetti & Lazarian, 2007, 2011). We emphasize, at the same time, that solenoidal turbulence, likely to be energetically dominant on large scales, could still play a reacceleration role through turbulent magnetic reconnection (e.g., Brunetii & Lazarian, 2016) or generation of small-scale slow-mode MHD waves that might interact resonantly with CRe (e.g., Lynn et al., 2014). In our study we do not depend on TA to produce a flat electron spectrum at the shock, but rather explore its potential role as an effective means of slowing energy loss downstream of the shock.

In the next section, we describe our numerical simulations incorporating both DSA and TA for shock-based models designed to explore this problem. In Section 3, our results are compared with the observations of the Toothbrush relic. A brief summary follows in Section 4.

2 NUMERICAL CALCULATIONS

2.1 DSA Simulations for 1D Spherical Shocks

According to cosmological simulations (e.g., Ryu et al., 2003; Vazza et al., 2009; Hong et al., 2014), the formation and evolution of cluster shocks can be quite complex and transient with time scale $\lesssim 1$ Gyr, but the overall morphologies of shock surfaces could be represented by partial surfaces of spherical bubbles blowing outward. As in Kang & Ryu (2015), we here attempt to follow for $\lesssim 0.2$ Gyr the evolution of a 1D spherical shock, which accounts for deceleration and adiabatic expansion behind the shock (see the Section 2.4 for the details).

In our simulations, the diffusion-convection equation for a relativistic electron population is solved in 1D spherical geometry:

[TABLE]

where $g_{\rm e}(r,p,t)=f_{\rm e}(r,p,t)p^{4}$ is the pitch-angle-averaged phase space distribution function of electrons, $r$ is the radial distance from the cluster center, and $y\equiv\ln(p/m_{\rm e}c)$, with the electron mass, $m_{\rm e}$, and the speed of light, $c$ (Skilling, 1975). The background flow velocity, $u(r,t)$, is obtained by solving the usual gas dynamic conservation equations in the test-particle limit where the nonthermal pressure is dynamically negligible.

The spatial diffusion coefficient for relativistic electrons is assumed to have the following Bohm-like form:

[TABLE]

where $\kappa^{*}=k_{\rm Bohm}\cdot m_{\rm e}c^{3}/(3eB)$ and $k_{\rm Bohm}\geq 1$, with the limiting value representing Bohm diffusion for relativistic particles.

The electron energy loss term, $b(p)=\dot{p}_{\rm Coul}+\dot{p}_{\rm sync+iC}$, takes account of Coulomb scattering, synchrotron emission, and iC scattering off the cosmic microwave background (CMB) radiation (e.g., Sarazin, 1999). For a thermal plasma with the number density $n_{\rm th}$, the Coulomb cooling rate is $\dot{p}_{\rm Coul}=3.3\times 10^{-29}n_{\rm th}[1+\ln(\gamma_{\rm e}/n_{\rm th})/75]$, while the synchrotron-iC cooling rate is $\dot{p}_{\rm sync+iC}=3.7\times 10^{-29}(\gamma_{\rm e}/10^{4})^{2}[(B/3.24\mu{\rm G})^{2}+(1+z)^{4}]$ in cgs units, where $z$ is the redshift. Hereafter, the Lorentz factor, $\gamma_{\rm e}=p/m_{\rm e}c$, will also be used for relativistic energy. Note that Coulomb cooling was not considered in our previous studies for DSA modeling of radio relics (e.g., Kang et al., 2012; Kang & Ryu, 2015; Kang, 2016a, b). However, since $\dot{p}_{\rm Coul}\gtrsim\dot{p}_{\rm sync+iC}$ for $\gamma_{\rm e}\lesssim 100$ in the cluster outskirts with $n_{\rm th}\approx 10^{-4}~{}{\rm cm^{-3}}$, while $t_{\rm Coul}\sim p_{\rm e}/\dot{p}_{\rm Coul}\lesssim$ Gyr for $\gamma_{\rm e}\lesssim 10$, Coulomb cooling can, in some cases, significantly affect the electron spectrum for $\gamma_{\rm e}\lesssim 10^{4}$ and also the ensuing radio emissivity at $0.1-1$ GHz.

The radiative cooling time due to synchrotron-iC losses is given by

[TABLE]

For $B=2.5~{}{\mu\rm G}$ and $z=0.225$, for example, $t_{\rm rad}\approx 8.2\times 10^{7}~{}{\rm yr}(\gamma_{\rm e}/10^{4})^{-1}$.

In order to explore the effects of stochastic acceleration by turbulence, TA, we include the momentum diffusion term and implement the Crank-Nicholson scheme for it in the momentum space into the existing CRASH numerical hydrodynamics code (Kang & Jones, 2006). Our simulations all assume a gas adiabatic index $\gamma_{g}=5/3$. Any nonthermal pressures from CRe and magnetic fields are dynamically insignificant in our models (see below), so they are neglected. The physical nature of the CRe momentum diffusion coefficient $D_{pp}$ is discussed in the following section.

2.2 Momentum Diffusion due to Turbulent Acceleration

We pointed out in the Introduction that the recent observations of van Weeren et al. (2016) showed that (1) the transverse FWHMs of the B1 Toothbrush component are 140 kpc at 150 MHz and 110 kpc at 610 MHz, and (2) the spectral index between the two frequencies increases from $\alpha_{150}^{610}\approx 0.8$ at the northern edge to 1.9 at approximately 200 kpc to the south, toward the cluster center. While the systematic spectral steepening suggests postshock electron cooling, these widths are much broader than the cooling length given in Equation (1). In effect, the spectral steepening due to radiative cooling in the postshock region is inconsistent with the observed profiles of radio fluxes and spectral index in this region, unless the effect of cooling is somehow substantially reduced (see Section 3).

In response, we explore a scenario in which the postshock electrons gain energy from turbulent waves via Fermi II acceleration, TA, thus mitigating spectral steepening downstream. Turbulence accelerates particles stochastically; that is, if the characteristic momentum shift in a collision is $\Delta p$ and the characteristic scattering time interval is $\Delta t$, then the resulting momentum diffusion coefficient is $D_{pp}\sim(\Delta p)^{2}/\Delta t$. Since scattering events in turbulence typically lead to $\Delta p\propto p$, a convenient general form is

[TABLE]

where $\tau_{\rm acc}\sim(1/4)\langle(p/\Delta p)^{2}\Delta t\rangle$ is an effective acceleration time scale. If $\tau_{\rm acc}$ is independent of momentum, this form with the factor 4 inserted into Equation (3) leads to $\tau_{\rm acc}=\langle p\rangle/(\partial\langle p\rangle/\partial t)$, where $\langle p\rangle$ is the mean momentum of the distribution, $f_{\rm e}(r,p,t)$.

Generally speaking, TA in an ICM context can include nonresonant scattering off compressive hydrodynamical (acoustic) turbulence (e.g., Ptuskin, 1988), as well as gyro-resonant scattering off Alfvenic turbulence (e.g., Fujita et al., 2015) and Landau (also known as Cerenkov or “transit time damping”, TTD) resonance off compressive MHD turbulence (with accompanying micro-instabilities to maintain particle isotropy; (e.g., Brunetti & Lazarian, 2007, 2011; Lynn et al., 2014)). Resonant acceleration will most often be faster than nonresonant acceleration (e.g., Brunetti & Lazarian, 2007; Miniati, 2015). Alfvenic gyro-resonance involves turbulent wavelengths comparable to particle Larmour radii, which in ICM conditions for the CRe energies of interest will be sub-astronomical unit scale. While solenoidal turbulence may very well dominate the turbulence of interest (e.g., Porter et al., 2015) and, in the form of Alfven waves, probably cascades to sufficiently small scales (e.g., Kowal & Lazarian, 2010), it should become highly anisotropic on small scales in ICM settings and thus very inefficient in resonant scattering of CRe (e.g., Yan & Lazarian, 2002)111We mention for completeness a proposed alternate scenario in which solenoidal turbulence leads to fast magnetic reconnection and produces a hybrid, first-second order reacceleration process (Brunetii & Lazarian, 2016).. Fast-mode, compressive MHD turbulence should remain isotropic to dissipative scales, however. So, even though the magnetic energy in the waves of this mode will be relatively less, they can be much more effective accelerators.

On these grounds, we adopt for our exploratory calculations a simple TA model based on TTD resonance with compressive, isotropic fast-mode MHD turbulence. Assuming that in the medium $\beta_{p}=P/P_{B}\gg 1$, where $P$ is the plasma thermal pressure and $P_{B}=B^{2}/(8\pi)$ is the magnetic pressure, we can then roughly express the acceleration time, $\tau_{\rm acc}$, as

[TABLE]

Here, $a$ is the acoustic wave speed, $W_{f}$ is the total energy density in fast-mode turbulence (mostly contained in compressive “potential energy,” but also including transverse magnetic fields essential for resonant scattering). The term $\langle k\rangle$ measures the power-spectrum-weighted mean wavenumber of the fast-mode turbulence (e.g., Brunetti & Lazarian, 2007). For a power spectrum, ${\cal P}_{f}(k)\propto k^{-\alpha}$, over the range $2\pi/L_{0}\leq k\leq 2\pi/\ell_{d}$, with $3/2\leq\alpha\leq 2$ (e.g., Brunetti & Lazarian, 2011) and $\langle 1/k\rangle=(L_{0}/2\pi)\mathcal{H(\alpha})$ with $\sqrt{\ell_{d}/L_{0}}\leq\mathcal{H}\leq 1/\ln{L_{0}/\ell_{d}}$. We can roughly estimate an outer scale, $L_{0}\sim 100$ kpc behind the shock of interest. The fast-mode dissipation scale, $\ell_{d}$, is uncertain and dependent on plasma collisionality, but it is likely to be less than $\sim 1$ kpc (e.g., Schekochihin & Crowley, 2006; Brunetti & Lazarian, 2007, 2011). Putting these together, we can estimate $1/(\langle k\rangle c)\sim 10^{4}$ yrs. With an acoustic speed, $a\sim 10^{3}~{}{\rm km~{}s^{-1}}$, and an estimate $W_{f}\sim(1/10)P$ for shock-enhanced fast-mode turbulence in the immediate postshock flow, we obtain a rough estimate of $\tau_{\rm acc,0}\sim 100$ Myr near the shock. As a simple model allowing for decay of this turbulence behind the shock, we apply the form $W_{f}\propto\exp\left[{{-(r_{\rm s}-r)}/r_{\rm dec}}\right]$ with $r_{\rm s}>r$, where $r_{\rm s}$ is the radius of the spherical shock. So the TA time scale increases behind the shock as

[TABLE]

with, in most of our simulations, $r_{\rm dec}\approx 100$ kpc.

2.3 DSA Solutions at the Shock

Since the time scale for DSA at the shock is much shorter than the cooling time scale for radio-emitting electrons ($\sim 100$ Myr), we assume that electrons are accelerated almost instantaneously to the maximum energy at the shock front. On the other hand, the minimum diffusion length scale to obtain converged solutions in simulations for Equation (3) is much smaller than the typical downstream cooling length of $\sim 100$ kpc. Taking advantage of such disparate scales, we adopt analytic solutions for the electron spectrum at the shock location as $f(r_{\rm s},p)=f_{\rm inj}(p)$ or $f_{\rm reacc}(p)$, while Equation (3) is solved outside the shock. Here, $f_{\rm inj}(p)$ represents the electrons injected in situ and accelerated at the shock, while $f_{\rm reacc}(p)$ represents the reaccelerated electrons preexisting in the preshock region. So, basically we follow the energy losses and TA of electrons behind the shock, while the DSA analytic solutions are applied to the zone containing the shock. Note that shocks in CRASH are true discontinuities and tracked on sub-grid scales (Kang & Jones, 2006). Since we do not need to resolve the diffusive shock precursor or follow the DSA process in detail, this scheme allows us to use a much coarser grid, reducing dramatically the required computation time.

The electron population injected in situ from the background plasma and accelerated by DSA at the shock is modeled as

[TABLE]

where $f_{N}$, $q$, $p_{\rm inj}$, and $p_{\rm eq}$ are the normalization factor, the standard test-particle DSA power-law slope, the injection momentum, and the cutoff momentum, respectively. The injection momentum roughly identifies particles with gyro-radii large enough to allow a significant fraction of them to recross the physical shock from downstream rather than being advected downstream (e.g., Gieseler et al., 2000; Kang et al., 2002; Caprioli et al., 2015). So particles with $p>p_{\rm inj}$ are assumed to participate in the Fermi I acceleration process. According to the hybrid simulations by Caprioli & Spitkovsky (2014), $p_{\rm inj}\approx(3-3.5)p_{\rm th,p}$ for protons at quasi-parallel shocks, where $p_{\rm th,p}=\sqrt{2m_{p}kT_{2}}$ is the proton thermal momentum and $k$ is the Boltzmann constant.

The electron injection to the DSA Fermi I process from the thermal pool is thought to be very inefficient, since the momentum of thermal electrons ($p_{\rm th,e}=\sqrt{2m_{\rm e}kT_{2}}$) is much smaller than $p_{\rm inj}$. Recent particle-in-cell (PIC) simulations of quasi-perpendicular shocks by Guo et al. (2014), however, showed that some of the incoming electrons are specularly reflected at the shock ramp and accelerated via multiple cycles of shock drift acceleration (SDA), resulting in a suprathermal, power-law-like tail. Those suprathermal electrons are expected to be injected to the full Fermi I acceleration and eventually accelerated to highly relativistic energies. Such a hybrid process combining specular reflection with SDA and DSA between the shock ramp and upstream waves is found to be effective at both quasi-perpendicular and quasi-parallel collisionless shocks (Park et al., 2015; Sunberg et al., 2016). However, the injection momentum for electrons is not well constrained, since the development of the full DSA power-law spectrum extending to $p/m_{\rm e}c\gg 1$ has not been established in the simulations due to severe computational requirements for these PIC plasma simulations. Here, we adopt a simple model in which the electron injection depends on the shock strength as $p_{\rm inj}\approx(6.4/\sigma)m_{p}u_{\rm s}$, in effect resulting in $p_{\rm inj}\sim 150p_{\rm th,e}$. For a smaller compression ratio, the ratio, $p_{\rm inj}/m_{p}u_{\rm s}$, is larger, so the injection becomes less efficient.

The factor $f_{N}$ in Equation (9) depends on the suprathermal electron population with $p\sim p_{\rm inj}$ in the background plasma. We assume that the background electrons are energized via kinetic plasma processes at the shock and form a suprathermal tail represented by a $\kappa$ distribution of $\kappa=1.6-2.5$, rather than a Maxwellian distribution. The $\kappa$ distribution is well motivated in collisionless plasmas such as those in ICMs, where nonequilibrium interactions can easily dominate for the distribution of suprathermal particles (Pierrard & Lazar, 2010). It has a power-law-like high-energy tail, which asymptotes to the Maxwellian distribution for large $\kappa$. The relatively large population of suprathermal particles enhances the injection fraction compared to the Maxwellian form (Kang et al., 2014). This enhancement is larger for smaller $\kappa$. The injection efficiency at the shock is also less sensitive to the shock Mach number, compared to that from the Maxwellian distribution. Note, however, that the suprathermal electron population and the injection rate do not affect significantly the shapes of the radio-emitting electron energy spectrum and the ensuing radio synchrotron spectrum, so the adopted models for $p_{\rm inj}$ and the $\kappa$ distribution do not influence the main conclusions of this study.

The cutoff momentum in Equation (9) can be estimated from the condition that the DSA acceleration rate is equal to the synchrotron/iC loss rate:

[TABLE]

where $u_{\rm s}$ is the shock speed and $B_{\rm e}^{2}=B^{2}+B_{\rm rad}(z)^{2}$ represents the effective magnetic field strength that accounts for both synchrotron and iC losses (Kang, 2011). For typical parameters with $u_{\rm s}\sim 3\times 10^{3}~{}{\rm km~{}s^{-1}}$, $B_{1}\sim 1~{}{\mu\rm G}$, and $k_{\rm Bohm}\sim 1$, the cutoff momentum becomes $p_{\rm eq}/m_{\rm e}c\sim 10^{8}$, but the exact value is not important, as long as $p_{\rm eq}/m_{\rm e}c\gg 10^{4}$.

If there is a preexisting, upstream electron population, $f_{\rm pre}(p)$, the accelerated population at the shock is given by

[TABLE]

(Drury, 1983). In previous studies, the DSA of preexisting CR particles is commonly referred to as “reacceleration” (e.g., Kang et al., 2012; Pinzke et al., 2013), so we label $f_{\rm reacc}$ as the “DSA reaccelerated” component. In contrast, $f_{\rm inj}$ in Equation (9) represents the DSA of the background suprathermal particles injected in situ at the shock, so we label it as the “DSA injected” component. We emphasize that our DSA reacceleration models involve irreversible acceleration of preexisting CRe, in contrast to the adiabatic compression models of Enßlin et al. (1998).

In our simulations, the preshock electron population is assumed to have a power-law spectrum with exponential cutoff as follows:

[TABLE]

where the slope $s$ is chosen to match the observed radio spectral index. As mentioned in the Introduction, we adopt a large cutoff Lorentz factor, $\gamma_{e,c}=p_{e,c}/m_{\rm e}c=10^{4}-10^{5}$, in the weak-shock models, while $\gamma_{e,c}=300$ in the strong-shock models (see also Table 1). The normalization factor, $f_{o}$, is arbitrary in the simulations, since the CR pressure is dynamically insignificant (that is, in the test-particle limit). Yet, it would be useful to parameterize it with the ratio of the CRe to the gas pressure in the preshock region, $N\equiv P_{\rm CRe,1}/P_{1}\propto f_{o}$ for a given set of $s$ and $p_{e,c}$. In the models considered here, typically $N\sim(0.05-0.5)\%$ matches the amplitude of observed radio flux in the Toothbrush relic.

In our DSA simulations, $p_{e,c}$ is assumed for simplicity to stay constant in the preshock region for the duration of the simulations ($\sim 200$ Myr). This is probably unrealistic for high-energy electrons with $\gamma_{e}>10^{4}$ (see Equation (5)), unless preexisting electrons are accelerated continuously in the preshock region, for instance by turbulence.

2.4 Model Parameters

2.4.1 Observed Properties of the Toothbrush Relic

Before outlining our simulation model parameters, we briefly review our target, the Toothbrush radio relic. The relic has a linear morphology aligned roughly east-west with multiple components that, together, resemble the head and handle of a toothbrush (van Weeren et al., 2012) on respectively the west and east ends. Our focus is on the head component (labeled as B1 in Figure 4 of van Weeren et al. (2012)), whose “bristles” point southward and whose northern edge seems to coincide with the shock location detected in X-ray observations. van Weeren et al. (2016) estimated rather similar preshock and postshock temperatures, $kT_{1}=8.3_{-2.4}^{+3.2}$ keV and $kT_{2}=8.2_{-0.9}^{+0.7}$ keV, respectively, indicating that $kT_{1}$ is more uncertain from their data. From the slope change in the X-ray surface brightness across the putative shock in the component B1, they estimated a low shock Mach number, $M_{\rm X}\sim 1.2$. On the other hand, $M_{\rm rad}\approx 2.8$ is required to explain the radio spectral index ($\alpha_{\rm s}\approx 0.8$) at the northern edge of B1 as a consequence of the DSA of CRe electrons injected locally from the thermal plasma.

2.4.2 Shock Dynamics

We assume for simplicity, but one step beyond a planar shock model, that the shock dynamics can be approximated initially by a self-similar blast wave that propagates through an isothermal ICM with the density profile of $n_{\rm th}=10^{-4}~{}{\rm cm^{-3}}(r/0.8{\rm Mpc})^{-2}$. Then, the shock radius and velocity evolve roughly as $r_{\rm s}\propto t^{2/3}$ and $u_{\rm s}\propto t^{-1/3}$, respectively, where $t$ is the time since the nominal point explosion for the spherical blast wave (e.g., Ryu & Vishniac, 1991). The shock Mach number decreases in time as the spherical shock expands in the simulations. For this self-similar shock, the downstream flow speed in the upstream rest frame decreases toward the cluster center as $u(r)\propto(r/r_{\rm s})$, so the postshock flow speed with respect to the shock front increases downstream away from the shock. We acknowledge that the actual shock dynamics in the simulations deviate slightly from such behaviors, since the model shocks are not strong, although this should not influence our conclusions.

Table 1 summarizes model parameters for the DSA simulations considered in this study. Considering the observed ranges for both $kT_{1}$ and $kT_{2}$, we vary the preshock temperature as $kT_{1}=3.0-5.2$ keV. At the onset of the simulations, the shock is specified by the initial Mach number, $M_{s,i}=1.7-3.6$, which sets the initial shock speed as $u_{\rm s,i}=M_{\rm s,i}\cdot 150~{}{\rm km~{}s^{-1}}\sqrt{T_{1}/10^{6}{\rm K}}$, and is located at $r_{\rm s,i}\approx 0.8$ Mpc from the cluster center. This can be regarded as the time when the relic-generating shock encounters the preshock region containing preexisting electrons, that is, the birth of the radio relic.

We define the “shock age,” $t_{\rm age}\equiv t-t_{\rm onset}$, as the time since the onset of our simulations. We find that the downstream radio flux profiles and the integrated spectrum become compatible with the observations at the “time of observation,” $t_{\rm age}\sim 140-150$ Myr, typically when the shock is located at $r_{\rm s}\approx 1.1-1.2$ Mpc. The fourth and fifth columns of Table 1 show the shock Mach number, $M_{\rm s,o}$, and the postshock temperature, $kT_{\rm 2,o}$, at the time of observation.

In this study, we examine if the various proposed DSA-based models can explain the observed radio flux profiles reported by van Weeren et al. (2016), which, as we pointed out, depend strongly on the electron cooling length behind the shock. Therefore the magnetic field strength, which impacts electron cooling, is another key parameter. The sixth column of Table 1 shows the preshock magnetic field strength, $B_{1}=1-1.5~{}{\mu\rm G}$, which is assumed to be uniform in the upstream region. The postshock magnetic field strength is modeled as $B_{2}(t)=B_{1}\sqrt{1/3+2\sigma(t)^{2}/3}\approx 2.5-2.7~{}{\mu\rm G}$, which decreases slightly as the shock compression ratio, $\sigma(t)$, decreases in time in response to shock evolution. For the downstream region ($r<r_{\rm s}$), we assume a simple model in which the magnetic field strength scales with the gas pressure as $B_{\rm dn}(r,t)=B_{2}(t)\cdot[P(r,t)/P_{2}(t)]^{1/2}$, where $P_{2}(t)$ is the gas pressure immediately behind the shock.

2.4.3 DSA Model Parameters

As mentioned in the Introduction, the discrepancy between the observationally inferred values of $M_{\rm X}$ and $M_{\rm rad}$ could be resolved if we adopt a preexisting electron population with the “right values” of $s$ and $p_{e,c}$. Alternatively, we can explain the observed radio spectral index with a shock with $M_{\rm rad}$, assuming that $M_{\rm X}$ and $M_{\rm rad}$ may represent different parts of nonuniform shock surfaces. Our study considers both of these possible scenarios: (1) in the weak-shock models a shock with $M_{\rm s}\lesssim 2$ encounters a preshock region of a flat preexisting CRe population with $\gamma_{e,c}>10^{4}$, and (2) in the strong-shock models a shock with $M_{\rm s}\approx 3.0$ accelerates low-energy seed electrons ($\gamma_{\rm e}\sim 30$), either shock-generated suprathermal electrons or preexisting fossil CRe.

In general, we find in these experiments that the models in which postshock electrons cool without turbulent reacceleration cannot explain the broad widths of the observed radio flux profiles, independent of the assumed shock strength and CRe sources, as shown in the next section. Consequently, we also explore models that include postshock TA with the characteristic acceleration time scale of $\tau_{\rm acc}\sim 100$ Myr, which, as argued in Section 2.2, is justifiable in this context and also is comparable to expected postshock electron cooling times.

To facilitate the analyses below, we comment briefly on the model naming convention in Table 1. The first character, W or S, refers to weak-shock or strong-shock models, respectively, while the number after the first letter corresponds to the initial Mach number, $M_{s,i}$. This is followed by a sequence label (a, b, c, d) as the preexisting CRe cutoff, $\gamma_{e,c}$, or TA time, $\tau_{\rm acc}$, parameters vary. If there is no postshock TA, the letter “N” is appended at the end.

In the weak-shock models, we adopt the initial shock Mach number, $M_{s,i}=1.7-2.0$, and set $s=4.4$ as the power-law slope for preexisting CRe. In order to see the dependence of emissions on the cutoff energy in the preexisting electron spectrum, we consider a wide range of $\gamma_{e,c}=10^{4}-10^{5}$ in the W1.7a, b, c and W2.0a, b, c models (column 9 of Table 1). In model W2.0d, an enhanced, postshock turbulent reacceleration with shorter $\tau_{\rm acc,0}$ is considered. For all W1.7 and W2.0 models, the in situ injection from the background plasma is turned off in order to focus on the “DSA reacceleration” of preexisting CRe.

In the case of the strong-shock scenario, the S3.6a model includes only the “DSA injection” from a suprathermal $\kappa$ distribution of $\kappa=1.6$, while the S3.6b model incorporates only the “DSA reacceleration” of the preexisting CRe population with $s=4.6$ and $\gamma_{e,c}=300$. For the S3.6b model, the simulation results remain similar for different values of cutoff energy, $\gamma_{e,c}$, as long as $\gamma_{e,c}>p_{\rm inj}/m_{\rm e}c\approx 30$. In the S3.6c model, the decay of turbulence is turned off ($r_{dec}\rightarrow\infty$), so the momentum diffusion coefficient is assumed to be uniform behind the shock; that is, $D_{pp}=p^{2}/(4\tau_{\rm acc,0})$.

The upper panels of Figure 1 show the preexisting electron spectrum, $f_{\rm pre}$ (red and black dotted lines), and the analytic solutions for the shock spectra, $f_{\rm reacc}$ and $f_{\rm inj}$, given in Equations (11) and (9), respectively. Here, the normalization for $f_{\rm pre}$ corresponds to $N\simeq 0.01$ for W1.7a and W2.0a and $N\simeq 0.001$ for S3.6b. For the W1.7 and W2.0 models, at the shock $g_{\rm e}(r_{\rm s},p)=p^{4}f_{\rm reacc}(p)$ is used, since the in situ injection from the background plasma is suppressed. For these models, the slope of $f_{\rm reacc}(p)$ at the shock position is the preshock, $s$, for $p<p_{\rm e,c}$, while it becomes the DSA value, $q$, for $p>p_{\rm e,c}$.

In the upper right panel of Figure 1, the black dot-dashed line illustrates the $\kappa$ distribution of $\kappa=1.6$ for $p<p_{\rm inj}$, while the black solid and red long-dashed lines show $f_{\rm inj}(p)$ and $f_{\rm reacc}(p)$, respectively, for $p\geq p_{\rm inj}$. As shown here, the normalization factor $f_{N}$ for $f_{\rm inj}(p)$ is specified by the $\kappa$ distribution. In all S3.6 models, the DSA slope $q$ is flatter than $s$, so both $f_{\rm inj}$ and $f_{\rm reacc}$ have power-law spectra with the slope $q$, extending to $p_{\rm eq}/m_{\rm e}c\sim 10^{8}$, independent of $\gamma_{e,c}$. As a result, preexisting low-energy CRe just provide seeds to the DSA process and enhance the injection, but do not affect the shape of the postshock electron spectrum for $p\gg p_{\rm inj}$ (i.e., the black solid line for S3.6a versus the red dashed line for S3.6b in the upper right panel).

Regarding the shock-generated suprathermal electron population and its posited non-Maxwellian, $\kappa$-distribution form, the $\kappa$ index is not universal, since it depends on a local balance of nonequilibrium processes. If we adopt a steeper $\kappa$ distribution, with, for example, $\kappa=2.5$ (blue dot-dashed line in Figure 1), then the amplitude of the injected electron flux at $p_{\rm inj}$ will be smaller, and so the ensuing radio flux will be reduced from the models shown here (S3.6a and S3.6aN).

3 RESULTS OF DSA SIMULATIONS

3.1 Radial Profiles of Radio Emissivity

Figure 2 shows the evolution of the synchrotron volume emissivity at 150 MHz, $j_{150}(r)$, and the associated spectral index between 150 and 610 MHz, $\alpha_{150}^{610}(r)$, determined from $j_{150}(r)$ and $j_{610}(r)$. The shock is located at $r_{s,i}\approx 0.8$ Mpc at the start of the simulations, $t_{\rm age}=0$. In the case of the W1.7 and W2.0 models, this can be regarded as the moment when the shock begins to accelerate preexisting electrons and become radio-bright. The figure shows that in the models with postshock TA (thick lines) the spectral steepening is significantly delayed relative to the models without TA (thin lines). Only the models with TA seem to produce $\alpha_{150}^{610}$ profiles broad enough to be compatible with the observed profile, which increases from $\alpha_{150}^{610}\approx 0.8$ to $\alpha_{150}^{610}\approx 2.0$ over $\sim 200$ kpc across the relic width. For the W1.7a and W2.0a models, the emissivity increases by an order of magnitude (a factor $8-12$) from upstream to downstream across the shock. Note that the subsequent, postshock emissivity decreases faster with time in the S3.6a model with only DSA injection from the background plasma, compared to the W1.7 and W2.0 models with the DSA reacceleration of the preexisting CRe. This is because for the particular injection model adopted here, the injection rate depends on $u_{\rm s}$ and $M_{\rm s}$, both of which decrease in time as the shock propagates.

3.2 Radio Surface Brightness Profiles

The radio surface brightness, $I_{\nu}$, is calculated by adopting the spherical wedge volume of radio-emitting electrons, specified with the two extension angles relative to the sky plane, $\psi_{1}$ and $\psi_{2}$, as shown in Figure 2 of Kang (2016a):

[TABLE]

where $R$ is the distance behind the projected shock edge in the plane of the sky (measured from the shock toward the cluster center), $r$ is the radial distance outward from the cluster center, and $h_{1}=r\sin{\psi_{1}}$ and $h_{2}=r\sin{\psi_{2}}$ are the path lengths along line of sight beyond and in front of the sky plane, respectively. (See Figure 1 of Kang (2015) for the geometrical meaning of $R$.)

Figure 3 shows the profiles of $I_{150}(R)$ and $\alpha_{150}^{610}(R)$, now calculated from $I_{150}(R)$ and $I_{610}(R)$, at the shock age of $t_{\rm age}=142-148$ Myr. The adopted values of $\psi_{1}$ and $\psi_{2}$ are given in the lower panels. In the weak-shock models with $M_{s,o}\approx 1.6-1.9$, a high-cutoff Lorentz factor, $\gamma_{e,c}\gtrsim 4\times 10^{4}$, is required to match $\alpha_{150}^{610}\approx 0.8$ at the shock position. From the geometric consideration only (that is, the line-of-sight length through the model relic), the first inflection point in the $I(R)$ profile occurs at $r_{\rm s}(1-\cos\psi_{1})\approx 38$ kpc for the shock radius $r_{\rm s}\approx 1.1$ Mpc and $\psi_{1}=15^{\circ}$, and the second inflection point occurs at $r_{\rm s}(1-\cos\psi_{2})\approx 87$ kpc for $\psi_{2}=23^{\circ}$. The third inflection point at $d\approx 150-160$ kpc occurs at the postshock advection length, $\sim u_{2}t_{\rm age}$, which corresponds to the width of the postshock spherical shell.

Note that the normalization factor for $I_{150}$ is arbitrary, but it is the same for all three models with $M_{s,i}=1.7$ (upper left panel) and for the three models with $M_{s,i}=2.0$ (upper middle panel). But note that for the S2.0d model $I_{150}$ is reduced by a factor of 0.6, compared to the other three models. So, for example, the relative ratio of $I_{150}$ between W1.7aN (without TA) and W1.7a (with TA) is meaningful. In the case of the S3.6 models (upper right panel), on the other hand, the normalization factor is the same for S3.6aN, S3.6a, and S3.6c (with only DSA injection of shock-generated suprathermal electrons), but a different factor is used for S3.6b (with preexisting, seed CRe) in order to plot the four models together in the same panel.

The effects of postshock TA can be seen clearly in the spectral steepening of $\alpha_{150}^{610}$ in the lower panels. As shown in Figure 4 below, for instance, the S3.6aN model (black) produces a “too-steep” spectral profile compared to observations, while the S3.6c model (green) without turbulence decay ($r_{dec}\rightarrow\infty$) produces a “too-flat” spectral profile.

To compare to the observed radio flux density distribution, $S_{\nu}$, the intensity, $I_{\nu}$, should be convolved with telescope beams. In Figure 4, a Gaussian smoothing with 23.5 kpc width (equivalent to $6.\arcsec 5$ at the distance of the Toothbrush relic) is applied to calculate $S_{\nu}(R)$, while the spectral index $\alpha_{150}^{610}$ is then calculated from $S_{150}(R)$ and $S_{610}(R)$. The observational data of van Weeren et al. (2016) are shown with magenta dots. The observed flux density at 150 MHz covering the region of $6.\arcsec 5\times 70\arcsec$ at $R\approx 50$ kpc behind the shock is $S_{150}\approx 0.20$ Jy. The required amount of preexisting CRe to match this flux level corresponds to $N\approx 0.4-0.5\%$ for the W1.7a,b and W2.0a,b models, and $N\approx 0.05\%$ for the S3.6b model. In the S3.6a model (without preexisting CRe), the corresponding flux density is $S_{150}\approx 0.004$ Jy, five times smaller than the observed value. Considering that the $\kappa=1.6$ distribution is already quite flat and so $\kappa$ index cannot be reduced further, it could be difficult to increase significantly the flux density $S_{150}$ in the S3.6a model. In that regard, the S3.6b model with preexisting CRe is favored over the S3.6a model. Note that the synchrotron intensity scales with $I_{150}\propto B_{2}^{(s-1)/2}$, while the downstream magnetic field strength in these models is chosen to be $B_{2,o}\approx 2.5-2.7~{}{\mu\rm G}$ (see Table 1) in order to maximize the downstream cooling length given in Equation (1).

In the upper panels of Figure 4, different normalization factors are adopted for each model to obtain the best match with the observed flux level of $S_{150}$ roughly at the peak values near $30-50$ kpc. The same relative normalization factors are scaled for the higher frequency and applied to $S_{610}$ in the middle panels. The observed profile of $S_{150}$ indicates that the region of the Toothbrush relic beyond $R>150$ kpc might be contaminated by a contribution from the radio halo.

We find that for the W1.7 and W2.0 models, a preexisting electron population with $s=4.4$ and $\gamma_{e,c}\gtrsim 4\times 10^{4}$ is necessary to reproduce the observed spectral steepening profile across the relic width. Moreover, the results demonstrate that the six models with TA (W1.7a,b, W2.0a,b, and S3.6a,b) can reproduce the observed profiles of $S_{\nu}(R)$ and $\alpha_{150}^{610}(R)$ reasonably well, while, as noted previously, none of the models without TA (black solid lines) can reproduce the profile of $\alpha_{150}^{610}$. However, it is also important to realize that the models should not produce “excess” TA. In particular, also as noted previously, the W2.0d model (green) with $\tau_{\rm acc,0}=50$ Myr and the S3.6c model (green) without turbulence decay produce “too-flat” profiles of $\alpha_{150}^{610}$.

At the time of observation, $M_{\rm s,o}\approx 3.03$, in the S3.6 models, so $\alpha_{\rm s}\approx 0.74$, which is slightly flatter than the observed index of 0.8 at the leading edge of the Toothbrush relic. This, we argue, is still consistent, because the observed radio flux profiles are blended by a finite telescope beam. We also considered a model (not shown) with $M_{s,i}=3.3$ with $M_{\rm s,o}\approx 2.85$, so at the time of observation, $q\approx 4.6$ ($\alpha_{\rm s}\approx 0.78$). That model, however, produces a spectral index profile across the relic a bit too steep to be compatible with the observed profile.

3.3 Volume-Integrated CRe and Radio Spectra

In the case of pure in situ injection without TA, the postshock momentum distribution function is basically the same as the DSA power-law spectrum given in Equation (9) except for the increasingly lower exponential cutoff due to postshock radiative cooling. So the volume-integrated CRe energy spectrum, $F_{\rm e}(p)=\int f_{\rm e}(r,p)dV$, is expected to have a broken power law form, whose slope increases from $q$ to $q+1$ at the break momentum, $p_{\rm br}/m_{\rm e}c\approx 10^{4}(t_{\rm age}/100{\rm Myr})^{-1}(5~{}{\mu\rm G})^{2}/(B_{2}^{2}+B_{\rm rad}^{2})$. In the lower right panel of Figure 1, for instance, we can see that the volume-integrated electron spectrum, $G_{\rm e}(p)=p^{4}F_{\rm e}(p)$, steepens gradually near $p/m_{\rm e}c\sim 3\times 10^{3}$ in the S3.6aN model (without TA, green long-dashed line).

Of course such a simple picture for the steepening does not apply to the W1.7 and W2.0 models with the DSA reacceleration of preexisting electrons, since the spectrum at the shock, $f_{\rm reacc}$, is a broken power-law that steepens from $p^{-s}$ to $p^{-q}$ above $p_{\rm e,c}$. In these models, $G_{\rm e}(p)$ depends on the assumed value of $\gamma_{e,c}$ (see the black, red, and blue lines in the lower left and lower middle panels of Figure 1) as well as $\tau_{\rm acc}$. The models without TA are also shown as green long-dashed lines for comparison.

In the S3.6a model in the lower right panel of Figure 1, the suprathermal $\kappa$-like population for $p\gtrsim p_{\rm inj}\approx 30\ m_{\rm e}c$ provides seed electrons for the in situ injection into DSA and subsequent TA in the postshock flow. In fact, this results in an excess, low-energy CRe population in the range $30\lesssim p/m_{\rm e}c\lesssim 300$ for the models, compared to the S3.6aN model, as shown in the figure. This low-energy component depends on the details of kinetic plasma processes operating near the shock, which are not yet fully understood, and would not contribute significantly to the observed radio emission in the range of $0.15-10$ GHz. For the postshock magnetic field strength, $B_{2}\approx 2.5~{}{\mu\rm G}$, electrons with $6.9\times 10^{3}\leq p/m_{\rm e}c\leq 5.6\times 10^{4}$ make the peak contribution in this observation frequency range.

From the spectral shape of $G_{\rm e}(p)$, we expect that the ensuing volume-integrated radio spectrum, $J_{\nu}=\int j_{\nu}(r)dV$, should steepen gradually toward high frequencies. Moreover, the form depends on $p_{e,c}$ and $\tau_{\rm acc}$ in the W1.7 and W2.0 models and on $p_{\rm br}$ and $\tau_{\rm acc}$ in the S3.6 models.

Figure 5 shows the volume-integrated radio spectrum, $\nu J_{\nu}$, for the W1.7a, W2.0a, S3.6a, and S3.6b models at three different shock ages to demonstrate how the spectrum evolves in time. For the models without TA (W1.7aN, W2.0aN, S3.6aN, and S3.6bN), the spectrum is shown only at the first epoch (the green long-dashed lines). In each panel, the normalization factor for the vertical scale is chosen so that the simulated curves match the observation data around 2 GHz. For the models without TA, the normalization factor is 1.6 times larger than for the corresponding models with TA. Note that the open squares (except at 4.85 and 8.35 GHz) are data for the B1 component of the Toothbrush relic in Table A1 of Stroe et al. (2016). Kierdorf et al. (2016) presented the sum of B1 + B2 + B3 flux at 4.8 and 8.35 GHz in their Table 5. Considering that the average ratio of the $B1/(B1+B2+B3)$ fluxes near 2 GHz is about 0.71 according to Tables 3 and A1 of Stroe et al. (2016), we lower the fluxes at 4.85 and 8.35 GHz in the Kierdorf’s data by the same factor. Basu et al. (2016) showed that the Sunyaev-Zel’dovich (SZ) decrement in the observed radio flux can be significant above 10 GHz for radio relics. We adopt their estimates for the SZ contamination factor for the Toothbrush relic given in their Table 1. Then the SZ correction factors, $F$, for the fluxes at 16 and 30 GHz are about 1.1 and 1.8, respectively. Two solid black filled circles correspond to the flux levels so-corrected at the two highest frequencies.

Although the models without TA do not reproduce the observed profile of $\alpha_{150}^{610}(R)$, as shown in Figure 4, the W1.7aN and W2.0aN models seem to fit the observed $J_{\nu}$ better than the W1.7a and W2.0a models. So this exercise teaches us that it is important to test any model against several different observed properties. Among the strong-shock models, S3.6a and S3.6b with TA seem to produce better fits to SZ-uncorrected $J_{\nu}$, while S3.6aN and S3.6bN without TA give the spectra more consistent with SZ-corrected $J_{\nu}$. In all models considered here, however, it seems challenging to explain the observed flux at 8.35 GHz.

In conclusion, adjustments of basic parameters can allow both of the weak-shock and strong-shock models to explain the observational data for the Toothbrush B1 component reasonably well. In the weak-shock scenario, as we argued in the Introduction, however, it would be challenging to fulfill the requirement for a homogeneous, flat-spectrum preexisting electron population over a region 400 kpc in length and $300$ kpc in width, which is needed to explain the observed uniformity in the spectral index along the length of the relic. If the preexisting electrons cool by radiative and collisional losses non-uniformly, or if the preshock CRe have a span in “ages”, both the cutoff energy and thus the spectral index at the relic edge would be expected to vary along the relic length.

To explore such effects, we compare in Figure 6 the weak-shock models allowing different cutoff energies, $10^{4}\leq\gamma_{e,c}\leq 10^{5}$. In order to reproduce the observed profiles of both $\alpha_{150}^{610}$ and $\nu J_{\nu}$, $\gamma_{e,c}\gtrsim 8\times 10^{4}$ is required for the W1.7 and W2.0 models. Considering that the cooling times for electrons with $\gamma_{e,c}=8\times 10^{4}$ in microgauss fields are only $\sim 13$ Myr, it would be very challenging to explain a constant $\gamma_{e,c}$ within the required preshock region.

In the right-hand panels of Figure 6, the W2.0d model (green long-dashed lines) shows that the “enhanced” TA with $\tau_{\rm acc,0}=50$ Myr would be too efficient to explain the observed profile of $\alpha_{150}^{610}(R)$. The model produces too many low-energy electrons with $\gamma_{\rm e}<10^{4}$, compared to high-energy electrons with $\gamma_{\rm e}\gtrsim 10^{4}$. This implies that the path to a model consistent with the observations cannot involve the adoption of smaller $\gamma_{e,c}$ combined with more rapid TA (smaller $\tau_{\rm acc}$).

Our results indicate that the strong-shock model with $M_{\rm s}\approx 3$ is favored. That could mean that the observed X-ray and radio Mach numbers represent different parts of a nonuniform shock surface (see the discussion in the Introduction). However, we should point out that the predicted $J_{\nu}$ values for the S3.6a and S3.6b models deviate from the observed curvature at 8.35 GHz (Figure 5). Finally, as noted earlier, in order to explain the rareness of detected radio relics in merging clusters, radio relics might be generated preferentially when shocks encounter regions of preexisting low-energy CRe (i.e., the S3.6b model).

4 SUMMARY

In this study, we reexamine the merger-driven shock model for radio relics, in which relativistic electrons are accelerated via DSA at the periphery of galaxy clusters. To that end, we perform time-dependent DSA simulations of one-dimensional, spherical shocks, and we compare the results with observed features of the Toothbrush relic reported by Stroe et al. (2016) and van Weeren et al. (2016). In addition to DSA, energy losses by Coulomb scattering, synchrotron emission, and iC scattering off the CMB radiation, and, significantly, TA by compressive MHD/plasma mode downstream of the shock are included in the simulations.

Considering apparently incompatible shock Mach numbers from X-ray ($M_{\rm X}\approx 1.2-1.5$) and radio ($M_{\rm rad}\approx 2.8$) observations of the Toothbrush relic, two possible scenarios are considered (see Table 1 for details): (1) weak-shock models in which a preexisting flat-spectrum electron population with high cutoff energy is accelerated by a weak shock with $M_{\rm s}\approx 1.6-1.9$, and (2) strong-shock models in which low-energy seed CRe, either shock-generated suprathermal electrons or preexisting soft-spectrum electrons, are accelerated by a strong shock with $M_{\rm s}\approx 3.0$.

The main results are summarized as follows:

1. In order to reproduce the broad profile of the spectral index behind the head (component B1) of the Toothbrush relic, TA with $\tau_{\rm acc}\approx 100$ Myr should be included to delay the spectral aging in the postshock region. This level of TA is strong but plausible in ICM postshock flows.

2. The strong-shock models with $M_{\rm s}\approx 3.0$, either with a $\kappa$-like distribution of suprathermal electrons (the S3.6a model) or with low-energy preexisting CRe with $p/m_{\rm e}c\lesssim 300$ (the S3.6b model), are more feasible than the weak-shock models. These models could explain the observed uniform spectral index profile along the relic edge over 400 kpc in relic length (component B1). Further, the S3.6b model may be preferred because (1) it can reproduce the observed flux density with a small fraction ($N\approx 0.05\%$) of preexisting CRe, and (2) it can explain the low occurrence ($\lesssim 10\%$) of giant radio relics among merging clusters, where otherwise “suitable” shocks are expected to be common. These low-energy fossil electrons could represent the leftovers either previously accelerated within the ICM by shock or turbulence or ejected from AGNs into the ICM, since their cooling times are long, $t_{\rm rad}>3.5$ Gyr with $B\sim 1~{}{\mu\rm G}$ for $\gamma_{e}<300$. The the S3.6a model, in which a $\kappa=1.6$ suprathermal distribution is adopted, the predicted flux density is about five times smaller than the observed level.

3. For the weak-shock models with $M_{\rm s}\approx 1.6-1.9$, a flat ($s\approx 4.4$) preexisting electron population with seemingly unrealistically high-energy cutoff ($\gamma_{e,c}\gtrsim 8\times 10^{4}$) is required to reproduce the observational data (the W1.7a and W2.0a models). It would be challenging to generate and maintain such a flat-spectrum preexisting population with a uniform value of $\gamma_{e,c}$ over the upstream region of 400 kpc in length and $300$ kpc in width, since the cooling time is short, $\tau_{\rm rad}\sim 10$ Myr for electrons with $\gamma_{\rm e}\sim 10^{5}$ in a $1~{}{\mu\rm G}$ level magnetic field.

H.K. was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2014R1A1A2057940). D.R. was supported by the National Research Foundation of Korea through grants 2014M1A7A1A03029872 and 2016R1A5A1013277. T.W.J. was supported by the US National Science Foundation through grant AST1211595. The authors thank R. J. van Weeren for providing the radio flux data for the Toothbrush relic published in van Weeren et al. (2016).