© ESO, 2015

1. Introduction

SN 1006 is a young bilateral supernova remnant (SNR) with two opposite bright limbs on the detected morphologies in the radio, non-thermal X-rays, and γ-rays. The hard X-rays from the NE and SW limbs (north is up and east to the left) are predominately non-thermal with featureless power law spectra (Koyama et al. 1995; Bamba et al. 2003; Long et al. 2003; Miceli et al. 2012), and the soft thermal X-rays interior of the remnant are associated with the shocked ejecta (Acero et al. 2007; Miceli et al. 2009). Based on the results of the deep observations by Chandra (Cassam-Chenaï et al. 2008) and XMM-Newton (Miceli et al. 2012) on SN 1006, the hard X-rays in the band 2−4.5 keV originate in the shocked interstellar medium, and the morphologies show distinct bumps along the shell. For instance, a big bump is located on the NE limb, and several small ones on the SW arc. Moreover, the radio observations with the Australia Telescope Compact Array (ATCA) and the Very Large Array (VLA; Cassam-Chenaï et al. 2008) obtain a morphology similar to the one in the hard X-rays. In the TeV γ-rays, the high energy stereoscopic system (HESS) also obtains a consistent bipolar morphology with those either in the hard X-rays or in the radio band (Acero et al. 2010).

It is widely accepted that SN 1006 is a Type-Ia SNR resulting from the thermonuclear supernova on 1006 AD. The remnant has an angular size of 30′, which corresponds to a radius of ~9.5 pc for a distance of 2.18 ± 0.08 kpc derived by Winkler et al. (2003). The environmental matter of SN 1006 is generally tenuous, but inhomogeneous components still exist based on the multiband observations. Acero et al. (2007) investigated the gas density around SN 1006 by modeling the post-shock thermal X-ray emission. They find that the ambient medium for the remnant is tenuous with a density of ~0.05 cm^-3 around the southeastern (SE) rim. It seems that this density applies to the rest of the remnant except the northwestern (NW) part, where a bright H_α filament is detected, and it can be deduced that the shock in the NW rim is interacting with dense material (e.g., Raymond et al. 2007). Moreover, Miceli et al. (2014) performed spatially resolved spectral analysis of the SW bright limb based on the deep observations with XMM-Newton and the archive HI data, and an isolated cloud is derived on the SW region of the remnant.

Studies of the radiative properties of SN 1006 can give insight into the particle acceleration process in SNR shocks. Reynoso et al. (2013) provide studies of radio polarization on SN 1006 based upon the observations with VLA and ATCA, and the results also confirm that the particle acceleration and the generation of magnetic turbulence are efficient for quasi-parallel shocks. Theoretically, the diffusive shock acceleration is widely accepted as the main process for accelerating electrons and ions to relativistic energies at the shocks of SNRs (Drury 1983; Malkov & Drury 2001). Recently, Caprioli & Spitkovsky (2014a) have studied the particle acceleration and the amplification of magnetic field at non-relativistic shocks based on hybrid simulations. They find the acceleration depends heavily on the angle between the shock normal and the preshock magnetic field. The acceleration of ions is efficient for parallel or quasi-parallel strong shocks. For quasi-perpendicular shocks, ions can be accelerated via the shock drift acceleration, but the acceleration is rather ineffective because the particles experiencing the shock-drift acceleration are advected downstream.

The morphology of a SNR relates heavily to the ambient environment that is influenced by the progenitor system (see Wang & Han 2012, for a review of the progenitors of Type Ia supernovae). A non-uniform medium has been used to explain the detected asymmetries of SNRs through either analytical treatments (Hnatyk & Petruk 1999) or numerical simulations (e.g., Orlando et al. 2007; Fang et al. 2014; Yang et al. 2015). Moreover, Guo et al. (2012) employed 2D magnetohydrodynamic (MHD) simulations to investigate the amplification of magnetic fields for a SNR propagating into a turbulent plasma with the 2D Kolmogorov-like power spectrum. When the shock interacts with the ambient turbulent medium, the front can become rippled, and the turbulent flow further amplify the magnetic field.

The dynamical and radiative properties of SN 1006 have been extensively studied through numerical simulations. Scheiter et al. (2010) performed 2D MHD simulations to obtain the synthetic synchrotron maps for SN 1006. In their simulations, a flat cloud parallel to the Galactic plane is assumed to produce the NW filament in the observed images both in optical wavelengths and in thermal X-rays. Bocchino et al. (2011) investigated the expansion of a SNR in a medium with a non-uniform magnetic field and obtained the resulting non-thermal radio maps, which can be compared with the detected images for SN 1006. Their results show that the quasi-parallel injection efficiency is approved, and a gradient of the magnetic field toward the Galactic plane is derived by comparing the resulting morphology of the synchrotron radiation with the image detected at 1 GHz.

Our 3D MHD simulations were meant to study the matter distribution of the environment for the SNR 1006 by comparing the synchrotron morphology of non-thermal X-ray emission with the observed ones. The background plasma is assumed to be turbulent with a relatively small amplitude to represent the real environment of the remnant; moreover, a cavity with a lower density is used to reproduce the bump in the NE limb, and the other three cavities are adopted to explain the observed bumps in the SW arc in the non-modification scenario. Moreover, the influence of the shock modification on the synchrotron morphology of the remnant due to the efficient shock acceleration of cosmic rays is investigated. For SN 1006, the observed morphology of synchrotron emission can be generally reproduced using the model. The organization of this paper is as follows. In Sect. 2, we describe the numerical model and the simulation setup. The results are presented in Sect. 3. Some summary and discussions are given in Sect. 4.

2. Numerical model

2.1. Numerical code

