A 3D MHD simulation of SN 1006: a polarized emission study for the turbulent case

Abstract

Three dimensional magnetohydrodynamical simulations were carried out in order to perform a new polarization study of the radio emission of the supernova remnant SN 1006. These simulations consider that the remnant expands into a turbulent interstellar medium (including both magnetic field and figuredensity perturbations). Based on the referenced-polar angle technique, a statistical study was done on observational and numerical magnetic field position-angle distributions. Our results show that a turbulent medium with an adiabatic index of 1.3 can reproduce the polarization properties of the SN 1006 remnant. This statistical study reveals itself as a useful tool for obtaining the orientation of the ambient magnetic field, previous to be swept up by the main supernova remnant shock.

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

1 INTRODUCTION

SN 1006 is a young supernova remnant (SNR) that has gained interest due to the large observational data available in a wide wavelength range, therefore serving as an ideal laboratory when trying to understand the observed morphology in each frequency. This remnant is classified as part of the bilateral SNR group, whose main characteristic is the presence of two bright and opposite arcs in radio-frequencies.

The morphology and emission of this type of SNRs are largely determined by the environment in which they evolve, sometimes having an irregular or clumpy background and/or with the presence of density and magnetic fields gradients.

The interstellar medium (ISM) is known to be turbulent, as inferred from observational data (Lee & Jokipii 1976) and local information on the interplanetary medium (Jokipii, Lerche & Schommer 1969). Such turbulent medium, with a Kolmogorov-like power spectrum that spans over a large range of spatial scales (Guo et al. 2012), gives rise to Rayleigh-Taylor (RT) instabilities during the fast expansion of the SNR, which, in turn, affects the radio emission and is imprinted in the measured Stokes parameters.

Jun & Jones (1999) performed 2D magnetohydrodynamic (MHD) simulations of the evolution of a young Type Ia SNR interacting with an interstellar cloud. They found that the interaction produced RT instabilities causing amplification of the magnetic field, and subsequently changing the synchrotron brightness. This work emphasized the importance of including a turbulent medium when trying to match observations. In this direction, Balsara, Benjamin & Cox (2001) performed more realistic 3D simulations of an SNR evolving in a turbulent background. They found that the variability of the emission during the early phase of SNR evolution gives information about the turbulent environment. At the same time, that the interaction with such a turbulent environment generates an enhanced turbulent post-shock region that has observational consequences, as it affects the acceleration of relativistic particles.

Guo et al. (2012) employed 2D MHD simulations to study the magnetic field amplification mechanisms during the propagation of the SNR blast wave. They found that, for high-resolution simulations, magnetic field growth at small scales occurs efficiently in two distinct regions: one related to the shock amplification, and the other, to the RT-instabilities at the contact discontinuity between the shocked ejecta and the shocked ISM.

Fang & Zhang (2012) carried out a study similar to that of Guo et al. (2012) but considering different adiabatic indexes to account for both the diffusive shock acceleration and the escape of the accelerated particles from the shock. They found that a smaller effective adiabatic index produces a larger magnetic field enhancement. They extended their study considering a scenario where the expanding ejecta interacts with a large, dense clump and found, in addition, to an increase in the complexity of the magnetic field, the development of a non-thermal emitting filament. In a subsequent study, Fang, Yu & Zhang (2014) performed 3D MHD simulations considering a turbulent medium with a Kolmogorov-like power spectrum with parameters resembling the SNR RX J0852.0−4622. Their X-ray and γ-ray synthetic maps showed the expected rippled morphology with a broken circular shape along the shell, nicely reproducing the observations.

In general, the diffusive shock acceleration is driven by magnetic fields strong enough to generate non-thermal emission. Therefore, any structure or dynamic interaction that results in the amplification of the magnetic field has important consequences on the observed brightness.

According to the non-thermal radio-emission morphology, SN 1006 has been classified as a bilateral or barrel-shaped SNR (Kesteven & Caswell 1987), showing two main bright arcs towards the NE and SW. The non-thermal X-ray emission is also predominant in the NE and SW rims (Cassam-Chenaï et al. 2008).

