An MHD study of SN 1006 and determination of the ambient magnetic field direction

Abstract

In this work, we employ a magnetohydrodynamic (MHD) numerical code to reproduce the morphology observed for SN 1006 in radio synchrotron and thermal X-ray emission. We introduce a density discontinuity, in the form of a flat cloud parallel to the Galactic plane, in order to explain the NW filament observed in optical wavelengths and in thermal X-rays. We compare our models with observations. We also perform a test that contrasts the radio emitting bright limbs of the supernova remnant against the central region, finding additional support to our results. Our main conclusion is that the most probable direction of the ambient magnetic field is on average perpendicular to the Galactic plane.

MHD, radiation mechanisms: general, methods: numerical, supernovae: individual: SN 1006, ISM: supernova remnants

1 INTRODUCTION

SN 1006 is a young supernova remnant (SNR) whose morphology in radio continuum can be described as a partial shell with a radius of 15 arcmin. Since most of the emission arises in two main bright arcs towards the NE and SW, this SNR has been classified as bilateral or barrel-shaped (BSNR; Kesteven & Caswell 1987; Gaensler 1998). Another characteristic of this source is that the thermal X-ray emission is uniformly distributed throughout the SNR, and is associated with the ejected material (see Cassam-Chenaï et al. 2008; Miceli et al. 2009). Recently, Acero, Ballet & Decourchelle (2007) carried out observations of this object with XMM–Newton, confirming the presence of a bright filament in the NW (see Dubner et al. 2002), suggesting that the expanding shell is colliding into a denser medium. Non-thermal X-ray emission is predominant in the NE and SW almost coincident with the synchrotron radio emission (Cassam-Chenaï et al. 2008).

Kesteven & Caswell (1987) propose two possible mechanisms for the production of the radio arcs: an intrinsic one, due, for example, to the interaction of internal jets with the SNR shell, and an extrinsic one, dominated by the characteristics of the surrounding interstellar medium (ISM).

Several theoretical works have lately studied the extrinsic mechanism for the formation of two opposed radio arcs in BSNRs through particle injection models together with analytical self-similar (for recent works see Cassam-Chenaï et al. 2008; Petruk et al. 2009) or magnetohydrodynamic (MHD) (see Orlando et al. 2007) numerical simulations.

Orlando et al. (2007) explore, through MHD simulations, different scenarios for explaining the asymmetries between the arcs in BSNRs. They propose two main mechanisms: the SNR expands in a non-uniform medium with an almost uniform magnetic field, and the SNR evolves in a uniform ISM with a non-uniform magnetic field. In either case, the density gradient or the magnetic field gradient is not aligned with the line of sight. To generate the synthetic synchrotron maps, they take into account three possible ways of particle injection: isotropic, quasi-parallel and quasi-perpendicular. The quasi-parallel (or perpendicular) nature of the injection is determined by the angle between the shock normal and the unshocked magnetic field (φ_Bn). According to Orlando et al. (2007) only the quasi-perpendicular and isotropic injection models reproduce radio images similar to the observed ones (in agreement with Fulbright & Reynolds 1990). However, it is important to know the direction of the unperturbed interstellar magnetic field (ISMF), in order to decide which injection mechanism is predominant.

Rothenflug et al. (2004) carried out an analysis of the radio-continuum and non-thermal X-ray emission of the remnant of SN 1006. By using an R_π/3 criterion and based mainly on the non-thermal X-ray emission distribution, Rothenflug et al. (2004) conclude that the ISMF around SN 1006 is mainly perpendicular to the NE and SW bright radio arcs (i.e. the ISMF is parallel to the Galactic plane).

A different conclusion was obtained by Petruk et al. (2009), who employed the analytical Sedov density and pressure profiles, where the aforementioned injection particles and self-similar magnetic field configuration are assumed to depend only on the shock obliquity through the compression factor in order to model the radio-continuum emission from SNRs. They analyse the orientation of the ISMF based on the azimuthal profile of the radio surface brightness on both observed and simulated radio maps of SN 1006. Applying this method, they found that the models that best reproduce both the observed radio morphology and the azimuthal radio brightness profile for the chosen aspect angles (angle between the ISMF and the line of sight) are the isotropic and quasi-perpendicular particle injection models. Therefore, they conclude that the most likely orientation of the ambient magnetic field at the galactic latitude of SN 1006 is parallel to the radio arcs, i.e. perpendicular to the Galactic plane, in contrast with the expected direction.

