NONTHERMAL RADIATION FROM RELATIVISTIC ELECTRONS ACCELERATED AT SPHERICALLY EXPANDING SHOCKS

Abstract

We study the evolution of the energy spectrum of cosmic-ray electrons accelerated at spherically expanding shocks with low Mach numbers and the ensuing spectral signatures imprinted in radio synchrotron emission. Time-dependent simulations of diffusive shock acceleration (DSA) of electrons in the test-particle limit have been performed for spherical shocks with parameters relevant for typical shocks in the intracluster medium. The electron and radiation spectra at the shock location can be described properly by the test-particle DSA predictions with instantaneous shock parameters. However, the volume integrated spectra of both electrons and radiation deviate significantly from the test-particle power-laws, because the shock compression ratio and the flux of injected electrons at the shock gradually decrease as the shock slows down in time. So one needs to be cautious about interpreting observed radio spectra of evolving shocks based on simple DSA models in the test-particle regime.

1. INTRODUCTION

X-ray emitting baryonic gas in galaxy clusters is thought to be heated via collisioness shocks that are induced during the formation of the large scale structure in the Universe (e.g., Ryu et al. 2003; Kang et al. 2007; Vazza et al. 2009; Skillman et al. 2011). The gas infalling toward nonlinear structures go through strong accretion shocks first, and later may encounter much weaker shocks produced by mergers and chaotic flow motions inside the hot intracluster medium (ICM) (Ryu et al. 2003). In fact, about a dozen of these weak ICM shocks have been detected through X-ray observations as bow shocks associated with merging activities (e.g., the so-called bullet cluster 1E 065756, Markevitch et al. 2002; Akamatsu & Kawahara 2013; Ogrean & Brüggen 2013).

In addition, quiet a few ICM shocks have been identified as “radio relic shocks” in the outskirts of galaxy clusters through radio synchrotron radiation from cosmic-ray (CR) electrons that are thought to be accelerated via diffusive shock acceleration (DSA) (van Weeren et al. 2010, 2011; Nuza et al. 2012). Radio relics are interpreted as synchrotron emitting structures that contain relativistic electrons with Lorentz factor γe ~ 104 gyrating in mean magnetic fields B ~ 5−7 μG in the ICM (Kang et al. 2012; Feretti et al. 2012; Pinzke et al. 2013; Brunetti & Jones 2014). The electrons are expected to be accelerated at weak shocks with the sonic Mach number, Ms ~ 2 − 4, inferred from the spectral index αinj ＞ 0.6 of synchrotron emission immediately behind the shock. Note that the spectral index of the volume-integrated synchrotron radiation is normally Aν = αinj + 0.5, since electrons cool via synchrotron and inverse-Compton (iC) losses behind the shock (e.g., Kang 2011).

It is well accepted that the DSA theory can explain how nonthermal particles are produced through their interactions with MHD waves in the converging flows across astrophysical shocks (Bell 1978; Drury 1983; Malkov & Drury 2001). In the test particle limit where the CR pressure is dynamically insignificant, it predicts that the CR energy spectrum at the shock position has a power-law energy spectrum, N(E) ∝ E−s, where s = (σ + 2)/(σ − 1) and σ = ρ2/ρ1 is the shock compression ratio. (Hereafter, we use the subscripts ‘1’ and ‘2’ to denote the conditions upstream and downstream of shock, respectively.) Hence, the DSA model can provide a simple and natural explanation for the spatially resolved synchrotron radiation spectrum at the shock location, jν(xs) ∝ ν−αinj , where αinj = (s − 1)/2.

On the other hand, pre-acceleration of thermal electrons to suprathermal energies and the subsequent injection into the DSA process still remains very uncertain, especially at low Mach number shocks in high beta (β = Pg/PB) plasmas (e.g., Kang et al. 2014). Recently, there have been some serious efforts to understand complex plasma instabilities and wave-particle interactions operating at nonrelativistic shocks with the parameters relevant for weak cluster shocks, through Particle-in-Cell (PIC) and hybrid plasma simulations (e.g., Riquelme & Spitkovsky 2011; Caprioli & Spitkovsky 2014a,b; Guo et al. 2014). They have demonstrated that self-excitation of MHD/plasma waves via various instabilities, such as Bell’s resonant and nonresonant instabilities and firehose instability, are crucial in the injection process. Moreover, it has been shown that protons can be injected efficiently at quasi-parallel shocks, while the electron injection occurs preferentially at quasi-perpendicular shocks. Here, we adopt a phenomenological injection model in which suprathermal electrons with momentum p > pinj are allowed to cross the shock transition and participate in the Fermi 1st-order acceleration, where the injection momentum, pinj, is several times the peak momentum of the postshock thermal protons (Kang et al. 2002). This subject will be discussed in detail in Section 2.2.