SN 1006 is located at a high Galactic latitude and is therefore thought to be evolving in a fairly homogenous interstellar medium. In Schneiter et al. (2010) and Schneiter et al. (2015), we made use of this assumption and presented studies of the system in 2D and 3D, respectively. In the latter work, we further introduced synthetic maps of the Stokes Q and U parameters, which allowed a better and more direct comparison with the observations (Reynoso, Hughes & Moffett 2013).

Recently, West, Safi-Harb & Ferrand (2016a) and West et al. (2016b) (see also, West et al. 2016a,b) studied radio emission for a sample of Galactic axisymmetric SNRs, including SN 1006. They concluded that the quasi-perpendicular acceleration mechanism can successfully fit the observed morphology of the remnants in their sample, although SN 1006 seems to be the exception, since the quasi-parallel mechanism agrees better with observations.

In the this work, we use the same SNR initialization as in the two previous works of Schneiter et al. (2010) and Schneiter et al. (2015) but relaxing the assumption of a homogeneous ISM by considering that the SNR evolves into a turbulent medium as it was shown in Yu et al. (2015) (see also Fang et al. 2014).

The goal of this work is to propose an alternative method to determine the position angle of the ambient magnetic field, based on the polar-referenced angle method employed by Reynoso et al. (2013) and Schneiter et al. (2015).

The numerical model is presented in Section 2. For completeness, we recall the basic model employed by Schneiter et al. (2015) in Section 2.1 and in Section 2.2 we explain how the turbulence is introduced. Section 3 explains how synthetic radio maps were obtained and the results are presented in Section 4. The discussion and final conclusions are given in Section 5.

2 NUMERICAL MODEL

The numerical model employed is basically the same as in Schneiter et al. (2010, 2015), with the additional assumption of a turbulent background, where both the density and magnetic fluctuations have a 3D Kolmogorov-like power spectrum, as will be explained further down.

To simulate the evolution of SN 1006, we employed the Mezcal code (De Colle & Raga 2006; De Colle, Raga & Esquivel 2008; De Colle et al. 2012). This code solves the full set of ideal MHD equations in a Cartesian geometry with an adaptive mesh, and includes a cooling function to account for radiative losses (De Colle & Raga 2006).

The computational domain is a cube of 24 pc per side, which we will denote as (x, y, z), and is discretized on a five level binary grid with a maximum resolution of 4.7 × 10⁻² pc. All the outer boundaries were set to outflow condition.

2.1 SNR initial conditions

In what follows, we describe the physical parameters and setup used to simulate the SNR, which are the same as those presented in Schneiter et al. (2015).

A supernova explosion is initialized by the deposition of E₀ = 2.05 × 10⁵¹ erg in a radius of R₀ = 0.65 pc located at the centre of the computational domain. The energy is distributed such that 95 per cent of it is kinetic and the remaining 5 per cent is thermal.

The ejected mass was distributed in two parts: an inner homogeneous sphere of radius r_c containing 4/7ths of the total mass (M_★ = 1.4 M_⊙) with a density ρ_c, and an outer shell containing the remaining 3/7ths of the mass following a power law (ρ ∝ r⁻⁷) as in Jun & Norman (1996a). The velocity has an increasing linear profile with r, which reaches a value of v₀ at r = R₀. The parameters ρ_c, r_c and v₀ are functions of E₀, M_★ and R₀, and were computed using equations (1)–(3) of Jun & Norman (1996a).

2.2 The turbulent background

Jun & Jones (1999) suggested that the turbulent structures of the brightest radio emission correlate well with the magnetic field in general. For this reason, we performed new simulations with a more realistic magnetized turbulent ISM. In what follows, we will make use of equations (10)–(15) of Yu et al. (2015), which we reproduce below. To introduce the turbulent background, we assumed a 3D Kolgomorov-like power spectrum fluctuation for both the density and the magnetic field, similar to that of Jokipii (1987) (see also Fang et al. 2014; Yu et al. 2015). The power spectrum follows:

\begin{equation} P \propto \frac{1}{1+(k L_c)^{11/3}}. \end{equation}

(1)