In this work, we follow a similar approach as that of Petruk et al. (2009), but employing a 2D axisymmetric MHD numerical simulation and solving for the magnetic field rather than making a priori assumptions on its final configuration. This work has a twofold objective: to explain the observed morphology in radio and X-rays, and to infer a possible structure for the ambient magnetic field, which can decide which injection particle mechanism is the most important to explain the non-thermal emission (radio and X-ray) of SNRs.

2 NUMERICAL SIMULATION

2.1 Code description

We carry out our simulations using a 2D Eulerian code (De Colle F. 2005; De Colle & Raga 2005, 2006), which solves the ideal MHD equations, including radiative cooling.

The equations are as follows.

1

2

3

4

where I is the identity matrix, p_tot=p_gas+B²/2 is the total pressure (thermal and magnetic), Λ(T) is the radiative cooling and the remaining symbols have the usual meaning. Equations (1–4) represent the mass, moment, energy and magnetic flux conservation, respectively.

As radiative loss term Λ(T), we use a tabulated coronal cooling function (Dalgarno & McCray 1972). The radiative cooling function is switched off at temperatures below 10⁴ K.

The constrained transport (CT) method (e.g. Tóth 2000) is used to conserve ∇·B= 0 to matching accuracy.

The MHD Riemann solver uses a second-order Runge–Kutta method for the time integration and a spatially second-order reconstruction of the primitive variables at the interfaces (except in shocks where it reduces at the first order). The fluxes are calculated using the Harten, Lax and Van Leer method (Harten, Lax & Van Leer 1983).

2.2 Initial conditions

All calculations were carried out on a 2D axisymmetrical grid with a physical size of 12 and 24 pc in the r- and z-directions, respectively, with a spatial resolution of 1.5 × 10⁻² pc. A uniform ISM of n₀= 5 × 10⁻² cm⁻³ (Acero et al. 2007) with a uniform magnetic field of 2 μG parallel to the z-axis.

The initial conditions consist of a sphere with radius 0.65 pc, located at the position (r₀, z₀) = (0, 12) pc, containing a total ejecta mass of M_★= 1.4 M_⊙, which is consistent with a typical Type Ia supernova. An inner region with a constant density containing 4/7 of the ejecta mass was imposed, while the outer region, containing the rest of the mass, has a density profile following a power-law density distribution ρ∝r⁻⁷ (see Jun, Jones & Norman 1996). This initial density distribution is adequate for Type Ia SNe.

To estimate the initial explosion energy we employed the following equation (Truelove & McKee 1999):

5

where E₅₁ is the initial energy in units of 10⁵¹ erg, D_2.2=D(kpc)/2.2 kpc is the distance to the object, R_b is the SNR observed radius in units of pc and t_yr is the SNR age in years. For the case of SN 1006, E₅₁= 2.05 considering R_b= 9.5 pc and t_yr= 1000. In the simulations, 95 per cent of the energy was set as kinetic energy and the remaining as thermal energy.

Three models were used for the simulations. In model M1, the remnant evolves in a uniform ISM with number density n₀. In the other models, and in order to simulate the thermal X-ray emission as well as the deformation of the NW region of SN 1006, we included an ISM-density discontinuity at 8.2 pc from the SNR centre, with a density three times higher than the ISM density. Two different directions for the ISMFs were considered by changing the position of the density discontinuity as shown in Fig. 1. In model M2 (M3) the density discontinuity is orthogonal (parallel) to the symmetry axis (we recall that the symmetry here is referred to the numerical simulation and not the bright radio and X-rays arcs of SN 1006). Due to the axisymmetry, the actual density discontinuity configuration is a ‘cylindrical cloud’ around the SNR.

Figure 1

Scheme of the initial numerical setup. The left-hand panel shows the setup for model M3, where the magnetic field is parallel to the density discontinuity, and the right-hand panel displays the setup for M2 where the magnetic field is perpendicular to the discontinuity. Notice that due to the axisymmetry assumption, the magnetic field has to be parallel to the symmetry axis and, for model M3, the density discontinuity is actually a cylinder about the SNR. Model M1 does not include the density discontinuity.