The morphology of observed radio relics varies rather widely from diffuse patches to well-defined, thin elongated structures. For example, the so-called Sausage relic in CIZA J2242.8+5301 is located at a distance of ~ 1.5 Mpc from the cluster center and has an arc-like shape of ~ 55 kpc in width and ~ 1 − 2 Mpc in length (van Weeren et al. 2010). So it would be natural to assume, at least for cases similar to the Sausage relic, that radio relics consist of some portions of a spherical shell of radiating electrons projected onto the sky plane. The physical width of this shell is mainly determined by the cooling length of electrons, ∆l ~ u2trad(γe), where u2 is the postshock flow speed and trad is the electron cooling time. But the projection of a partial shell could involve some geometrical complexities such as the extent of the partial shell and the viewing angle (e.g., van Weeren et al. 2010, 2011; Kang et al. 2012).

For several observed radio relics the projected spatial profiles of the radio flux, Sν, and the spectral index, αν, have been explained by the electron energy spectrum Ne(r, γe) cooling radiatively behind onedimensional planar shocks (e.g., van Weeren et al. 2010, 2011; Kang et al. 2012). This may be justified, since the width (~ 50 − 100 kpc) of those relics is much smaller than their length (~ 1 − 2 Mpc). Note that a simple projection of a spherical shell on the sky plane was attempted in these previous papers, although the DSA simulation results for a planar shock was adopted for Ne(r, γe).

In this work we consider electron acceleration at spherical shocks with Ms ~ 2.5−4.5, which expand into a hot uniform ICM. A self-similar solution for SedovTaylor blast wave is adopted as the initial states of the postshock gasdynamic quantities (i.e., ρ, u, P) at the beginning of the simulations. Later on the shock structure deviates from the Sedov-Taylor similarity solution, because the sonic Mach numbers of our model shocks are not large enough. Nevertheless, the shock slows down approximately as us(t) ∝ t−3/5. This model may represent a spherically decelerating shock, but it is not meant to mimic realistic cluster shocks. Note that the flow speed, u(r, t), is the only dynamical information that is fed into the diffusion convection equation for the evolution of the electron distribution function, fe(r, p, t). As the shock slows down and weakens, the shock compression ratio and the injection flux of CR electrons should decrease in time. Consequently, the amplitude of electron spectrum at the injection momentum decreases and the electron energy spectrum steepens gradually.

Considering that the synchrotron/iC cooling time scale (trad ~ 108 yr) for the electrons of our interest (γe ~ 104 ) is shorter than the dynamical time scale of typical cluster shocks (tdyn ~ 109 yr), the electron energy spectrum at the shock location is expected to follow the test-particle DSA predictions with the instantaneous shock parameters. But the spatial profile of the electron distribution function, fe(r, p), and the volume-integrated distribution function, Fe(p) = 4πfe(r, p) r2dr, would be affected by the time-dependent evolution of the spherically decelerating shock. Our primary objective here is to study any signatures of such time-dependence through DSA simulations of a heuristic example (i.e., a spherical blast wave).

In the next section we describe the numerical method and some basic physics of the DSA theory. The simulation results of a planar shock and a spherical shock will be compared, and spherical shocks with different magnetic field models will be discussed in Section 3. A brief summary will be given in Section 4.

 

2. DSA NUMERICAL SIMULATIONS

2.1. 1D Spherical CRASH Code

We consider DSA of CR electrons at nonrelativistic, gasdynamical shocks in one-dimensional (1D) spherical geometry. So we solve the following timedependent diffusion-convection equation for the pitchangle-averaged phase space distribution function for CR electrons, ge(r, p, t) = fe(r, p, t) p4:

where u(r, t) is the flow velocity, y = ln(p/mec), me is the electron mass, c is the speed of light, and D(r, p) is the spatial diffusion coefficient (Skilling 1975). For D(r, p), we adopt a Bohm-like diffusion coefficient with a weaker non-relativistic momentum dependence

The cooling term b(p) = −dp/dt accounts for electron synchrotron and iC losses:

in cgs units, where e is the electron charge. The ‘effective’ magnetic field strength takes account for radiative losses due to both synchrotron and iC processes, where Brad = 3.24 μG(1+z)2 corresponds to the cosmic background radiation (CBR) at redshift z (Schlickeiser 2002). In this study, we set z = 0.2 as a reference epoch and so Brad = 4.7 μG. Then the cooling time scale for electrons is given as

For typical ICM conditions relevant for radio relics (i.e., nH ~ 10−4 cm−3 and B ~ 1 μG), the Coulomb loss dominates over the other losses for relativistic electrons with γe 100, while the synchrotron and iC losses are the main cooling processes for γe 100 (Sarazin 1999). Nonthermal bremsstrahlung loss is much less important and can be ignored for relativistic electrons, since its cooling time scale is longer than the Hubble time (Petrosian 2001).