To generate the turbulence, we assume a coherence length of L_c = 3 pc and a wavenumber |$k=\frac{2 \pi }{L}$|⁠. Our simulations are computed with N_m = 903 wave modes and L varies between L_min = Δx and L_max = L_sim, being Δx and L_sim the cell and the computational domain sizes, respectively.

The initial magnetic field is given by

\begin{equation} {\boldsymbol B}(x,y,z)={\boldsymbol B}_0+{\delta \boldsymbol B}, \end{equation}

(2)

where |${\boldsymbol B_0}= B_0(x) \hat{y}$| and the magnetic field perturbation is

\begin{equation} {\delta \boldsymbol B}=\Re \bigg [\sum _{n=1}^{N_m} A(k_n)\left( \cos \alpha _n \hat{x}^{\prime } + i \sin \alpha _n \hat{y}^{\prime } \right) \exp (i k_n z_n^{\prime })\bigg ], \end{equation}

(3)

with

\begin{equation} A^2(k_n)=\sigma _{\rm B}^2 \frac{\Delta V_{\rm n}}{1+(k_nL_c)^{11/3}}\sum _{n=1}^{N_m} \left[\frac{\Delta V_{\rm n}}{1+(k_nL_c)^{11/3}} \right]^{-1}, \end{equation}

(4)

where |$\sigma ^2_{\rm B}$| is the wave variance of the magnetic field and the normalization factor, ΔV_n, is given by

\begin{equation} \Delta V_{\rm n}=4 \pi k_n^2 \Delta k_n. \end{equation}

(5)

Both systems (x΄, y΄, z΄) and (x, y, z) are related by

\begin{equation*} {\left(\begin{array}{ccc}x^{\prime } \\ y^{\prime } \\ z^{\prime } \end{array}\right)} = {\left(\begin{array}{ccc}\cos \theta _n \cos \phi _n& \cos \theta _n \sin \phi _n & -\sin \theta _n \\ -\sin \phi _n & \cos \phi _n & 0 \\ \sin \theta _n \cos \phi _n & \sin \theta _n \sin \phi _n & \cos \theta _n \end{array}\right)} {\left(\begin{array}{ccc}x \\ y \\ z \end{array}\right)} \end{equation*}

where θ_n and ϕ_n represent the direction of propagation of the wave mode n with the wavenumber k_n and polarization α_n.

For the density fluctuation, we employed the same log-normal distribution as in Giacalone & Jokipii (2007):

\begin{equation} n(x,y,z)=n_0 \exp ({f_0+\delta f}), \end{equation}

(6)

where f₀ is a constant and δf is the density perturbation with a wave variance |$\sigma ^2_{\rm d}$|⁠. Both δf and |${\delta \boldsymbol B}$| have the same 3D Kolmogorov-like power spectral index.

The turbulent environment temperature and number density were set to T₀ = 10⁴K and n₀ = 5 × 10⁻²cm⁻³, respectively. The wave variances for magnetic field and density were set as |$\sigma ^2_{\rm B}=0.25\, B_0^2$| and |$\sigma ^2_{\rm d}=0.4 n_0^2$|⁠.

As mentioned above, one of the purposes of this work is to be able to study the effect of a perturbation on some observed parameters. To compare with the observations, the synthetic maps were performed for an integration time corresponding to 1000 yr (the age of SN 1006).

3 SYNTHETIC EMISSION MAPS

In Table 1, we list the runs that were carried out. With these runs, we want to analyse the effects of both the reduction of the value of the adiabatic index γ and/or the inclusion of a turbulent background with a Kolgomorov-like power spectrum for both density and magnetic field. These runs are labeled as Rip, where i is 1 for the case of simulations with an adiabatic index γ = 5/3 and 2 for γ = 1.3. The label ‘p’ indicates the presence or absence of turbulence in the background magnetic field and density: p = a indicates no turbulence in either variable; p = b includes turbulence in both; p = c only considers a turbulent |$\boldsymbol B$|⁠; and p = d only considers a turbulent ρ.

Table 1.

Runs carried out in this paper with their corresponding hypothesis: adiabatic index, turbulent magnetic field and/or turbulent density distribution.

run	γ	δ|$\boldsymbol B$|	δρ
R1a	5/3	no	no
R1b	5/3	yes	yes
R2a	1.3	no	no
R2b	1.3	yes	yes
R2c	1.3	yes	no
R2d	1.3	no	yes