The propagation of the shock front for SN 1006 in the magnetized plasma with the turbulent density and the turbulent magnetic field is simulated using the time-dependent ideal MHD equations of mass, momentum, and energy conservation in the 3D Cartesian coordinate system. In the simulation, we do not consider the radiative cooling, the thermal conduction, the nonlinear effects due to acceleration of charged particles, and the magnetic field amplification by the accelerated particles. The ideal MHD equations are $\begin{array}{ccc}\frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{v}\mathrm{\right)}& \mathrm{=}& \mathrm{0}\mathit{,}\\ \frac{\mathit{\partial \rho }{v}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{(}\mathit{\rho }{vv}\mathrm{-}{BB}\mathrm{)}\mathrm{+}\mathrm{\nabla }{\mathit{P}}^{\mathrm{\ast }}& \mathrm{=}& \mathrm{0}\mathit{,}\\ \frac{\mathit{\partial E}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathrm{\right(}\mathit{E}\mathrm{+}{\mathit{P}}^{\mathrm{\ast }}\mathrm{)}{v}\mathrm{-}{B}\mathrm{(}{v}\mathrm{·}{B}\mathrm{\left)}\mathrm{\right)}& \mathrm{=}& \mathrm{0}\mathit{,}\\ \frac{\mathit{\partial }{B}}{\mathit{\partial t}}\mathrm{-}\mathrm{\nabla }\mathrm{\times }\mathrm{\left(}{v}\mathrm{\times }\mathit{B}\mathrm{\right)}& \mathrm{=}& \mathrm{0}\mathit{,}\end{array}$where P^∗ = P + B²/ 2 is the total pressure (thermal pressure, P, and magnetic pressure, B²/ 2), and E is the total energy density $\begin{array}{ccc}\mathit{E}\mathrm{=}\frac{\mathit{P}}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{\mathit{v}}^{\mathrm{2}}\mathrm{+}\frac{{\mathit{B}}^{\mathrm{2}}}{\mathrm{2}}\mathit{,}& & \end{array}$t is the time, ρ = 1.27m_Hn is the mass density with an assumption of a 10:1 H:He ratio, m_H is the mass of a hydrogen atom, n the number density of ions, γ = 5/3 is the adiabatic index for nonrelativistic gas in the non-modification case, v the gas velocity, and B the magnetic field. The temperature of the plasma involved in the simulation is always higher than 10⁴ K. Moreover, a mean atomic mass μ = 0.6 is used for the ionized medium, which means that ρ = μm_H(n + n_e), where n_e is the number density of electrons.

The 3D MHD simulation of the SNR shock front propagating through the turbulent plasma is performed with the PLUTO code (Mignone et al. 2007). In the calculation, the MHD equations are solved with the cell-centered finite-volume scheme with high-order Godunov methods. The magnetic field is evolved with the constrained transport technique to preserve the divergence-free condition, and the fluxes are calculated using the HLL Riemann solver (Miyoshi & Kusano 2005). Moreover, the simulations are performed on a cubic Cartesian domain that extends between − 12 pc to 12 pc for each of the x, y, and z directions. The whole numerical grid has 512³ zones, and a resolution of ~0.05 pc is obtained in each of the directions.

2.2. Simulation setup

Shortly after the supernova explosion with a time t, the ejecta consists of an inner core with a uniform density within a radius r_c and an outer part with a density following the power-law distribution (Fraschetti et al. 2010), ${\mathit{\rho }}_{\mathrm{ej}}\mathrm{\left(}\mathit{t,r}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ {\mathit{\rho }}_{\mathrm{c}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{(}\mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{c}}{\mathrm{)}}^{\mathrm{-}\mathit{s}}& \mathrm{if}\mathit{r}\mathit{>}{\mathit{r}}_{\mathrm{c}}\mathit{,}\\ {\mathit{\rho }}_{\mathrm{c}}\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{if}\mathit{r}\mathit{<}{\mathit{r}}_{\mathrm{c}}\mathit{.}\end{array}$(5)The index s is usually adopt to be 7 for a Type-Ia remnant. Moreover, we assume that the ejecta expands freely to the circumstellar medium, and the velocity for the matter of the ejecta at r is ${v}\mathrm{=}{r}\mathit{/}\mathit{t}\mathit{.}$(6)With an initial kinetic energy E_ej, a mass M_ej for the ejecta, and a mass fraction η = 3/7 of the outer part of the ejecta for Type-Ia SNRs (Colgate & Mckee 1969), r_c can be calculated with (Jun & Norman 1996) ${\mathit{r}}_{\mathrm{c}}\mathrm{=}{\mathit{R}}_{\mathrm{ej}}{\left[\mathrm{1}\mathrm{-}\frac{\mathit{\eta }\mathrm{(}\mathrm{3}\mathrm{-}\mathit{s}\mathrm{)}{\mathit{M}}_{\mathrm{ej}}}{\mathrm{4}\mathit{\pi }{\mathit{\rho }}_{\mathrm{0}}{\mathit{R}}_{\mathrm{ej}}^{\mathrm{3}}}\right]}^{\mathrm{1}\mathit{/}\mathrm{(}\mathrm{3}\mathrm{-}\mathit{s}\mathrm{)}}\mathit{,}$(7)where R_ej = 0.8 pc is adopted as the initial radius of the ejecta, and ρ₀ is the density of the uniform component of the background medium. The density of the inner core ρ_c and the velocity at the outer boundary of the ejecta v₀ can be obtained with (Jun & Norman 1996; Fang & Zhang 2012) ${\mathit{\rho }}_{\mathrm{c}}\mathrm{=}\frac{\mathrm{3}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\eta }\mathrm{)}{\mathit{M}}_{\mathrm{ej}}}{\mathrm{4}\mathit{\pi }{\mathit{r}}_{\mathrm{c}}^{\mathrm{3}}}\mathit{,}$(8)and ${\mathit{v}}_{\mathrm{0}}\mathrm{=}{\mathit{E}}_{\mathrm{ej}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\begin{array}{c}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\frac{\mathrm{2}\mathit{\pi }{\mathit{\rho }}_{\mathrm{c}}{\mathit{r}}_{\mathrm{c}}^{\mathrm{5}}}{\mathrm{5}{\mathit{R}}_{\mathrm{ej}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{2}\mathit{\pi }{\mathit{\rho }}_{\mathrm{0}}{\mathit{R}}_{\mathrm{ej}}^{\mathrm{3}}\mathrm{[}\mathrm{1}\mathrm{-}\mathrm{(}{\mathit{R}}_{\mathrm{ej}}\mathit{/}{\mathit{r}}_{\mathrm{c}}{\mathrm{)}}^{\mathit{s}\mathrm{-}\mathrm{5}}\mathrm{]}}{\mathrm{5}\mathrm{-}\mathit{s}}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(9)To simulate the real environment around SN 1006, the background plasma, in which the shock front propagates, is assumed to be turbulent with a relatively smaller amplitude. In the simulation, the 3D Kolmogorov-like power spectrum, which is consistent with the observations of interstellar turbulence (Armstrong et al. 1995; Chepurnov & Lazarian 2010), is adopted for the turbulent density and magnetic field, i.e., the power spectrum follows $\mathit{P}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{\propto }\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathrm{(}\mathit{k}{\mathit{L}}_{\mathrm{c}}{\mathrm{)}}^{\mathrm{11}\mathit{/}\mathrm{3}}}\mathit{,}$(10)where k is the magnitude of the wavevector, and L_c = 3 pc is the coherence length. The magnetic field consists of a uniform background and a turbulent component; i.e., ${B}\mathrm{\left(}\mathit{x,y,z}\mathrm{\right)}\mathrm{=}{{B}}_{\mathrm{0}}\mathrm{+}\mathit{\delta }{B}\mathrm{\left(}\mathit{x,y,z}\mathrm{\right)}\mathit{,}$(11)and the turbulent component can be calculated with (Giacalone & Jokippi 1999) $\mathit{\delta }{B}\mathrm{\left(}\mathit{x,y,z}\mathrm{\right)}\mathrm{=}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{m}}}\mathit{A}\mathrm{\left(}{\mathit{k}}_{\mathit{n}}\mathrm{\right)}\left[\cos \mathrm{\left(}{\mathit{\alpha }}_{\mathit{n}}\mathrm{\right)}{{x̂}}_{\mathit{n}}^{\mathrm{\prime }}\mathrm{+}\mathit{i}\sin \mathrm{\left(}{\mathit{\alpha }}_{\mathit{n}}\mathrm{\right)}{{ŷ}}_{\mathit{n}}^{\mathrm{\prime }}\right]\exp \mathrm{\left(}\mathit{i}{\mathit{k}}_{\mathit{n}}{\mathit{z}}_{\mathit{n}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{,}$(12)where $\left(\begin{array}{c}\\ {\mathit{x}}^{\mathrm{\prime }}\\ {\mathit{y}}^{\mathrm{\prime }}\\ {\mathit{z}}^{\mathrm{\prime }}\end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ \cos {\mathit{\theta }}_{\mathit{n}}\cos {\mathit{\phi }}_{\mathit{n}}& \cos {\mathit{\theta }}_{\mathit{n}}\sin {\mathit{\phi }}_{\mathit{n}}& \mathrm{-}\sin {\mathit{\theta }}_{\mathit{n}}\\ \mathrm{-}\sin {\mathit{\phi }}_{\mathit{n}}& \cos {\mathit{\phi }}_{\mathit{n}}& \mathrm{0}\\ \sin {\mathit{\theta }}_{\mathit{n}}\cos {\mathit{\phi }}_{\mathit{n}}& \sin {\mathit{\theta }}_{\mathit{n}}\sin {\mathit{\phi }}_{\mathit{n}}& \cos {\mathit{\theta }}_{\mathit{n}}\end{array}\right)\left(\begin{array}{c}\\ \mathit{x}\\ \mathit{y}\\ \mathit{z}\end{array}\right)\mathit{,}$(13)θ_n and φ_n represent the direction of propagation of the wave mode n with the wave number k_n, polarization α_n and phase β_n. The amplitude of the wave mode n is (Giacalone & Jokippi 1999) ${\mathit{A}}^{\mathrm{2}}\mathrm{\left(}{\mathit{k}}_{\mathit{n}}\mathrm{\right)}\mathrm{=}{\mathit{\sigma }}_{\mathrm{B}}^{\mathrm{2}}\frac{\mathrm{\Delta }{\mathit{V}}_{\mathit{n}}}{\mathrm{1}\mathrm{+}\mathrm{(}{\mathit{k}}_{\mathit{n}}{\mathit{L}}_{\mathit{c}}{\mathrm{)}}^{\mathrm{11}\mathit{/}\mathrm{3}}}{\left[\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{m}}}\frac{\mathrm{\Delta }{\mathit{V}}_{\mathit{n}}}{\mathrm{1}\mathrm{+}\mathrm{(}{\mathit{k}}_{\mathit{n}}{\mathit{L}}_{\mathit{c}}{\mathrm{)}}^{\mathrm{11}\mathit{/}\mathrm{3}}}\right]}^{-1}\mathit{,}$(14)where ${\mathit{\sigma }}_{\mathrm{B}}^{\mathrm{2}}$ is the wave variance and the normalization factor $\mathrm{\Delta }{\mathit{V}}_{\mathit{n}}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{k}}_{\mathit{n}}^{\mathrm{2}}\mathrm{\Delta }{\mathit{k}}_{\mathit{n}}$. In the simulation, N_m = 1000 discrete wave modes are summed with random phases and polarizations to represent the isotropic turbulence. The maximum wavelength of these modes is chosen to be λ_max = 15 pc, whereas the minimum wavelength λ_min = 0.1 pc, which is twice as large as the grid-cell size in the simulation.