Since we focus on the synchrotron emission from CR electrons here, we do not need to consider the acceleration of CR protons. Moreover, the CR proton pressure is expected to be dynamically insignificant at weak shocks with Ms 5, so we can just follow the hydrodynamic evolution of the shock without following DSA of CR protons. This significantly alleviates the requirements for computational resources, since much wider ranges of particle energy and diffusion length/time scales should be included, if both proton and electron populations were to be calculated. For instance, the maximum energy of electrons are limited to ~ 100 TeV for the problem considered here, while that of protons can go up much higher energies by several orders of magnitude.

At weak shocks in the test-particle limit, the CR feedback becomes negligible and the background flow, u(r, t), is governed by the usual gasdynamic conservation equations in 1D spherical coordinates (Kang & Jones 2006):

where P is the gas pressure, eg = P/[ρ(γg − 1)] + u2/2 is the total energy of the gas per unit mass and the rest of variables have their usual meanings. The gas adiabatic index is assumed to be γg = 5/3, since the background gas flow is nonrelativistic.

In order to optimize the shock tracking scheme, a comoving frame that expands with the instantaneous shock speed is adopted. In such a frame, a spherically expanding shock can be made to be stationary by adopting comoving variables which factor out a uniform expansion. The details of the CRASH (CosmicRay Amr SHock) code in a comoving spherical grid can be found in Kang & Jones (2006).

2.2. Electron Pre-heating and Injection

Injection of suprathermal particles into the Fermi process remains one of the outstanding problems in the DSA theory. In the so-called thermal leakage injection model, suprathermal particles with p pinj ~ 3pth,p are expected to re-cross the shock upstream and participate in the DSA process, because the shock thickness is of the order of the gyroradius of postshock thermal protons (e.g., Kang et al. 2002). Here, is the most probable momentum of thermal protons with postshock temperature T2 and kB is the Boltzmann constant. On the other hand, recent hybrid and PIC simulations of collisioness shocks have revealed a somewhat different picture, in which incoming protons and electrons are specularly reflected at the shock surface and scattered by plasma waves excited in the foreshock region, leading to suprathermal populations (e.g., Guo & Giacalone 2013; Guo et al. 2014; Caprioli et al. 2015). The trajectories of these suprathermal particles are non-diffusive and confined to the region mostly upstream of the shock transition.

If one assumes that downstream electrons and protons have Maxwellian distributions defined with the same kinetic temperature, Te ≈ Tp, then thermal electrons have much higher speeds but much smaller momenta compared to thermal protons (i.e., pth,e = (me/mp) 1/2pth,p). So electrons must be pre-accelerated from their thermal momenta to pinj ~ 3pth,p ~ 130pth,e, before they can begin to take part in the DSA process. In fact, pre-heating of electrons above the thermal distribution has been widely observed in space and laboratory plasmas (e.g., Vasyliunas 1968). Such suprathermal non-Maxwellian distributions can be described by the combination of a Maxwellian-like core and a power-law tail, which is known as the κ-distribution:

where is the thermal peak momentum, the mass of the particle is m = me for electrons and m = mp for protons, and Γ(x) is the Gamma function (e.g., Pierrard & Lazar 2010). In the limit of large κ, it asymptotes to the Maxwellian distribution.

The power-law index of the κ-distributions for electrons and protons, κe and κp, respectively, should depend on plasma and shock parameters, i.e., β, Ms, MA = us/vA (where is the Alfvén speed), and the shock obliquity angle, ΘBn. For example, the electron distributions measured near interplanetary shocks can be fitted with the κ-distributions with κe ~ 2 − 5, while the proton distributions prefer a somewhat larger κp 10 (e.g., Pierrard & Lazar 2010).

As an illustration, Figure 1 compares the electron κ- distributions, fe,κ(p) with κe = 1.6 − 3.0 for p ≤ pinj and the electron power-law distribution, fe,κ(pinj) · (p/pinj) −q for p > pinj. Here the shock parameters adopted are us = 3.2 × 103 km s−1, Ms = 3, and the DSA power-law slope, q = 4.5. The Maxwellian distributions for electrons and protons are also plotted for comparison. These κ-distributions may represent the pre-heated electron populations, while the power-law distributions may represent the shock accelerated electron populations. We note that pre-heated electrons are assumed to be injected at p = pinj with the amplitude, fe,κ(pinj), which is much larger than the corresponding value for the electron Maxwellian distribution, fe,M(pinj). Here we adopt pinj ≈ 5.34mpu2, which becomes pinj ≈ 3.3pth,p for Ms = 3. The filled circle on the proton Maxwellian distribution marks the amplitude of proton distribution, fp,M(pinj), at the injection momentum, which is larger than fe,κ(pinj) by several orders of magnitude.