Synthetic radio emission maps were obtained from the numerical results. In order to compare them with the observations, the computational domain (denoted as xyz system), they were rotated with respect to the ‘image’ system (x_iy_iz_i system). At the beginning both systems coincide. The x_iy_i plane of the ‘image’ system is the plane of the sky, being |$\hat{y}_i$| the direction towards the ‘North’, |$-\hat{x}_i$| the direction towards the ‘East’ and the line of sight (hereafter LoS) the |$\hat{z}$| direction. Then, synthetic maps were obtained after performing three rotations on the computational domain. First, a rotation around |$\hat{y}_i$| was applied (the ‘North’) with an angle |$\varphi _{y_i}$|⁠. Secondly, a rotation in |$\varphi _{x_i}$| was carried out around the |$\hat{x}_i$|-direction. Finally, the computational domain was rotated by |$\varphi _{z_i}$| around the |$\hat{z}_i$|-direction (the LoS). To get the best agreement with the observations, after several tests these angles were set as −15°, −30° and 60° for |$\varphi _{x_i}$|⁠, |$\varphi _{y_i}$| and |$\varphi _{z_i}$|⁠, respectively. In this way, the projection of the unperturbed magnetic field |${\boldsymbol B_0}$| was tilted by 60° with respect to the |$\hat{y}_i$| direction.

3.1 Synchrotron emissivity

For each point (x_i, y_i, z_i) of the SNR, we can obtain the synchrotron specific intensity as (see Cécere et al. 2016, and references therein):

\begin{equation} j(x_i,y_i,z_i,\nu ) \propto K \rho \ v^{4\alpha }\ B_{\perp }^{\alpha +1} \nu ^{-\alpha }, \end{equation}

(7)

where ν is the observed frequency, B_⊥ is the component of the magnetic field perpendicular to the LoS, α is the spectral index that was set to 0.6 for this object, and ρ and v are the density and velocity of the gas, respectively. The coefficient K includes the obliquity dependence, being either proportional to sin ²Θ_Bs for the quasi-perpendicular case or to cos ²Θ_Bs for the quasi-parallel case. The angle Θ_Bs is the angle between the shock normal and the post-shock magnetic field (Fulbright & Reynolds 1990). For the case of SN 1006, Bocchino et al. (2011) and Reynoso et al. (2013) showed that the observed morphology of this remnant can be explained by considering the quasi-parallel case. Schneiter et al. (2015) carried out a polarization study of this remnant, based on 3D MHD simulations, finding that only the quasi-parallel case can successfully reproduced the observed morphology of the Stokes parameter Q (see fig. 3 of Schneiter et al. 2015). For these reasons only the quasi-parallel case was considered in this work. The synthetic synchrotron emission maps are obtained by integrating j(x_i, y_i, z_i, ν) along the LoS or z_i-axis, i.e.

\begin{equation} I(x_i,y_i,\nu )=\int _{\mathrm{LoS}} j(x_i,y_i,z_i,\nu ) dz_i, \end{equation}

(8)

being x_i and y_i the coordinates in the plane of the sky.

3.2 Stokes parameters and position-angle distribution maps

With the purpose of comparing the numerical results with the observations, synthetic maps of the Stokes parameters Q and U were calculated as follows (Jun & Norman 1996b; Clarke, Burns & Norman 1989; Schneiter et al. 2015):

\begin{equation} Q(x_i,y_i,\nu )=\int _{\mathrm{LoS}} f_0 j(x_i,y_i, z_i,\nu ) \cos \left[ 2\phi (x_i,y_i,z_i)\right] {\rm d}z_i, \end{equation}

(9)

\begin{equation} U(x_i,y_i,\nu )=\int _{\mathrm{LoS}} f_0 j(x_i,y_i, z_i,\nu ) \sin \left[ 2\phi (x_i,y_i,z_i)\right] {\rm d}z_i, \end{equation}

(10)

where ϕ(x_i, y_i, z_i) is the position angle of the local electric field in the plane of the sky and f₀ is the degree of linear polarization, which is a function of the spectral index α:

\begin{equation} f_0=\frac{\alpha +1}{\alpha + 5/3}. \end{equation}

(11)

The position angle of the local magnetic field ϕ_B(x_i, y_i, z_i) is known from the simulations and the position angle of the local electric field ϕ(x_i, y_i, z_i) can be obtained from ϕ_B(x_i, y_i, z_i) by applying a |$\frac{\pi }{2}$| rotation and correcting for Faraday rotation, as

\begin{equation} \phi (x_i,y_i,z_i)=\phi _{\mathrm{B}}(x_i,y_i,z_i)-\frac{\pi }{2}+\Delta \chi _{\mathrm{F}}. \end{equation}

(12)

The Faraday correction term Δχ_F is given by

\begin{equation} \Delta \chi _{\mathrm{F}}=\mathrm{RM}\ \lambda ^2, \end{equation}

(13)

where RM is the rotation measure (in units of rad m⁻²) and λ is the wavelength of the observations (given in metres). In this study, only the external RM arising in the foreground ISM is considered. Bandiera & Petruk (2016) explored the effects of internal RM on the polarized synchrotron emission of SNRs. For the case of SN 1006, we have carried out a rough estimation of the internal RM, obtaining a value close to 10 per cent of the reported one by Reynoso et al. (2013) who obtained a RM of 12 rad m⁻² for this remnant.

The linearly polarized intensity is computed as

\begin{equation} I_{\rm P}(x_i,y_i,\nu )= \sqrt{Q(x_i,y_i,\nu )^2+U(x_i,y_i,\nu )^2} \end{equation}

(14)

and the maps of the polarization angle distribution were computed as follows:

\begin{equation} \chi (x_i,y_i)=\frac{1}{2}\tan ^{-1}(U(x_i,y_i,\nu )/Q(x_i,y_i,\nu )) \end{equation}

(15)

Similarly, the distribution of the magnetic field orientation can be calculated as

\begin{equation} \chi _{\mathrm{B}}(x_i,y_i)=\frac{1}{2}\tan ^{-1}(U_{\mathrm{B}}(x_i,y_i,\nu )/Q_{\mathrm{B}}(x_i,y_i,\nu )), \end{equation}

(16)

where U_B(x_i, y_i, ν) and Q_B(x_i, y_i, ν) are calculated with the equations (9) and (10) and replacing ϕ(x_i, y_i, ν) with ϕ_B(x_i, y_i, ν).

4 RESULTS

4.1 Synthetic polarization maps

Fig. 1 shows a comparison between the observed map of the linearly polarized intensity of SN 1006 (upper panel), and the synthetic polarization maps obtained for runs with γ = 5/3 (run R1b, middle panel) and 1.3 (run R2b, bottom panel) by using equation (14). The observed map was constructed by combining data obtained with the Australia Telescope Compact Array and Very Large Array at 1.4 GHz, as explained in Reynoso et al. (2013). The polarized emission image, as well as the distribution of the Q and U parameters, were convolved to a 15 arcsec beam. The two models used to construct the synthetic maps consider that the SNR is expanding into a turbulent medium with perturbations in both density and magnetic field. These synthetic maps are shown in arbitrary units. Two opposite bright, noisy arcs are observed in both synthetic maps, having a striking resemblance to the observations. Reducing the value of γ produces a 5 per cent smaller expansion radius of the remnant as measured on the linearly polarized intensity maps.

Figure 1.

Comparison of the linearly polarized intensity at 1.4 GHz between the observations (panel (a)) and synthetic maps of a SNR expanding into a turbulent medium with adiabatic index γ = 5/3 (model R1b, panel (b)) and γ = 1.3 (model R2b, panel(c))

Following the analysis developed by Schneiter et al. (2015), we compare the observational and synthetic maps of the Stokes parameter Q, obtained for the quasi-parallel case using equation (9). These maps are shown in Fig. 2. A good overall agreement between observations and numerical results is obtained, since the opposite bright arcs are well reproduced. Also, note that these arcs are slightly more elongated than the ones obtained for the non-turbulent ISM (Schneiter et al. 2015).

Figure 2.

Same as Fig. 1 but for maps of the Stokes parameter Q.