For the fluctuations of the density, a log-normal probability is adopted (Giacalone & Jokippi 2007), i.e., $\mathit{n}\mathrm{\left(}\mathit{x,y,z}\mathrm{\right)}\mathrm{=}{\mathit{n}}_{\mathrm{0}}\exp \mathrm{(}{\mathit{f}}_{\mathrm{0}}\mathrm{+}\mathit{\delta f}\mathrm{)}$(15)where δf has the 3D Kolmogorov-like power spectrum as δB with a wave variance ${\mathit{\sigma }}_{\mathit{f}}^{\mathrm{2}}$.

Based on the detected image in the hard X-rays with Chandra (Cassam-Chenaï et al. 2008), the bump at the NE arc of SN 1006 has an azimuth extension of ~50° (10° ≤ θ ≤ 60°), and each of the three smaller bumps on the SW limb has an extension of ~12°. To reproduce the observed bumps in the NE and SW limbs, one large and three smaller spherical cavities with lower densities are adopted in the background plasma. The temperatures of the plasma in these cavities are determined by assuming that they are in pressure equilibrium with the outer medium. The first cavity with a lower density n_c1 is located in the NE part of the simulation domain at (x_c1, y_c1, z_c1) with a radius of R_c1, and the other cavities have a density of n_ci at (x_ci, y_ci, z_ci) with a radius of R_ci (i = 2,3,4).