Figure 1.Schematic model demonstrating the pre-heated suprathermal distributions and the power-law distributions of electrons accelerated at a weak shock of us = 3.2 × 103 km s−1 and Ms = 3. The Maxwellian distribution (black solid line) and κ-distributions with κe = 1.6 (red dotted), 2.0 (blue dashed), 2.5 (green dot-dashed), and 3.0 (magenta dot-long dashed) are shown for electrons. For protons only the Maxwellian distribution (black dot-long dashed line) is shown. The vertical line demarcates the injection momentum, so fe(p) = fe,κ(p) for p < pinj and fe(p) = fe,κ(pinj) · (p/pinj)−q for p ≥ pinj. The filled circle marks the amplitude of fp,M (pinj) at the injection boundary for the proton Maxwellian distribution.

For the case shown in Figure 1, the ratio of CR electron to proton numbers can be approximated by Ke/p ≈ fe,κ(pinj)/fp,M(pinj), which becomes Ke/p ~ 1/300 for κe = 1.6. So Figure 1 demonstrates that the ratio becomes Ke/p ~ 10−3 − 10−2, depending on the value of pinj, if the electrons are pre-accelerated to κ-distributions with κe 2.

Since the main goal of this study is to explore how the electron energy spectrum and the synchrotron emission are affected by the dynamical evolution of spherically expanding shocks (i.e., us(t) and rs(t)), the specific value of the amplitude fe,κ(pinj) for a given set of models parameters (i.e., κe and pinj/pth,p) is not important. So we simply assume that the injection momentum decreases with the shock speed as pinj(t) ≈ 5.34mpus(t)/σ, and the electrons are injected with the value Ke/p = 1 in the simulations described below. Note that we do not attempt to compare the theoretically estimated radio flux, Sν, with the observed radio flux of any specific radio relics.

2.3. Shock Accelerated CR Electron Spectrum

In the test-particle regime of DSA, the distribution of CR distribution function at the shock position can be approximated, once it reaches equilibrium, by a power-law spectrum with super-exponential cutoff,

where the power-law slope is q = 3σ/(σ − 1) and finj is the amplitude at p = pinj (Kang 2011). The cutoff momentum can be derived from the equilibrium condition that the DSA momentum gains per cycle are equal to the synchrotron/iC losses per cycle (Kang 2011):

The corresponding Lorentz factor for typical cluster shock parameters is then

Postshock electrons cool radiatively while advecting downstream, so the cutoff of the electron momentum spectrum decreases as one moves away from the shock. As a result, the volume integrated electron energy spectrum, Fe(p), becomes steeper than the power-law in Equation (9) by one power of the momentum above the ‘break momentum’ (p > pe,br). At the shock age t, the ‘break Lorentz factor’ can be estimated from the condition tage = trad = p/b(p):

The synchrotron emission from mono-energetic electrons with γe peaks around the frequency νpeak ≈ 0.3(3eB/4πmec), so it will be useful to have the following relation:

Then the synchrotron radiation spectrum emitted by the electron population with the power-law given in Equation (9) also has a power-law form of jν ∝ ν−αinj with .

For diffuse radio structures such as radio halos and radio relics with B ~ 5 μG and the radio flux density, Sν ~ 5 mJy, the optical depth for synchrotron emission is small for the frequencies considered here (ν > 10 MHz) (Jones et al. 1974; Lang 1999). So the synchrotron self-absorption can be ignored in the calculation of synchrotron radiation spectrum of the model shocks.

 

3. DSA SIMULATION RESULTS

3.1. Planar versus Spherical Shocks

In this section we compare the electron acceleration at a spherical shock with that at a planar shock with similar properties. First, the planar shock is characterized by the following parameters: the shock speed, us = 3.5 × 103 km s−1, the sonic Mach number, Ms = 4.5, the preshock density, nH,1 = 10−3 cm−3, and the preshock and postshock magnetic field strength, B1 = 2 μG (for x > xs) and B2 = 7 μG (for x < xs).

In the upper panel of Figure 2, we show the spatial profiles of the electron distribution function, ge(x, γe) for different energies, log γe = 3.3, 3.7, 4.0, 4.3, 4.5, and 4.8. Note that the shock location is shifted slightly from left to right at later time epochs, tage = (1.8, 4.4, 8.0) × 107 yr. The downstream width of ge, L(γe), for low energy electrons with log γe ≲ 3.3 increases simply as ladv = (us/σ)tage ≈ 100kpc(tage/108 yr). For high energy electrons, on the other hand, L(γe) decreases with γe, since the radiative cooling time scales inversely with electron energy as in Equation (4). As can be seen in the figure, for tage > 4.4 × 107 yr, for instance, the width L(γe) for γe ≥ 104 becomes constant in time (tage > trad) as L(γe) ≈ (us/σ)trad(γe) (see blue dashed, green dot-dashed, magenta dot-long dashed, and cyan dash-long dashed lines).