2.3 Synchrotron emission simulation

From the numerical results, we calculated the synthetic radio emission maps by integrating, along the line of sight, the radio emissivity given by (see Ginzburg & Syrovatskii 1965)

6

with K being the normalization of the electron distribution, B_⊥ the magnetic field component perpendicular to the line of sight, α the synchrotron spectral index [which is related with ζ, the index of the electron energy distribution N(E) ∝E^ζ, by α= (ζ− 1)/2] whose value is 0.6 for SN 1006 and ν the radio frequency.

The bilateral morphology observed in SNRs, such as SN 1006, can be generated by the fact that the efficiency of the particle acceleration at SNR shock front has a systematic dependence on the obliquity angle φ_Bn between the normal shock and the ISMF. The literature bifurcates in two main lines that argue for and against the quasi-parallel (φ_Bn= 0°) or the quasi-perpendicular (φ_Bn= 90°) models as the most efficient mechanism for accelerating particles at the remnant shock front. The quasi-parallel mechanism is associated with the classical diffusive shock acceleration (Blandford & Eichler 1987), while the quasi-perpendicular one (also called ‘shock drift mechanism’) can be a faster process because of particle acceleration by electric fields along the shock front (see e.g. Jokipii 1987). To include all most representative possible cases, in this work we study the isotropic (no φ_Bn dependence), the quasi-parallel (φ_Bn= 0°) and the quasi-perpendicular (φ_Bn= 90°) particle injection models. Following Orlando et al. (2007) (also Fulbright & Reynolds 1990), the isotropic case was obtained with equation (6) and K being a constant, whereas for the quasi-parallel and quasi-perpendicular cases, the factor K is proportional to cos²φ_Bs and sin²φ_Bs1, respectively. In order to identify the obtained simulated maps, we have added to M1, M2 and M3, the words ‘iso’, ‘par’ and ‘per’, for the case of isotropic, quasi-parallel and quasi-perpendicular particle injection.

2.4 Thermal X-ray emission

In order to compare with observations, the X-ray emission coefficient j_ν(n, T) was integrated along the line of sight to generate synthetic emission maps.

The emission coefficient was calculated as j_ν(n, T) =n²ψ(T) (considering the low-density regime), where n and T are the electronic density and temperature distribution obtained from our numerical simulations, while ψ(T) is a function that changes smoothly with T. ψ(T) was computed over the energy range [0.3–5] keV using the chianti atomic data base2 and its associated idl software (Dere et al. 1997; Landi & Phillips 2006). The interstellar absorption was added employing the Morrison & McCammon (1983) model, and taking into account a column density of N_H= 7 × 10²⁰ cm⁻² (Dubner et al. 2002; Acero et al. 2007). It was also assumed that the gas is in coronal ionization equilibrium (IEQ; Mazzotta et al. 1998).3

3 RESULTS AND DISCUSSION

Fig. 2 displays a comparison of the synthetic radio-continuum maps among all models, where in M1 the SNR is expanding in a uniform ISM (left-hand column), in M2 the SNR collides with a density discontinuity (at a distance of 8.2 pc from the centre of the SNR) perpendicular to the symmetry axis of the simulation emulating a magnetic field perpendicular to the Galactic plane, and in M3 the discontinuity is set at the same radial distance (8.2 pc) but parallel to the axis of symmetry, emulating a magnetic field parallel to the Galactic plane (the higher density flat cloud, shown in Fig. 1, is always parallel to the Galactic plane). The collision between the SNR blast wave and the density discontinuity occurs at t≃ 750 yr (in the simulation), which is in agreement with previous estimates (e.g. Acero et al. 2007).

Figure 2

Synthetic synchrotron maps for models M1 to M3 (from left to right) and with quasi-parallel (top row), quasi-perpendicular (middle row) and isotropic (bottom row) relativistic electron density distributions. The logarithmic colour-scale represents the synchrotron flux in arbitrary units (the same scale is used for all models), and the horizontal and vertical length-scales are given in units of pc. The direction of the ambient magnetic field is always along the symmetry axis towards +x. The vertical bright lines observed in M3 are artefacts introduced by the process for obtaining the simulated radio maps from a 2D axisymmetrical simulation.