The supernova explodes at the center of the simulation domain with the explosion energy E_ej = 1.5 × 10⁵¹ erg, the mass of the ejecta M_ej = 1.4 M_⊙. The background plasma has a temperature of T_ISM = 10⁴ K, a turbulent density with n₀ = 0.035 cm^-3, f₀ = 0, δf = 0.41, and a turbulent magnetic field with B₀ = 3μG, B_x = 0, B_y = −B₀cos(π/ 4), B_z = B₀sin(π/ 4), σ_B=0.1. Moreover, the centers of the four spherical cavities are located at the plane x_c = 0 pc, and the other parameters are shown in Table 1.

2.3. Synchrotron emission

The relativistic electrons accelerated by the forward shock in SNRs can emit photons with energies from radio to X-rays through synchrotron radiation. Assuming these electrons have a power law distribution with a cutoff around the maximum energy E_max, i.e., ${\mathit{N}}_{\mathrm{e}}\mathrm{=}{\mathit{K}}_{\mathrm{e}}{\mathit{E}}^{\mathrm{-}{\mathit{s}}_{\mathrm{e}}}\exp \left(\mathrm{-}\frac{\mathit{E}}{{\mathit{E}}_{\max }}\right)\mathit{,}$(16)where s_e is the index, the emissivity of the synchrotron radiation for the electrons can be approximated by (Katsuda et al. 2010) ${\mathit{i}}_{\mathrm{syn}}\mathrm{\left(}{\mathit{E}}_{\mathrm{ph}}\mathrm{\right)}\mathrm{=}{\mathit{C}}_{\mathit{j}}\mathrm{\left(}{\mathit{s}}_{\mathrm{e}}\mathrm{\right)}{\mathit{K}}_{\mathrm{e}}{\mathit{B}}_{\mathrm{\perp }}^{\mathrm{(}{\mathit{s}}_{\mathrm{e}}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{/}\mathrm{2}}{\mathit{E}}_{\mathrm{ph}}^{\mathrm{-}\mathrm{(}{\mathit{s}}_{\mathrm{e}}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{/}\mathrm{2}}\exp \left(\mathrm{-}\sqrt{\frac{{\mathit{E}}_{\mathrm{ph}}}{{\mathit{E}}_{\mathrm{cut}}}}\right)\mathit{,}$(17)where K_e is a factor describing the amount of the relativistic electrons, C_j is a constant, B_⊥ is the magnetic field perpendicular to the line of sight, E_ph is the energy of the photon, and ${\mathit{E}}_{\mathrm{cut}}\mathrm{\propto }{\mathit{E}}_{\max }^{\mathrm{2}}{\mathit{B}}_{\mathrm{\perp }}$ is the exponential cutoff energy of the synchrotron radiation. The above assumption on the energy distribution of the electrons does not include any post-shock evolution of these relativistic particles. Moreover, we ignore the variation in the E_cut along the remnant outline, which have been indicated by the X-ray observations (Rothenflug et al. 2004; Miceli et al. 2014). Non-thermal X-rays from a SNR are usually characterized by a cutoff frequency that corresponds to the maximum energy of the electrons accelerated by the SNR shock. Based on the theory of the diffusive shock acceleration, E_max is determined by the shock speed, age, and various energy loss mechanisms (Reynolds 1998). The effect of the maximum energy, which varies both with time and with the surface of the remnant, on the synchrotron spectrum must be taken into account to accurately synthesize the X-ray emission from the SNR (Orlando et al. 2011).

The studies from the hybrid simulations by Caprioli & Spitkovsky (2014a) show that the acceleration efficiency of ions at nonrelativistic shocks depend heavily on the angle between the shock normal and the direction of the preshock magnetic field. The acceleration of the ions is efficient for parallel and quasi-parallel shocks (Caprioli & Spitkovsky 2014a). Assuming thalethe electrons share the same feature as the ions, the acceleration of the electrons follow the quasi-parallel model. The normalization of the relativistic electrons is multiplied by a factor of cos²θ_Bn, where θ_Bn is the angle between the preshock magnetic field and the shock normal. In our simulations, the amplitude of the turbulent component of the ambient magnetic field is relatively small, therefore calculating θ_Bn does not take the turbulent component of the magnetic field into account.

Fig. 1 Spatial distribution of the density (ρ) from the 3D simulation with an age of 1000 yr in the plane z = 0 (left panel), x = 0 (middle panel), and y = 0 (right panel). The mass density is in units of 1.67 × 10-24 g cm-3.

For the spherical shock of a SNR, the accelerated particles can escape from the acceleration process beyond the contact discontinuity since the magnetic field strength is negligible within the discontinuity; moreover, the relativistic electrons suffer serious energy loss due to the synchrotron radiation behind the shock. As a result, the high-energy electrons distributed closely to the SNR shock. According to the thermal injection theory on the diffusive acceleration regime, more charged particles can be injected into the acceleration with a higher density of the plasma immediately downstream of the shock (e.g., Blasi et al. 2005). It is difficult to accurately determine the number density of the accelerated particles at a location downstream of the shock because it depends on the density distribution of the background medium swept by the shock and the detail transportation of these relativistic particles, especially for the shock disturbed by inhomogeneous medium with a complex periphery. In this paper, the relativistic electrons are supposed to distribute in the thin area along the forward shock, and we simply assume that K_e in the SNR is proportional to the local density n at a certain age of the remnant since generally more charged particles can be injected in the acceleration process with a higher density at the shock. Therefore, when taking this effect and the quasi-parallel acceleration into account, K_e ∝ ncos²θ_Bn is adopted in the simulation. The cutoff energy of the detected non-thermal X-rays from SN 1006 varies azimuthly along the non-thermal limb, and a mean value of ~0.2 keV is derived based on the observations with Chandra (Allen et al. 2008) and XMM-Newton (Rothenflug et al. 2004; Miceli et al. 2014). We simply assume the cutoff energy E_cut = 0.2 keV for all directions although the measured cut-off energy varies both radially and along the remnant outline (Rothenflug et al. 2004; Allen et al. 2008; Miceli et al. 2014).

Magnetic field can be greatly amplified at SNR shocks caused by the acceleration of charged particles. The hybrid simulations also indicated that the amplification of the magnetic field correlates with the efficiency of the ion acceleration. Caprioli & Spitkovsky (2014b) found that the upstream magnetic field for non-relativistic shocks could be significantly amplified, and the strength of the amplified magnetic field generated via the resonant streaming instability have the relation (Lagage & Cesarsky 1983; Amato & Blasi 2006; Caprioli & Spitkovsky 2014b), $\mathrm{⟨}{\mathit{B}}_{\mathrm{tot}}\mathit{/}{\mathit{B}}_{\mathrm{0}}{\mathrm{⟩}}^{\mathrm{2}}\mathrm{\approx }\mathrm{3}{\mathit{\zeta }}_{\mathrm{cr}}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{M}\end{array}$, where $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{M}\end{array}$ is the Mach number of the upstream medium in the shock reference frame, ζ_cr is the normalized pressure in the accelerated particles, and ζ_cr is related to the acceleration efficiency of the non-thermal particles. The generation of the magnetic field turbulence cannot be consistently tackled by the MHD simulations, which ignore the acceleration of the charged particles, and we simply assume that the strength of amplified magnetic field is proportional to the square root of the local number density of the accelerated particles in calculating the emissivity of the synchrotron radiation, i.e., ${\mathit{B}}_{\mathrm{\perp }}\mathrm{\propto }{\mathit{K}}_{\mathrm{e}}^{\mathrm{1}\mathit{/}\mathrm{2}}$.

3. Simulation results

Figure 1 shows the spatial distribution of the density in the simulation domain after a time of 1000 yr. The slices in the plane z = 0, x = 0, and y = 0 are shown in the left, middle, and right panels, respectively. In the beginning of the simulation, the density and the pressure of the ejecta are high enough to drive a strong shock in the turbulent background. The ambient plasma is continually heated and compressed by the shock. Moreover, a reverse shock is produced since the shocked circumstellar medium pushes back on the ejecta (Truelove & McKee 1999). The inner ejecta expands freely until it is swept by the reverse shock, and both the density and the pressure of the unshocked ejecta decrease quickly during the expansion. The forward and the reverse shocks are separated by the contact discontinuity, where the Rayleigh-Taylor instability is triggered. Figure 2 indicates the distribution of the averaged density on the radial distance to the center of the explosion.

The forward shock is located at a radius of ~9.5 pc, which is consistent with the detected extension of the remnant for a distance of ~2.18 kpc (Winkler et al. 2003). As a consequence of the turbulent background, the shape of the forward shock is slightly distorted. In the forward shock, the inhomogeneous distribution of the density is obtained because the materials swept by the shock in various directions are different from each other. In some places, knots and filaments are produced since more material has been swept by the shock than in the other locations. In the middle panel of Fig. 1, the distribution of the density shows one bump in the NE rim of the remnant and three smaller ones in the SW limb because the shock speed is not very decelerated as the shock sweeps the lower density cavities. The radius of the forward shock from the simulation is indicated in Fig. 3 as a function of the azimuth angle (θ). In the ranges 10°<θ< 50° and 220°<θ< 270°, the radius significantly exceeds the averaged one owing to the existence of the cavities, and the other dispersion is the result of the turbulent property of the ambient medium. The resulting azimuth extensions of the bumps on the limbs are consistent with the detected ones due to the NE breakout with an extension of ~50° and with the three ones in the SW limb with an averaged extension of ~12°.

Fig. 2 Density profile as a function of the radial distance to the center of the remnant.

Fig. 3 Shock radius as a function of the azimuth angle when the LOS is aligned with the x axis.

Fig. 4 Strength of the magnetic field as a function of the radial distance to the center of the remnant.

Fig. 5 Spatial distribution of the pressure in units of 1.67 × 10-8 dyn cm-2 from the 3D simulation with an age of 1000 yr. The left, middle, and right panels are for the plane z = 0, x = 0, and y = 0, respectively.

The magnetic filed strength in the remnant as a function of the radial distance is shown in Fig. 4. The magnetic field is amplified at the location of the contact discontinuity due to the Rayleigh-Taylor instability. Moreover, the magnetic field both in the shocked ejecta within the contact discontinuity and in the interior of the remnant becomes negligible because of the expansion of the ejecta, and the similar feature can also be seen in Fig. 5 for the spatial distribution of the pressure.

Fig. 6 Synthesis morphologies of the non-thermal synchrotron radiation with energies from 2−4.5 keV when the LoS is aligned with the z (left panel), x (middle panel), and y (right panel) axis. The intensity is in arbitrary units.

The distribution of the electrons accelerated by the forward shock is assumed to follow the quasi-parallel scenario, i.e., N_e ∝ ncos²θ_Bn. Because the energetic electrons can escape from the shock owing to weak confinement of the magnetic field and because these electrons suffer strong energy loss due to the synchrotron radiation, we suppose that the electrons emitting the high-energy synchrotron emission are distributed close to the shock front; i.e., the energetic electrons are distributed in the locations with pressure p> 10^-8 dyn cm^-2. Furthermore, assuming that the strength of the amplified magnetic field is proportional to the number density of the accelerated electrons, we show the resulting synchrotron morphologies of the hard X-rays with energies from 2−4.5 keV in the amplified magnetic field when the line of sight (LoS) is aligned with the z, x, and y axes in the left, middle, and right panels in Fig. 6, respectively. The image with the LoS along the x direction shows a bilateral shape due to the quasi-parallel scenario of the acceleration process, in which the efficiency of the acceleration is high at locations with the shock normal parallel or antiparallel to the background magnetic field. However, the non-thermal X-rays in the NW and SE parts of the remnant are faint because of the inefficient particle acceleration. Although the morphology shows a general circular shape, some disturbances still exist owing to the turbulent feature of the background plasma. Besides the NE and SW bumps, the morphologies of the remnant both in radio and in X-rays have the knotty and filamentary structures, and it deviates from the circular shape since the SNR forward shock is transported into the turbulent medium.

The images in the non-thermal X-rays detected with XMM-Newton (Rothenflug et al. 2004; Miceli et al. 2012) and Chandra (e.g., Long et al. 2003; Cassam-Chenaï et al. 2008; Winkler et al. 2014) have distinct bumps along the SNR shell. The image when the LoS is aligned with the x direction show a two-arc structure, which is consistent with the observed morphology in the energy band 2−4.5 keV by Chandra and XMM-Newton. The largest bump located at the NE arc of the remnant can be reproduced in the simulation by imposing a spherical region with a density of 0.4n₀ and a radius of 2.5 pc. On the opposite limb, the smaller bumps can be explained with the three spherical cavities with a density of 0.2n₀ swept by the shock front. The radial profiles of the intensities for the azimuth range 0°<θ< 90° and 180°<θ< 270° of the synchrotron radiation are shown in Figs. 7 and 8, respectively. The intensity peaks at a radius of around 9 pc, and the other peaks beyond this distance result from the lower density cavities.

Fig. 7 Radial profile for the intensity of the synchrotron radiation with the azimuth range 0°<θ< 90° when the LoS is aligned with the x axis. The intensity is in arbitrary unit.

Fig. 8 Radial profile for the intensity of the synchrotron radiation with the azimuth range 180°<θ< 270°. The others are the same as in Fig. 7.

The radius of the cavity determines the azimuth extension of the bump on the limb of SN 1006, and the distance of the bump edge to the center of the remnant depends heavily on the density of the cavity. The simulation in this paper shows that the radius of the NE cavity is about 2.5 pc. The resulting morphologies of the remnant for n_c1 = 0.2n₀ and 0.8n₀ are shown in Figs. 9 and 10, respectively. Moreover, Fig. 11 indicates the azimuth distribution of the radius for the NE limb of the remnant with n_c1 = 0.2n₀, 0.4n₀, and 0.8n₀. The morphology with n_c1 = 0.4n₀ is more like the hard X-ray images detected with Chandra or XMM-Newton than like n_c1 = 0.2n₀ and 0.8n₀.

The expansion rates of the NW limb of SN 1006 based on the Chandra observations have been proposed by Katsuda et al. (2009). For the NE limb, the proper motion has an averaged value of ~ yr^-1, which is consistent with that derived from the radio observations for the entire remnant with a value of ~ yr^-1 (Moffett et al. 1993). The derived shock velocities for the NE limb based on the Chandra observations range from 4000 to 6000 km s^-1 as a function of the azimuth angle for a distance of 2.2 kpc. Recently, Winkler et al. (2014) have given the measurement of the X-ray proper motion for the entire periphery of SN 1006. The derived expansion velocity in the NW limb is consistent with Katsuda et al. (2009), and the results show that the velocity varied significantly with azimuth along the entire rim of the remnant. Especially, the highest velocity, which is about 1.5 times the averaged one along the NE periphery, was found to be located at the SE limb. The velocities of the forward shock distributed over the azimuth angle from the simulation are illustrated in Fig. 12. Along the NE periphery with the azimuth range 0°<θ< 90°, the forward shock expands with velocities around 5000 km s^-1 with dispersions due to the existence of the cavity and the turbulent medium, and the variation in the velocity over the azimuth angle is generally consistent with the results derived from the Chandra observations (Katsuda et al. 2009; Winkler et al. 2014). For the SW limb, the measured velocities in Winkler et al. (2014) for the expansion of the outer most X-rays vary significantly in the range from ~4000 to ~8000 km s^-1 with the azimuth angle 180°<θ< 270°. In the range of 200°<θ< 240°, the variation tendency obtained from the simulation is generally consistent with the observation, although the exact values have some differences. Moreover, a value of >7000 km s^-1 is derived for the shock velocity along the SE periphery and at some places on the SW limb, and these expansion rates can be explained as X-rays emitted from the ejecta overtaking the forward shock. For the NW periphery, the observations show that the expansion velocities are almost <4000 km s^-1, which results from a denser medium located in the NW region around the remnant. However, because the NW limb is not luminous in non-thermal X-rays/radio, we do not take the denser medium into account in this paper.

Fig. 9 Synthesis morphologies of the non-thermal synchrotron radiation for nc1 = 0.2n0. Others are the same as in Fig. 6.

Fig. 10 Synthesis morphologies of the non-thermal synchrotron radiation for nc1 = 0.8n0. Others are the same as in Fig. 6.

Figure 13 shows the simulated morphology for the remnant evolved in a uniform ambient medium. The bilateral shape with luminous NE and SW limbs is also held with the quasi-parallel scenario for the acceleration process when the LoS is aligned with the x axis. With the uniform background plasma, the distribution of the intensity that changes regularly and smoothly along the NE/SW limb is inconsistent with the observed ones that have the knotty and filamentary structures and deviations from the circular shape in the radio and X-rays.

A SNR shock can be modified due to the acceleration of the cosmic rays, and the distance between the contact discontinuity and the forward shock is significantly decreased in the case of efficient acceleration (Blondin & Ellison 2001; Wang 2011; Warren & Blondin 2013). As a result, a fraction of the ejecta can approach closer or even overtake the forward shock, and the remnant outline is disturbed by the hydrodynamic instabilities. We also investigate the influence of the shock modification on the resulting morphology of the remnant by using an artificially lower adiabatic index as in Warren & Blondin (2013). Figure 14 illustrates the resulting distribution of density of the remnant evolved in a uniform medium with an adiabatic index of 1.2, which is also adopted in Warren & Blondin (2013) to study the effect of the cosmic ray acceleration on the morphology of the forward shock. The higher compressibility increases the density downstream of the forward shock compared with γ = 5/3, and the Rayleigh-Taylor fingers can extend more closely to the forward shock to influence the edge of the shock as the remnant evolved in the medium. Assuming the accelerated electrons are distributed close to the forward shock in the regions with the pressure above 0.45, the resulting morphologies of the synchrotron emission along the z, x, and y directions are shown in the left, middle, and right panels in Fig. 15, respectively. Along the z and y directions, bubbles are visible in the morphologies since the hydrodynamical instabilities overtake the forward shock. Although some of the fingers can extend beyond the shock boundary to produce bubbles visible on the forward shock as in Warren & Blondin (2013), the extensions of these bubbles are usually not as large as the bump on the NW limb of SN 1006. The big bump at the NE rim of SN 1006 cannot therefore be explained as these bubbles. However, the small bumps on the SW limb can also be interpreted as the bubbles when the Rayleigh-Taylor fingers effectively overtake the forward shock in the case of the efficient particle acceleration at the forward shock.

4. Summary and discussion

Based on the observations in the non-thermal X-ray and radio bands, the SNR SN 1006 indicates a generally bilateral shape with two opposite bright limbs. This feature suggests that the circumstellar medium around the remnant is more uniform than the other SNRs with complex shapes, such as RX J0852.0-4622 and RX J1713.7-3946. However, some inhomogeneities still exist along the shell of SN 1006, especially the distinct bumps that are obvious in both the radio and the X-ray morphologies. Besides the bumps, the inhomogeneities cannot be explained by isotropic supernova ejecta evolved in a uniform medium, in which the circular shape of the remnant is not disturbed and the knotty and filamentary morphology cannot be reproduced.

Fig. 11 Shock radius as a function of the azimuth angle in the range 0°<θ< 100° for nc1 = 0.2n0 (square), 0.4n0 (dot), and 0.8n0 (triangle).

Fig. 12 Velocity of the forward shock as a function of the azimuth angle for the range 0°<θ< 90° (upper panel) and 180°<θ< 270° (lower panel).

To give a more realistic explanation to the detected hard X-ray/radio morphology of SN 1006, we assumed that the ambient medium is turbulent with a relatively small amplitude to ensure that the uniformity, which is derived from the bilateral shape of the remnant, of the background plasma is not significantly disturbed. The observed deviation of the circular shape of the hard X-ray morphologies with knotty and filamentary structure can be explained as the SNR evolved in the turbulent medium, although we did not fit the exact deviations from the observations by exploring the parameter space on the turbulent background plasma. Furthermore, several spherical cavities with lower densities were adopted in the background plasma to reproduce the distinct bumps in the detected morphologies of the synchrotron emission.

Fig. 13 Synthesis morphologies of the non-thermal synchrotron radiation with energies from 2−4.5 keV when the LoS is aligned with the z (left panel), x (middle panel), and y (right panel) axes in the scenario of the remnant evolved in a uniform medium. The intensity is in arbitrary units.

Fig. 14 Spatial distribution of the density (ρ) from the 3D simulation with an age of 1000 yr for γ = 1.2. The others are the same as Fig. 1.

Fig. 15 Synthesis morphologies of the non-thermal synchrotron radiation with energies from 2−4.5 keV when the LoS is aligned with the z (left panel), x (middle panel), and y (right panel) axis in the scenario of the remnant evolved in a uniform medium with γ = 1.2. The intensity is in arbitrary units.

In this paper, the knotty and filamentary structure of the remnant was generated as the forward shock wave expanding into the inhomogeneous medium. Some other mechanisms that can lead to the structure have been proposed. For example, Orlando et al. (2012) investigated the influence of the ejecta clumping on the density distribution of the post-shock region for a Type-Ia SNR. They find that the reduction of the distance between the forward shock and the contact discontinuity is also produced (even more efficiently) by the ejecta fragmentation (Orlando et al. 2012). Moreover, the growth of the RT instabilities is greatly enhanced by the interactions among the clumps, and the enhanced instabilities send the ejecta closer to or even beyond the forward shock. Then the outline of the remnant is disturbed, and a knot and filamentary structure are produced.

The ejecta clumping can lead to the occurrence of several protrusions due to the clumps overtaking the forward shock, and the model in Orlando et al. (2012) successfully explained the occurrence of protrusions in some SNRs, such as SN 1006 (e.g., Rakowski et al. 2011) and Tycho’s SNR (Hwang & Gotthelf 1997; Velazquez et al. 1998). However, the bump on the NE limb in the morphology of SN 1006 has a relatively large scale, and it seems that they should not be produced by the RT instability fingers overtaking the forward shock. In this paper, the big bump at the NE rim of the remnant is formed by the cavity in the background medium. The smaller bumps like those on the SW limb of SN 1006 can be explained either by cavities swept by the forward shock or by the bubbles resulting from the Rayleigh-Taylor fingers that overtake the forward shock in the scenario of efficient cosmic ray acceleration. The physical origin of these cavities is not clear.

The recent hybrid simulation shows that underdense cavities, which have scales comparable to the gyroradius of the highest-energy particles in the simulation, can be produced in the shock precursor of the quasi-parallel shocks with efficient particle acceleration due to the streaming instability (Caprioli & Spitkovsky 2014a). The turbulent medium with those cavities around the remnant imposed in our simulation can be considered as an approximation of lots of the cavities upstream of the quasi-parallel shock due to the particle acceleration. Moreover, lower density cavities around a Type-Ia remnant can also be generated as blowing of the wind from the progenitor system (e.g., Kashi & Soker 2011; Shen et al. 2012; Broersen et al. 2014).

The hydrodynamic instabilities growing at the contact discontinuity are related to the numerical resolution, and the plasma properties in the mixing region between the contact discontinuity and the forward shock are influenced by these instabilities (Blondin & Ellison 2001). To increase the precision of the resolution in a 3D simulation is very expensive, so we adopted 512³ zones to simulate the evolution of the remnant in the medium with a resolution of ~0.05 pc obtained in the simulations. In the non-modification scenario with γ = 5/3, the instabilities cannot extend very close to the forward shock, and the resulting morphology of emission is less influenced. Moreover, even in the modification case, the bubbles on the surface of the remnant stemming from the instabilities overtaking the forward shock are unlikely to have enough extensions to explain the big bump located on the NW limb of SN 1006. Therefore, our results are not influenced by the resolution adopted in this paper.

The two assumptions on the acceleration of charged particles at SNR shocks and the amplification of the magnetic field were adopted to obtain the distribution of the energetic electrons in the SNR shell. The acceleration and the amplification cannot be self-consistently tackled with MHD simulations that ignore the acceleration process of charged particles by the shock. Both the recent hybrid simulations and the radio observations on SN 1006 support the quasi-parallel scenario for the acceleration at SNR shocks. Therefore, the quasi-parallel scenario was adopted, and two more assumptions on the strength of the amplified magnetic field and the spatial distribution of the energetic electrons producing the non-thermal X-rays were used to achieve the morphology of the synchrotron emission at any given energy.

We paid attention to the bumps at the NE/SW limb, which is bright in non-thermal X-rays, by simulating the remnant as an explosion evolved in a turbulent medium with either spherical cavities with the lower densities or modifications of the morphologies of the remnant due to the efficient cosmic ray acceleration. In the NW part of the remnant, a denser medium interacting with the shell exists based on the bright filament in soft X-rays (Dubner et al. 2002; Acero et al. 2007; Nikolić et al. 2013), and Scheiter et al. (2010) has investigated the influence of the NE cloud on the dynamical properties of the remnant by 2D MHD simulations. However, the NW region is not luminous in the non-thermal X-rays, and the dense cloud around the NW limb is not taken into account in this paper.

Acknowledgments

We acknowledge the referee for the valuable comments that helped improve this paper. This work is partially supported by the National Natural Science Foundation of China (NSFC 11433004, 11173020), by the Key Project of Chinese Ministry of Education (212160), by the Doctoral Fund of Ministry of Education of China (RFDP 20115301110005), by the West Light Foundation of the Chinese Academy of Sciences (W8090303), by the Open Research Program of the Key Laboratory for the Structure and Evolution of Celestial Objects, Chinese Academy of Sciences (OP201402), and by the Top Talents Program of Yunnan Province, China.

References

All Tables

All Figures

	Fig. 1 Spatial distribution of the density (ρ) from the 3D simulation with an age of 1000 yr in the plane z = 0 (left panel), x = 0 (middle panel), and y = 0 (right panel). The mass density is in units of 1.67 × 10-24 g cm-3.
In the text

	Fig. 2 Density profile as a function of the radial distance to the center of the remnant.
In the text

	Fig. 3 Shock radius as a function of the azimuth angle when the LOS is aligned with the x axis.
In the text

	Fig. 4 Strength of the magnetic field as a function of the radial distance to the center of the remnant.
In the text

	Fig. 5 Spatial distribution of the pressure in units of 1.67 × 10-8 dyn cm-2 from the 3D simulation with an age of 1000 yr. The left, middle, and right panels are for the plane z = 0, x = 0, and y = 0, respectively.
In the text

	Fig. 6 Synthesis morphologies of the non-thermal synchrotron radiation with energies from 2−4.5 keV when the LoS is aligned with the z (left panel), x (middle panel), and y (right panel) axis. The intensity is in arbitrary units.
In the text

	Fig. 7 Radial profile for the intensity of the synchrotron radiation with the azimuth range 0°<θ< 90° when the LoS is aligned with the x axis. The intensity is in arbitrary unit.
In the text

	Fig. 8 Radial profile for the intensity of the synchrotron radiation with the azimuth range 180°<θ< 270°. The others are the same as in Fig. 7.
In the text

	Fig. 9 Synthesis morphologies of the non-thermal synchrotron radiation for nc1 = 0.2n0. Others are the same as in Fig. 6.
In the text

	Fig. 10 Synthesis morphologies of the non-thermal synchrotron radiation for nc1 = 0.8n0. Others are the same as in Fig. 6.
In the text

	Fig. 11 Shock radius as a function of the azimuth angle in the range 0°<θ< 100° for nc1 = 0.2n0 (square), 0.4n0 (dot), and 0.8n0 (triangle).
In the text

	Fig. 12 Velocity of the forward shock as a function of the azimuth angle for the range 0°<θ< 90° (upper panel) and 180°<θ< 270° (lower panel).
In the text

	Fig. 13 Synthesis morphologies of the non-thermal synchrotron radiation with energies from 2−4.5 keV when the LoS is aligned with the z (left panel), x (middle panel), and y (right panel) axes in the scenario of the remnant evolved in a uniform medium. The intensity is in arbitrary units.
In the text

	Fig. 14 Spatial distribution of the density (ρ) from the 3D simulation with an age of 1000 yr for γ = 1.2. The others are the same as Fig. 1.
In the text

	Fig. 15 Synthesis morphologies of the non-thermal synchrotron radiation with energies from 2−4.5 keV when the LoS is aligned with the z (left panel), x (middle panel), and y (right panel) axis in the scenario of the remnant evolved in a uniform medium with γ = 1.2. The intensity is in arbitrary units.
In the text