Figure 2.Time evolution of the planar shock of us = 3.5 × 103 km s−1, Ms = 4.5, and B2 = 7 μG is shown for tage = (1.8, 4.4, 8.0)×107 yr (from left to right). For clarity the shock location is shifted to the right by arbitrary amount at different time epochs. Top: Spatial distributions of the electron distribution function, ge(x, γe), for log(γe) = 3.3 (black solid lines), 3.7 (red dotted), 4.0 (blue dashed), 4.3 (green dot-dashed), 4.5 (magenta dot-long dashed), and 4.8 (cyan dash-long dashed). Middle: Spatial distributions of the synchrotron emissivity, jν(x, ν) for log(ν) = 7.6 (black solid lines), 8.5 (red dotted), 9.0 (blue dashed), 9.6 (green dot-dashed), 10. (magenta dot-long dashed), and 10.7 (cyan dash-long dashed). Bottom: Spectral index αlog ν1−log ν2, where (log ν1,log ν2) is (7.6,8.5) (black solid line), (8.5,9.0) (red dotted), (9.0,9.6) (blue dashed), (9.6,10.) (green dot-dashed), and (10.,10.7) (magenta dot-long dashed).

From Equation (13), the electrons with γe are expected to have the peak of synchrotron emission at νpeak ≈ 1GHz(γe/104)2. The middle panel of Figure 2 shows the synchrotron emissivity, jν(x, ν), for the peak frequency that corresponds to the Lorentz factor chosen in the upper panel (i.e., log ν = 7.6, 8.5, 9.0, 9.6, 10., and 10.7), at the same epochs as in the upper panel. Note that the spatial distribution of synchrotron emission at νpeak is much broader than that of electrons with the corresponding γe. This is because more abundant lower energy electrons also contribute significantly to the synchrotron emissivity at νpeak. For example, the spatial profile of jν(1GHz) (blue dashed lines) is much broader and decreases much gradually, compared to that of ge of the corresponding γe = 104 (also blue dashed lines) with L(γe = 104) ≈ 40 kpc. So the downstream width of synchrotron emitting volume would be greater than that of the spatial distribution of ge(γe) for the value of γe related by Equation (13), that is, L(νpeak) > L(γe), if the width is constrained by the cooling time trad. Note that ge and jν are plotted in arbitrary units in Figure 2.

The bottom panel of Figure 2 shows the spectral index which is defined as

estimated between ν1 and ν2. The peak frequencies chosen in the middle panel are adopted for ν1 and ν2, so α7.6−8.5 (black solid line), α8.5−9.0 (red dotted), α9.0−9.6 (blue dashed), α9.6−10 (green dot-dashed), and α10−10.7 (magenta dot-long dashed) are plotted. The characteristic scale over which αlog ν1−log ν2 increases is governed mainly by the length scale of jν at the higher frequency ν2. For example, the spectral index α9.0−9.6(blue dashed line in the bottom panel) behaves similarly to the emission jν at log ν = 9.6 (green long dashed line in the middle panel).

For the spherical shock, we adopt a Sedov-Taylor similarity solution with the shock speed, us,i = uo(ti/to)−3/5, as the initial set up. Here ti is the initial time from which electrons are injected at the shock and uo = 2.3 × 103 km s−1 and to = 1.73 × 108 yr. For ti = 0.5to, the initial shock speed, us,i ≈ 3.5 × 103 km s−1 and Ms ≈ 4.5, which are similar to the parameters for the planar shock shown in Figure 2. Since the sonic Mach number is not large and the shock compression ratio is σ = ρ2/ρ1 ≈ 3.5, the shock evolution behaves slightly differently from the true Sedov-Taylor solution. So at later time the shock speed evolves only approximately as us(t) ≈ uo(t/to)−3/5. Again note that only the structure and evolution of the flow velocity u(r, t), is fed into the diffusion convection Equation (1) and the other gasdynamical variables do not influence the DSA process. For the comparison with the planar shock, the magnetic field strengths are assumed to be constant in both space and time with B1 = 2 µG and B2 = 7 µG. In the next section, we will consider several models with variable magnetic field strength, B(r, t), behind the shock.

The characteristics of the spherical shock model shown in Figure 3 can be compared directly with those of the planar shock shown in Figure 2. Here the results are shown at the time elapsed from the initial injection, tage ≡ (t−ti) = (1.7, 4.3, 8.7)×107 yr, during which the shock speed decreases from 3.5×103 km s−1 (Ms = 4.5) to 1.5 × 103 km s−1 (Ms = 3.1). The obvious feature, which is different from the planar shock case, is that the injected particle flux decreases as the spherical shock slows down in time. As a result, the amplitude of ge(r, γe) increases downstream away from the shock at a given time, especially in the case of low energy electrons, which simply advect downstream without much cooling. So low frequency (ν 360 MHz) radio emission may increase towards downstream as well. For high energy electrons, on the other hand, the decline of ge(r, γe) due to cooling dominates over the effect of decreasing injection flux. Thus radio emission at higher frequencies (ν 1 GHz) should decrease downstream behind the shock. Again note that ge and jν are plotted in arbitrary units in Figure 3.