By comparing the simulations with observations, we were able to discriminate some models, keeping M3par, M2per and M2iso (Fig. 3) as possible ones. Models M2par, M3iso and M3per were left aside because the flattened shape is parallel to the bright arcs, which is not observed in the SN 1006 case. The flattened shape at the NW and the formation of the two bright arcs are observed in M3par, M2per and M2iso cases. Models M2per and M2iso show a very similar brightness distribution as well as intensities, while model M3par shows an almost two orders of magnitude lower intensity for the maximum at the arcs in comparison to the maximum intensities obtained for the quasi-perpendicular and isotropic injection models. According to the position of SN 1006 in the Galaxy, one would expect a magnetic field perpendicular to the arcs (parallel to the Galactic plane), namely our models M2per and M2iso could be excluded from the analysis if this was the case.

Figure 3

Direct comparison between observation (at 1.5 GHz, bottom right) and models M3par (top-right), M2per (top-left) and M2iso (bottom-left). In order to improve the comparison, simulated synchrotron maps from models M2per and M2iso were tilted 130° clockwise in the plane of the sky, while the corresponding model M3par was tilted 40°, also clockwise in the plane of the sky.

Rothenflug et al. (2004) employ an R_π/3 criterion which is basically the ratio between the flux in the interior half of the SNR and the flux in the limbs. They obtained an R_π/3∼ 0.7 for the radio-continuum emission of SN 1006, and ∼0.3 and ∼0.1 for the non-thermal X-ray emission of the remnant in the [0.8–2.] and [2.0–4.5] keV bands, respectively. Rothenflug et al. (2004) conclude that the bright regions observed at the limbs are probably polar caps and not an equatorial belt. They also concluded that acceleration of relativistic electrons proceed preferentially where the magnetic field is parallel to the shock speed. We computed the ratio R_π/3 for synthetic radio-continuum maps for the three chosen models (see Fig. 3), and obtained ≃0.7, ≃0.72 and ≃0.14 for M2per, M2iso and M3par, respectively. These maps were generated considering an aspect angle of 90° (i.e. the ambient magnetic field lies on the plane of the sky). The M3par result exhibits a high contrast in radio emission between the central region and the polar cap regions, which does not agree with observations. This high contrast increases up to 3 if an aspect angle of 11° (i.e. the symmetry axis is tilted by 79° with respect to the plane of the sky) is considered when generating the simulated map. Such aspect angle was calculated by Petruk et al. (2009). The values of R_π/3 for models M2per and M2iso remain unchanged if this criterion is applied to synthetic maps of these models performed with the aspect angles obtained by Petruk et al. (2009) (64° and 70° for M2per and M2iso, respectively). In other words, the R_π/3 criterion is quite robust for both models.

Like Petruk et al. (2009), our model is based on the classical MHD Rankine–Hugoniot jump conditions for the parallel and perpendicular components (with respect to the normal shock) of the ISMF. This could not be valid if the pre-shock magnetic field was randomized by an efficient non-linear particle diffuse acceleration. In this case, the unperturbed ISMF can be amplified by orders of magnitude (Bell 2004). However, Petruk et al. (2009) found that, in spite of this, their predicted aspect angles remain almost unchanged for quasi-perpendicular and quasi-parallel injection models, after considering that the compression of a turbulent ISMF does not depend on the obliquity (Bohm limit Völk, Berezhko & Ksenofontov 2003) and that shocks at different initial obliquity under magnetic field turbulent amplification become perpendicular immediately upstream (Rakowski, Laming & Ghavamian 2008). The case of isotropic injection is discarded because it produces synthetic maps with no azimuthal flux variations, which is not observed (Petruk et al. 2009).

The simulated thermal X-ray map for model M2 is displayed in Fig. 4. The symmetry axis (⁠⁠) was clockwise tilted by 130° in the plane of the sky, in order to do a direct comparison with X-ray observations. This image was generated considering an aspect angle of 90°. A bright filament is formed at the top right region due to the collision with the density discontinuity. This collision does not affect the morphology of the radio emission. Rayleigh–Taylor features are observed in the internal region, close to the contact discontinuity. Due to the axisymmetry, the RT fingers become bright lines that cross the SNR face.4 In a 3D numerical simulation, the RT features would give a clumpy appearance to the ejecta emission close to the contact discontinuity interface (the contact discontinuity density distribution was shown in a 3D simulation by Miceli et al. 2009).