4.2 Polar-referenced angle and a statistical study

The polar-referenced angle distribution χ_r is given by

\begin{equation} \chi _r=\cos ^{-1}(\hat{r} \cdot {\hat{b}}_{\perp }), \end{equation}

(17)

where |$\hat{r}=(-x,y)/\sqrt{x^2+y^2}$| is the radial direction and |$\hat{b}_{\perp }=(-{\rm sin}(\chi _{\mathrm{B}}), {\rm cos}(\chi _{\mathrm{B}}))$| is the direction of the magnetic field perpendicular to the LoS, which is obtained through equation (16).

By using equation (17), polar-referenced angle maps are obtained for runs R1b and R2b, which are shown in Fig. 3. The blue colour indicates the region where the magnetic field is almost parallel to the radial direction, while red colour highlights the regions where it is mostly perpendicular.

Figure 3.

Synthetic polar-referenced angle maps obtained from numerical simulations of the SNR 1006 expanding into a turbulent medium and with adiabatic index γ = 5/3 (model R1b, panel (a)) and γ = 1.3 (model R2b, panel(b)) The axes are in pc and the linear colour scale is given in degrees.

Given that the polar-referenced angle technique is a useful tool for estimating the direction of the pre-shock ISM magnetic field (see Reynoso et al. 2013; Schneiter et al. 2015) we carried out a statistical study based on the analysis of the position-angle distribution in those regions where χ_r approaches zero.

The analysis was performed for both the numerical results and the observations in regions with synchrotron intensities larger than 10 per cent of the maximum, and where the condition χ_r ≤ 14° (the observational angle error) is satisfied. These criteria were chosen to assure a good signal-to-noise ratio of the sample. The results are shown in Fig. 4. To obtain the observational distribution, electric vectors were computed by combining the Q and U images with the miriad task IMPOL, and were later rotated by |$\frac{\pi }{2}$| to convert them into magnetic vectors. The observational curve has two maxima: a large one at ∼57° and a small one at ∼25°. Fig. 5 shows a fitting of the observational distribution by two Gaussian labeled as g1 and g2, whose parameters are given in Table 2. The mean value of Gaussian g1 is 58°, which is close to the inclination of the Galactic plane (60°) and in agreement with the maxima of the curves computed for the runs R1a and R2a (see Fig. 4).

Figure 4.

Comparison of the normalized distributions (with respect to their maxima) of the magnetic field position-angle as obtained from the observations and from the runs R1a and R2a (different γ and both corresponding to a SNR propagating into a uniform medium, i.e. with δ|$\boldsymbol B$| = δρ = 0; see Table 1), applying the polar-referenced angle selection.

Figure 5.

Two Gaussian fit of the observed position-angle distribution.

Table 2.

Two Gaussian fit of the observational position-angle distribution

Gaussian	Mean(°)	σ(°)
g1	58.0	9.0
g2	23.0	13.0

It is reasonable to expect that the resulting ambient magnetic field is parallel to the Galactic plane at the latitude of SN 1006 (Bocchino et al. 2011; Reynoso et al. 2013; Schneiter et al. 2015). The secondary peak of the observed position-angle distribution would indicate an extra magnetic field component, which is not considered in our simulations. In order to eliminate this secondary unmodelled component, we subtracted the Gaussian fit labeled g2 from the observed distribution. The resulting curve will be referred to as the modified observational distribution.

Fig. 6 compares the modified observational distribution with those resulting from the runs R1b and R2b. Clearly these numerical curves display now wider distributions compared with runs R1a and R2a, which is a consequence of including perturbations in |$\boldsymbol {B}$| and ρ. The maxima of numerical distributions remain close to 60°.

Figure 6.

Same as Fig. 4 but comparing the distribution obtained for observations (excluding the contribution of Gaussian g2) and runs R1b and R2b (with different γ and for a turbulent medium, i.e. including both δρ and δ|$\boldsymbol B$|⁠).