Figure 3.Time evolution of the spherical shock of us = 2.3 × 103 km s−1 (t/to)−3/5, Ms = 3.0(t/to) −3/5, and B2 = 7 µG (where to = 1.73 × 108 yr) is shown for tage = (1.7, 4.3, 8.7) × 107 yr (from left to right). Electrons are injected at the shock from ti = 0.5to when us,i = 3.5 × 103 km s−1. Top: Spatial distributions of ge(r, γe). Middle: Spatial distributions of jν(r, ν). Bottom: Spectral index αlog ν1−log ν2. The line types are the same as in Figure 2.

Similar to the planar shock model, the spatial distribution of synchrotron emission at νpeak is much broader than that of ge(γe) for the corresponding γe. For example, again the width of ge, L(γe = 104) ≈ 40 kpc (blue dashed line in the top panel), while the width of jν, L(ν ≈ 1GHz) ≈ 100 kpc (blue dashed line in the middle panel). But the spectral index α9.0−9.6 between 1 GHz and 4 GHz (blue dashed line in the bottom panel) follows the spatial profile of the emission jν at 4 GHz (green dot-dashed line in the middle panel). Thus the main difference between a planar shock and a decelerating spherical shock is the decline of the injected particle flux in time, leading to possible increase of low frequency radio emission downstream away from the shock.

Next we attempt to find any spectral signatures imprinted in the synchrotron emission spectrum in the case of the spherical shock case. In the left panels of Figure 4, we show, at the three epochs, the electron spectrum at the shock, ge(xs, p), and the volume integrated electron spectrum Ge(p) = p4Fe(p), where Fe(p) = fe(p, x)dx, and their slopes, q = −d ln fe(xs)/d ln p, and Q = −d ln Fe/d ln p. The slope from the simulation results is consistent with the test-particle slope q = 4.2 for σ = 3.5, and the slope for the volume integrated spectrum, Q = q+1, above the break momentum pe,br. Note that γe = p/mec for relativistic energies. As can be seen in Figure 4, the electron energy spectrum at the shock has reached steady state at these ages. The volume integrated spectrum shows the gradual decline of the break Lorentz factor, γe,br with the shock age.

Figure 4.Electron (left) and synchrotron (right) spectral properties for the planar shock shown in Figure 2 are plotted for tage = (1.8, 4.4, 8.0) × 107 yr (solid, dotted, and dashed lines, respectively). Upper left: electron distribution function at the shock position, ge(xs, p) = p4fe(xs, p), and integrated over the simulation volume, Ge(p) = p4Fe(p). Lower left: electron distribution slopes, q = −d ln fe(xs)/d ln p and Q = −d ln Fe/d ln p. Upper right: synchrotron radiation spectrum at the shock position, νjν(xs), and integrated over the simulation volume, νJν. Lower right: synchrotron spectral indexes, αν = −d ln jν(xs)/d ln ν and Aν = −d ln Jν/d ln ν.

The same kinds of quantities for the spherical shock are shown in the left panels of Figure 5. We can see the time evolution of the electron spectrum as the shock slows down in time. Since the shock Mach number decreases, the slope q for ge(rs, p) increases in time, leading to a nonlinear change of the volume integrated spectrum, Ge(p). So the slope Q no longer follows the simple q+1 steepening above pe,br. In addition to usual radiative cooling, both the reduced injected particle flux and steepening of ge(rs, p) due to the slowing-down and weakening of the shock affect the spectral shape of Ge(p).

Figure 5.Same as Figure 4 except that the results for the spherical shock shown in Figure 3 are plotted for tage = (1.7, 4.3, 8.7) × 107 yr (solid, dotted, and dashed lines, respectively).

In addition, we show the synchrotron emission spectrum at the shock, jν(xs) in Figure 4 and jν(rs) in Figure 5, and the volume integrated emission spectrum, Jν = jν (x)dx in Figure 4 and Jν = jν(r)4πr2dr in Figure 5. Their respective spectral indexes, αν = −d ln jν/d ln ν and Aν = −d ln Jν/d ln ν are shown in the lower right panels of Figures 4 and 5. For the planar shock with σ = 3.5, the test-particle spectral indexes are αν = 0.6 and Aν = 1.1 for ν > νbr, as can be seen in Figure 4. We note that the transition of the volume-integrated spectral index from Aν = αν to Aν = αν + 0.5 is more gradual, compared to the transition of the particle slope from Q = q to Q = q + 1. For the shock parameters chosen here, at the shock age of ~ 108 yr this transition occurs over the broad frequency range, 0.01 GHz ν 10 GHz.

For the spherical shock shown in Figure 5, on the other hand, the volume integrated spectrum, Jν, develops a concave curvature and its spectral index Aν deviates from the simple relation of Aν = αν + 0.5. Again, at the shock age of ∼ 108 yr the nonlinear curvature occurs over the broad frequency range, 0.01 GHz ν 10 GHz. Thus the simple relation, Ainteg ≈ αinj + 0.5, which is often invoked for the DSA origin of synchrotron emission spectra, should be applied with some caution in the case of evolving shocks. Again in Figures 4 and 5, ge, Ge, jν and Jν all are plotted in arbitrary units.