Figure 4

Simulated thermal X-ray emission corresponding to model M2. The logarithmic colour bar represents the X-ray flux in units of erg s⁻¹ cm⁻² st⁻¹. The horizontal and vertical axes are given in units of pc.

Fig. 5 displays the emission measure (EM) distribution obtained by solving for the integral ⁠, where dl lies along the line of sight. This figure was also clockwise tilted by 130° in the plane of the sky and the aspect angle is 90°. The simulated EM can be compared to the parameter norm defined as 10⁻¹⁴ EM Ω/4π, where Ω is the solid angle of the emitting gas (Acero et al. 2007). We obtained an average value for the bright X-ray filament region of 5 × 10⁻⁴ cm⁻⁵ which doubles the value reported by Acero et al. (2007).

Figure 5

Simulated EM maps corresponding to model M2. The logarithmic colour bar represents the EM in units of cm⁻⁵. The horizontal and vertical axes are given in units of pc.

4 CONCLUSIONS

In this work we present MHD axisymmetric numerical simulation of the remnant of SN 1006, taking into account different directions for the ISMF. Based on the numerical results, synthetic radio and thermal X-ray emission maps were constructed.

According to our results, and considering the morphology of the radio and X-rays synchrotron emission, the best fit is given by models M2per and M2iso, where the magnetic field direction is perpendicular to the arcs, and therefore perpendicular to the Galactic plane. This result is in agreement with those obtained through analytical models in Petruk et al. (2009). These results are reinforced when applying the R_π/3 criterion, since models M2per and M2iso yield a ratio in accordance with the observations (of the order of 0.7) even after correcting for the aspect angle given by Petruk et al. (2009). At the same time, the model that could explain the morphology observed with a magnetic field parallel to the Galactic plane (M3par) can be rejected with this test for synchrotron radiation.

Furthermore, our simulated X-ray maps show a bright filament such as that observed in the NW of SN 1006 (Acero et al. 2007), and the internal emission (from the ejected gas close to the contact discontinuity interface) revealed the RT features. Finally, our estimates of the EM are in good agreement with the values reported by Acero et al. (2007), confirming that the collision scenario, while capable of reproducing the observed thermal X-ray emission, does not affect the synchrotron emission.

In summary, our results show that around the SN 1006 remnant, the ISMF direction is mainly perpendicular to the Galactic plane, in agreement with the previous work (Petruk et al. 2009). Based in this fact, the best fit with radio observations of this remnant is achieved if the injection of the relativistic particles at the forward SNR shock wave is isotropic or quasi-perpendicular. This result could be different if non-linear diffuse acceleration process takes place at the SNR shock front. In this case, the model based on isotropic particle injection is no longer valid because the predicted synthetic radio maps would show a complete shell morphology for any position angle. Notwithstanding, models based on quasi-perpendicular particle injection remain valid, with no changes in the predicted aspect angle (Petruk et al. 2009), if this process is independent of obliquity angle φ_Bn (Bohm limit).

1

φ_Bs is the angle between the normal shock and the post-shock magnetic field.

2

The chianti data base and its associated idl procedures are freely available at http://www.solar.nrl.navy.mil/chianti.html, http://www.arcetri.astro.it/science/chianti/chianti.html, http://www.damtp.cam.ac.uk/users/astro/chianti/chianti.html

3

The IEQ is an approximation for the case of young SNRs such as the SN 1006 giving a lower limit for the X-ray flux.

4

If an aspect angle difference of 90° is employed, these bright lines and the bright filament would show an elliptical morphology.

Authors acknowledge an anonymous referee for her/his very useful suggestions and comments, which helped in improving the previous version of this manuscript. PFV acknowledges support from grants CONACyT 79744 and DGAPA IN119709. EMR is partially supported by grants UBACyT X482, PIP 114-200801-00428 (CONICET) and ANPCYT-PICT-2007-00902. Authors thank Enrique Palacios Boneta (cómputo-ICN) for maintaining and supporting the LINUX servers where the numerical simulations were carried out.