In Fig. 7, we explore the effect of selectively including the perturbations in B and/or ρ for γ = 1.3. The curve corresponding to the run R2d (δρ ≠ 0 and |$\delta \boldsymbol {B}$| =0) does not show any significant widening, in contrast with those corresponding to the runs R2b and R2c, both with |$\delta {\boldsymbol B}\,\ne 0$|⁠. In order to carry out a more quantitative study, a statistical analysis was performed, comparing the distributions of the position-angle obtained from observations (subtracting the contribution of Gaussian g2) and simulations. The results of this statistical analysis are summarized in Table 3. Comparing the values obtained for the standard deviation, skewness and kurtosis, we note that the models achieving a better fit with the observations are those that consider a perturbed magnetic field and a lower γ (runs R2b and R2c).

Figure 7.

Comparison of the position-angle distributions obtained for runs R2p, i.e. those considering γ = 1.3. In model R2a the medium is initialized as uniform at the beginning of the simulation, while in model R2b the SNR propagates into a turbulent medium (i.e. δB, δρ ≠ 0. In models R2c and R2d the medium is ‘partially’ turbulent, with δB ≠ 0, δρ = 0 and δρ ≠ 0, δB = 0, respectively.

Table 3.

Statistics of the position-angle distributions.

	Mean(°)	σ(°)	Skewness	Kurtosis
R1a	60	3.5	− 0.32	0.2
R1b	59	7.5	0.52	5.2
R2a	62	6.6	− 5.80	91.0
R2b	58	11.0	− 3.0	9.7
R2c	56	13.0	− 2.2	4.7
R2d	58	9.9	− 3.70	14.0
obs-g2	57	22.0	− 0.68	6.5

5 DISCUSSION AND CONCLUSIONS

A radio polarization study of SN 1006 was performed based on 3D MHD simulations, considering the expansion of the remnant into a turbulent ISM (including both magnetic field and density perturbations).

The inclusion of perturbations in magnetic field and density do not change the main conclusion of Schneiter et al. (2015), namely that the quasi-parallel case is the acceleration mechanism that better explains the observed morphology in the distribution of the Stokes parameter Q.

Based on the polar-referenced angle method, a study of the position angle distribution of the magnetic field was performed on both the observations and numerical results. This study reveals that the observational distribution is wide and has two maxima or components: a large one at 58°that coincides with the expected direction of the ambient magnetic field and a small component at 23°. The curves corresponding to runs R1a an R2a have a single peak and narrower distributions. The secondary peak of the observations could be due to local blow-outs produced during the expansion of parts of the main SNR shock front into low-density cavities, such as those explored by Yu et al. (2015). Since our simulations do not have this secondary component, we filtered it out from the observations and focused our comparison with the dominant component, as explained in the previous section.

The subsequent analysis was carried out by comparing the position-angle distributions obtained from numerical simulations with the observational distribution curve. A statistical study performed on these distributions reveals that models, which include a turbulent pre-shock magnetic field, successfully explain the observed polar-angle distribution. Furthermore, all of them share a maximum at a position angle of 60°, which is parallel to the Galactic plane and coincides with the expected direction of the Galactic magnetic field around SN1006 (Bocchino et al. 2011; Reynoso et al. 2013).

In summary, the statistical analysis based on the polar-referenced angle technique carried out in this work made it possible to recover the orientation of the ambient magnetic field. Based on previous observational and theoretical studies of SN1006 (Bocchino et al. 2011; Reynoso et al. 2013; Schneiter et al. 2015), we have only considered the quasi-parallel acceleration mechanism to estimate the synchrotron emission. However, it is important to mention that the same technique presented here can be applied with the quasi-perpendicular case. Finally, a good agreement between observations and numerical results is obtained if the simulations include magnetic field perturbations.

Acknowledgments

The authors thank the referee for reading this manuscript and for her/his useful comments. PFV, AE and FDC acknowledge financial support from DGAPA-PAPIIT (UNAM) grants IG100516, IN109715, IA103315, and IN117917. EMS, EMR, and DOG are members of the Carrera del Investigador Científico of CONICET, Argentina. MVS and AM-B are fellows of CONICET (Argentina) and CONACyT (México), respectively. EMS thanks Paulina Velázquez for her hospitality during his visit to Mexico City. We thank Enrique Palacios-Boneta (cómputo-ICN) for mantaining the Linux servers where our simulations were carried out.