3.2. Spherical Shocks with Different Magnetic Field Models

Self-excitation of MHD waves and amplification of magnetic fields via plasma instabilities are the integral parts of DSA at strong collisionless shocks. CRs streaming upstream in the shock precursor excite resonant Alfvén waves with wavelengths comparable with the CR proton gyroradii, and turbulent magnetic fields can be amplified into the nonlinear regime (Bell 1978; Lucek & Bell 2000). In addition, the nonresonant fastgrowing instability driven by the CR current can amplify the magnetic field by orders of magnitude (Bell 2004; Riquelme & Spitkovsky 2009). It was shown that, in the case of strong shocks (Ms > 10), magnetic ampli- fication (MFA) via streaming instabilities and Alfvénic drift can lead to significant nonlinear feedback on the shock structure and the predicted CR spectrum (e.g., Caprioli 2012; Kang 2012).

In a recent study of hybrid plasma simulations, Caprioli & Spitkovsky (2014b) has shown that the MFA factor increases with the Alfvenic Mach number as 〈δB/B〉2 ~ , where δB is the turbulent magnetic fields perpendicular to the mean back-ground magnetic fields. For the sonic Mach number in the range of 3 Ms 5, the CR proton acceleration efficiency varies as ~ 0.01 − 0.05 (Kang & Ryu 2013; Caprioli & Spitkovsky 2014a). For typical cluster shocks with MA ~ 10, the MFA factor due to the streaming stabilities is expected be rather small, 〈δB/B〉2 ~ 0.3 − 1.5. How magnetic fields may be amplified at weak shocks in high beta ICM plasmas has not yet been fully understood. As in the previous study of Kang et al. (2012), we adopt the postshock magnetic field strength B2 ~ 7 µG in order to produce the observed width of the Sausage relic (∆l ~ 50 kpc), while the preshock magnetic field strength is chosen to be B1 ~ 2 − 3 µG. Note that B1 cannot be constrained directly from observations, since we do not detect any radio emission upstream of radio relic shocks. Thus B1 can be adjusted so that B2 becomes about 7 µG after considering MFA via plasma instabilities or the compression of the perpendicular components of magnetic fields across the shock. In any case, the postshock electron energy spectrum and the synchrotron radiation spectrum do not depend sensitively on the value of B1 as long as it is a few microgauss.

Here we consider the following three magnetic field models:

MF1: B1 = 2 µG & B2 = 7 µG. MF2: B1 = 3 µG, B2 = & B(r) = B2 · (P(r)/P2)1/2 for r < rs. MF3: B1 = 3 µG, B2 = & B(r) = B2 · (ρ(r)/ρ2) for r < rs.

In MF2 and MF3 models, the postshock magnetic field energy is assumed to be dissipated, and so B(r) decreases with the gas pressure or density behind the shock. As shown in the bottom panels of Figure 6, the decay length scale is about 100-150 kpc in MF2 and MF3 models.

Figure 6.Spherical shock models with three different magnetic field profiles: MF1, MF2, and MF3. The results are shown for tage = 1.3×108 yr, when us = 2.5×103 km s−1 and Ms = 2.4. Spatial distributions of the electron distribution function, ge(r, γe), the synchrotron emissivity, jν(r, ν), the spectral index, αlog ν1−log ν2, and magnetic field strength, B(r) (from top to bottom panels). The line types for ge, jν, and αlog ν1−log ν2 are the same as in Figures 2 and 3. Note that ge and jν are plotted in arbitrary units.

As in the previous section, we consider a spherical blast-wave, us,i ≈ uo(ti/to)−3/5, where uo = 3.0 × 103 km s−1 and to = 1.3 × 108 yr (slightly faster than the shock adopted in Figure 3). The calculations were started at ti = 0.5to with us,i ≈ 4.5 × 103 km s−1 and Ms ≈ 4.3 and ended at tf = 1.5to with us,f ≈ 2.5 × 103 km s−1 and Ms ≈ 2.4. The results of these three models at the final time (tage = t − ti ≈ 1.3 × 108 yr) are compared in Figure 6. As in the spherical shock model shown in Figure 3, the spatial profile of ge for low energy electrons increases downstream away from the shock. The distributions of ge(γe ≈ 5 × 103) (red dotted lines) show some differences among the three models due to different synchrotron cooling rates. The spatial profiles of jν at low frequencies (e.g., 360 MHz in red dotted lines) show rather weak dependences on the B(r) model, while those at higher frequencies (e.g., 1 GHz in blue dashed lines) behave similarly in all three cases. Moreover, the steepening patterns of αν behind the shock show no significant variations among the three models.

Thus it may be difficult to differentiate the postshock magnetic field profile from the spatial variation of spatially resolved radio flux, Sν, or spectral index, αν. In other words, signatures imprinted on synchrotron emission due to different B(r) behind the shock would be rather subtle to detect. There are two reasons that may explain this prediction. Because the iC cooling with Brad = 4.7 µG provides the baseline for the cooling rate, b(p), the difference in the synchrotron cooling due to different B(r) has relatively minor effects. Also any sharp features in the distribution of fe(r, p) would be smeared out in the distribution jν(r), since the synchrotron emission at a given frequency comes from electrons with a somewhat broad range of γe.

Finally, in Figure 7 we compare the synchrotron radiation spectrum at the shock, jν(rs), and the volume integrated radiation spectrum, Jν, for the three magnetic field models at tage = (0.33, 0.67, 1.3) × 108 yr. Departures from the predictions for the test-particle planar shock are most severe in MF1 model, while it becomes relatively milder in the two models with decaying B(r) behind the shock. As expected, any variations in the spatial profiles of fe(r, p) and B(r) are averaged out in the volume integrated quantities such as Fe(p) and Jν.

Figure 7.Spherical shock models with three different magnetic field profiles shown in Figure 6. Top: synchrotron spectral distribution at the shock position, νjν(rs), and integrated over the simulation volume, νJν. Bottom: synchrotron spectral index, αν for jν(rs) and Aν for Jν . The results are shown for tage = (0.33, 0.67, 1.3) × 108 yr (solid, dotted, dashed lines, respectively). Note that jν and Jν are plotted in arbitrary units.

In conclusion, any features in the spectral curvature, dAν/dν, or the steepening of spectral index, dαν/dr in the postshock region, should be affected in principle by both the time-dependent evolution of the shock and the spatial variation of postshock magnetic fields. But the effect of decaying magnetic fields may be somewhat subtle to detect due to the dominance of iC scattering off the CBR.

 

4. SUMMARY

We have performed time-dependent DSA simulations for cosmic-ray (CR) electrons at a planar shock with a constant speed and decelerating spherical shocks similar to Sedov-Taylor blast waves. These shocks were chosen to have parameters relevant for weak shocks typically found in the outskirts of galaxy clusters: us ≈ (1.5−4.5)×103 km s−1, Ms ≈ 2.5−4.5, B1 ≈ 2−3 µG, and B2 ≈ 7 µG. For spherical shocks, we used a onedimensional spherical version of the CRASH code, in which a co-expanding spherical grid was adopted in order to optimize the shock tracking and adaptive mesh refinement features (Kang & Jones 2006). The diffusion convection equation for the pitch-angle-averaged phase space distribution function for CR electrons in the testparticle regime was solved along with the usual gasdynamic conservation equations without including the dynamical feedback of the CR proton pressure. The synchrotron emission from CR electrons accelerated at these model shocks was estimated by adopting several models for magnetic field profile, B(r), and the electron energy spectrum, Ne(r, γe), from the DSA simulations.

Since the time scales for electron acceleration and cooling are much shorter than the dynamical time scale, the electron energy spectrum at the shock, fe(rs, p, t), has reached the steady state in the planar shock model, as can be seen in Figure 4. The slope of the volumeintegrated momentum distribution, Fe(p), is simply Q = q + 1 above the break momentum (p > pe,br). However, the spectral index of the volume-integrated synchrotron emission spectrum, Jν, steepens gradually to Aν = αinj + 0.5 over the broad frequency range, ~ (0.1 − 10)νbr. They are consistent with the predictions of the DSA theory in the test-particle limit.

In the spherical shock cases, the electron energy spectrum at the shock has reached the steady state defined by the instantaneous shock parameters (i.e., us(t), rs(t) and Ms(t)), as shown in Figure 5. So the synchrotron radiation spectra at the shock location could be described properly by the test-particle DSA predictions. However, the volume integrated spectra of Fe(p, t) and Jν(ν, t) both evolve differently from those of the planar shock, depending on the time-dependent evolution of the shock parameters. As a result, the slopes, Q and Aν, exhibit some nonlinear signatures that come from decreasing particle flux and declining shock Mach number (see Figure 5). This suggests that, in the case of evolving shocks, one needs to be cautious about interpreting observed radio spectra by adopting simple DSA models in the test-particle regime.

In addition, we have considered spherical shock models with decaying postshock magnetic fields, in which B(r) decreases behind the shock with a scale height of 100-150 kpc, as shown in Figure 6. Because inverseCompton scattering off the cosmic background photons with Brad ≈ 4.7 µG (for z = 0.2 adopted here as a reference epoch) provides the baseline cooling rate, the electron energy spectra do not depend sensitively on the B(r) profile. We also find that the impacts of the different B(r) profile on the spatial profile of synchrotron emission, jν(r), and its spectral index, αν(r), is rather minimal, since electrons in a somewhat broad range of γe contribute emission to jν at a given frequency. Of course, any nonlinear features due to the spatial variations in fe(r, p) and B(r) are mostly averaged out and may leave only subtle signatures in the volume integrated spectrum, Jν, and its spectral index, Aν.