MNRAS000,1–8(2002) Preprint23January2017 CompiledusingMNRASLATEXstylefilev3.0 A 3D MHD simulation of SN 1006: a polarized emission study for the turbulent case P. F. Vela´zquez1(cid:63), E. M. Schneiter2,3, E. M. Reynoso4, A. Esquivel1, F. De Colle1, J. C. Toledo-Roy5, D.O. Go´mez 4,6, M.V. Sieyra2, A. Moranchel-Basurto7 1Instituto de Ciencias Nucleares, Universidad Nacional Auto´noma de M´exico, Ciudad de M´exico, M´exico 2Instituto de Astronom´ıa Te´orica y Experimental, Universidad Nacional de Co´rdoba, C´ordoba, Argentina 3Departamento de Materiales y Tecnologia, UNC, Co´rdoba, Argentina 7 4CONICET-Universidad de Buenos Aires. Instituto de Astronom´ıa y F´ısica del Espacio, Suc. 28, CP: 1428, Buenos Aires, Argentina 1 5Centro de Ciencias de la Complejidad, Universidad Nacional Aut´onoma de M´exico, Ciudad de M´exico, M´exico. 0 6 Departamento de F´ısica, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires, Buenos Aires, Argentina 2 7Escuela Superior de F´ısica y Matema´tica, Instituto Polit´ecnico Nacional, Ciudad de M´exico, M´exico. n a J Accepted.Received 9 1 ABSTRACT ] Three dimensional magnetohydrodynamical simulations were carried out in order to E perform a new polarization study of the radio emission of the supernova remnant SN H 1006. These simulations consider that the remnant expands into a turbulent interstel- . lar medium (including both magnetic field and density perturbations). Based on the h p referenced-polarangletechnique,astatisticalstudywasdoneonobservationalandnu- - merical magnetic field position-angle distributions. Our results show that a turbulent o medium with an adiabatic index of 1.3 can reproduce the polarization properties of r t theSN1006remnant.Thisstatisticalstudyrevealsitselfasausefultoolforobtaining s the orientation of the ambient magnetic field, previous to be swept up by the main a supernova remnant shock. [ 1 Key words: MHD–radiation mechanisms: general – methods : numerical – super- v novae: individual: SN1006 –ISM: supernova remnants 9 8 5 5 0 1 INTRODUCTION in turn affects the radio emission and is imprinted in the . measured Stokes parameters. 1 SN 1006 is a young supernova remnant (SNR) that has Jun&Jones(1999)performed2DMHDsimulationsof 0 gainedinterestduetothelargeobservationaldataavailable the evolution of a young Type Ia SNR interacting with an 7 inawidewavelengthrange,thereforeservingasanideallab- interstellarcloud.Theyfoundthattheinteractionproduced 1 oratorywhentryingtounderstandtheobservedmorphology : RT instabilities causing amplification of the magnetic field, v in each frequency. This remnant is classified as part of the andsubsequentlychangingthesynchrotronbrightness.This i bilateral SNR group, whose main characteristic is the pres- X work emphasized the importance of including a turbulent ence of two bright and opposite arcs in radio-frequencies. medium when trying to match observations. In this direc- r The morphology and emission of this type of SNRs is a tion,Balsaraetal.(2001)performedmorerealistic3Dsim- largelydeterminedbytheenvironmentinwhichtheyevolve, ulationsofaSNRevolvinginaturbulentbackground.They sometimeshavinganirregularorclumpybackgroundand/or found that the variability of the emission during the early with the presence of density and magnetic fields gradients. phase of SNR evolution gives information about the turbu- The interstellar medium is known to be turbulent, as lent environment. At the same time, that the interaction inferred from observational data (Lee & Jokipii 1976) and with such a turbulent environment generates an enhanced local information on the interplanetary medium (Jokipii turbulent post-shock region which has observational conse- et al. 1969). Such turbulent medium, with a Kolmogorov- quences,asitaffectstheaccelerationofrelativisticparticles. likepowerspectrumthatspansoveralargerangeofspatial Guo et al. (2012) employed 2D MHD simulations to scales (Guo et al. 2012), gives rise to Rayleigh-Taylor (RT) study the magnetic field amplification mechanisms during instabilities during the fast expansion of the SNR, which 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 (cid:63) Email:[email protected] related to the shock amplification, and the other, to the (cid:13)c 2002TheAuthors 2 Vel´azquez et al. RT-instabilities at the contact discontinuity between the shocked ejecta and the shocked interstellar medium. Fang&Zhang(2012)carriedoutastudysimilartothat (a) of Guo et al. (2012) but considering different adiabatic in- dexes to account for both the diffusive shock acceleration and the escape of the accelerated particles from the shock. Theyfoundthatasmallereffectiveadiabaticindexproduces a larger magnetic field enhancement. They extended their studyconsideringascenariowheretheexpandingejectain- teracts 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 sub- sequentstudy,Fangetal.(2014)performed3DMHDsimu- lationsconsideringaturbulentmediumwithaKolmogorov- like power spectrum with parameters resembling the SNR (b) RX J0852.0-4622. Their X-ray and γ-ray synthetic maps showedtheexpectedrippledmorphologywithabrokencir- cular shape along the shell, nicely reproducing the observa- tions. In general, the diffusive shock acceleration is driven by magneticfieldsstrongenoughtogeneratenon-thermalemis- sion. Therefore, any structure or dynamic interaction that resultsintheamplificationofthemagneticfieldhasimpor- tant consequences on the observed brightness. According to the non-thermal radio-emission morphol- ogy, 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 (c) (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 in- terstellar medium. In Schneiter et al. (2010) and Schneiter et al. (2015) we made use of this assumption and presented studiesofthesystemin2Dand3D,respectively.Inthelat- terworkwefurtherintroducedsyntheticmapsoftheStokes Q and U parameters wich allowed a better and more direct comparison with the observations (Reynoso et al. 2013). Recently, West et al. (2016b, a) (see also, West et al. 2016a,b)studiedradioemissionforasampleofGalacticax- isymmetric SNRs, including SN 1006. They concluded that thequasi-perpendicularaccelerationmechanismcansucess- Figure 1. Comparison of the linearly polarized intensity at 1.4 fully fit the observed morphology of the remnants in their GHzbetweentheobservations(panel(a))andsyntheticmapsofa sample, although SN 1006 seems to be the exception, since SNRexpandingintoaturbulentmediumwithadiabaticindexγ= the quasi-parallel mechanism agrees better with observa- 5/3(modelR1b,panel(b))andγ=1.3(modelR2b,panel(c)) tions. InthepresentworkweusethesameSNRinitialization as in the two previous works of Schneiter et al. (2010) and presentedinsection§4.Thediscussionandfinalconclusions Schneiteretal.(2015)butrelaxingtheassumptionofaho- are given in section §5. mogeneousISMbyconsideringthattheSNRevolvesintoa turbulent medium as it was shown in Yu et al. (2015) (see also Fang et al. 2014). 2 NUMERICAL MODEL The goal of this work is to propose an alternative methodtodeterminethepositionangleoftheambientmag- The numerical model employed is basically the same as in netic field, based on the polar-referenced angle method em- Schneiteretal.(2010,2015),withtheadditionalassumption ployedbyReynosoetal.(2013)andSchneiteretal.(2015). ofaturbulentbackground,whereboththedensityandmag- The numerical model is presented in section §2. For netic fluctuations have a 3D Kolmogorov-like power spec- completeness we recall the basic model employed by trum, as will be explained further down. Schneiter et al. (2015) in §2.1 and in section §2.2 we ex- To simulate the evolution of SN 1006 we employed the plain how the turbulence is introduced. Section §3 explains Mezcal code (De Colle & Raga 2006; De Colle et al. 2008, howsyntheticradiomapswereobtainedandtheresultsare 2012).ThiscodesolvesthefullsetofidealMHDequationsin MNRAS000,1–8(2002) A 3D MHD study of the SN1006: the turbulent case 3 a Cartesian geometry with an adaptive mesh, and includes a cooling function to account for radiative losses (De Colle (a) & Raga 2006). The computational domain is a cube of 24 pc per side, whichwewilldenoteas(x,y,z),andisdiscretizedonafive level binary grid with a maximum resolution of 4.7×10−2 pc. All the outer boundaries were set to outflow condition. 2.1 SNR initial conditions In what follows, we describe the physical parameters and setupusedtosimulatetheSNR,whicharethesameasthose presented in Schneiter et al. (2015). A supernova explosion is initialized by the deposition of E =2.05×1051 erg in a radius of R =0.65 pc located 0 0 (b) at the centre of the computational domain. The energy is distributedsuchthat95%ofitiskineticandtheremaining 5% is thermal. Theejectedmasswasdistributedintwoparts:aninner homogeneous sphere of radius r containing 4/7ths of the c total mass (M = 1.4M ) with a density ρ , and an outer ∗ (cid:12) c shell containing the remaining 3/7ths of the mass following a power law (ρ ∝ r−7) as in Jun & Norman (1996a). The velocityhasaincreasinglinearprofilewithr,whichreaches a value of v at r = R . The parameters ρ , r , and v 0 0 c c 0 are functions of E , M , and R , and were computed using 0 ∗ 0 equations (1)-(3) of Jun & Norman (1996a). (c) 2.2 The turbulent background Jun&Jones(1999)suggestedthattheturbulentstructures of the brightest radio emission correlate well with the mag- neticfieldingeneral.Forthisreasonweperformednewsim- ulationswithamorerealisticmagnetizedturbulentinterstel- lar medium. In what follows, we will make use of equations (10) to (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)(seealsoYuetal.2015;Fangetal.2014).Thepower spectrum follows: Figure2.SameasFigure1butformapsoftheStokesparameter Q. 1 P ∝ . (1) 1+(kL )11/3 c with To generate the turbulence we assume a coherence ulelnagtitohnsofarLecc=om3ppucteadnwdiathwNavmen=um90b3erwkav=e m2Lπod.eOs,uarnsdimL- A2(kn)=σB2 1+(k∆VLn)11/3 (cid:88)Nm(cid:20)1+(k∆VLn)11/3(cid:21)−1, (4) variesbetweenLmin =∆xandLmax =Lsim,being∆xand n c n=1 n c L the cell and the computational domain sizes, respec- sim whereσ2 isthewavevarianceofthemagneticfield,andthe tively. B normalization factor, ∆V , is given by: The initial magnetic field is given by: n ∆V =4πk2∆k . (5) n n n B(x,y,z)=B +δB, (2) 0 Both systems (x(cid:48),y(cid:48),z(cid:48)) and (x,y,z) are related by: whereB =B (x)yˆandthemagneticfieldperturbationis: 0 0 x(cid:48) cosθ cosφ cosθ sinφ −sinθ x n n n n n δB=(cid:60)(cid:20)(cid:88)NmA(kn)(cid:0)cosαnxˆ(cid:48)+isinαnyˆ(cid:48)(cid:1)exp(iknzn(cid:48))(cid:21), yz(cid:48)(cid:48)=sin−θsnincoφsnφn sincθonssφinnφn cos0θn yz n=1 (3) where θ and φ represent the direction of propagation of n n MNRAS000,1–8(2002) 4 Vel´azquez et al. run γ δB δρ (a) 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 Table 1. Runs carried out in this paper with their correspond- ing hypothesis: adiabatic index, turbulent magnetic field and/or turbulentdensitydistribution. of a turbulent background with a Kolgomorov-like power spectrum for both density and magnetic field. These runs (b) 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”indicatesthepresenceorabsenceofturbulencein the background magnetic field and density: p = a indicates no turbulence in either variable; p = b includes turbulence inboth;p=conlyconsidersaturbulentB;andp=donly considers a turbulent ρ. Synthetic radio emission maps were obtained from the numericalresults.Inordertocomparethemwiththeobser- vations,thecomputationaldomain(denotedasxyzsystem), theywererotatedwithrespecttothe“image”system(x y z i i i system). At the beginning both systems coincide. The x y i i plane of the “image” system is the plane of the sky, being yˆ the direction towards the “North”, −xˆ the direction to- i i Figure 3. Synthetic polar-referenced angle maps obtained from wards the “East”, and the line of sight (hereafter LoS) the numericalsimulationsoftheSNR1006expandingintoaturbulent zˆdirection. Then, synthetic maps were obtained after per- mediumandwithadiabaticindexγ=5/3(modelR1b,panel(a)) formingthreerotationsonthecomputationaldomain.First, and γ = 1.3 (model R2b, panel(b)) The axes are in pc and the linearcolorscaleisgivenindegrees. arotationaroundyˆi wasapplied(the“North”)withanan- gle ϕ . Second, a rotation in ϕ was carried out around yi xi thexˆ -direction.Finally,thecomputationaldomainwasro- i thewavemodenwiththewavenumberkn andpolarization tated by ϕzi around the zˆi-direction (the line of sight). To αn. get the best agreement with the observations, after several For the density fluctuation we employed the same log- tests these angles were set as −15◦, −30◦, and 60◦ for ϕ , xi normal distribution as in Giacalone & Jokipii (2007): ϕ , and ϕ , respectively . In this way, the projection of yi zi the unperturbed magnetic field B was tilted by 60◦ with 0 respect to the yˆ direction. n(x,y,z)=n exp(f +δf), (6) i 0 0 where f is a constant and δf is the density perturbation 0 withawavevarianceσ2.Bothδf andδBhavethesame3D 3.1 Synchrotron emissivity d Kolmogorov-like power spectral index. Foreachpoint(x ,y ,z )oftheSNRwecanobtainthesyn- The turbulent environment temperature and number i i i chrotron specific intensity as (see C´ecere et al. 2016, and density were set to T =104K and n =5×10−2cm−3, re- 0 0 references therein): spectively.Thewavevariancesformagneticfieldanddensity were set as σ2 =0.25B2 and σ2 =0.4n2. B 0 d 0 As mentioned above, one of the purposes of this work j(x ,y ,z ,ν)∝Kρ v4α Bα+1ν−α, (7) i i i ⊥ is to be able to study the effect of a perturbation on some observedparameters.Tocomparewiththeobservations,the where ν is the observed frequency, B⊥ is the component of synthetic maps were performed for an integration time cor- themagneticfieldperpendiculartotheLoS,αisthespectral responding to 1000 yr (the age of SN 1006). index -which was set to 0.6 for this object- and ρ and v are the density and velocity of the gas, respectively. The coeffi- cientK includestheobliquitydependence,beingeitherpro- portional to sin2Θ for the quasi-perpendicular case or to 3 SYNTHETIC EMISSION MAPS Bs cos2Θ forthequasi-parallelcase.TheangleΘ isthean- Bs Bs InTable1welisttherunsthatwerecarriedout.Withthese gle between the shock normal and the post-shock magnetic runs we want to analyze the effects of both the reduction field(Fulbright&Reynolds1990).ForthecaseofSN1006, of the value of the adiabatic index γ and/or the inclusion Bocchinoetal.(2011)andReynosoetal.(2013)showedthat MNRAS000,1–8(2002) A 3D MHD study of the SN1006: the turbulent case 5 The position angle of the local magnetic field φ (x ,y ,z ) B i i i is known from the simulations, and the position angle of the local electric field φ(x ,y ,z ) can be obtained from i i i φ (x ,y ,z ) by applying a π rotation and correcting for B i i i 2 Faraday rotation, as: π φ(x ,y ,z )=φ (x ,y ,z )− +∆χ . (12) i i i B i i i 2 F The Faraday correction term ∆χ is given by: F ∆χ =RM λ2, (13) F where RM is the rotation measure (in units of rad m−2) and λ is the wavelength of the observations (given in m). In the present study, only the external RM arising in the foreground ISM is considered. Bandiera & Petruk (2016) explored the effects of internal RM on the polarized syn- chrotronemissionofSNRs.ForthecaseofSN1006,wehave Figure 4.Comparisonofthenormalizeddistributions(withre- carriedoutaroughestimationoftheinternalRM,obtaining spect to their maxima) of the magnetic field position-angle as a value close to 10% of the reported one by Reynoso et al. obtained from the observations and from the runs R1a and R2a (2013)whoobtainedaRMof12rad m−2 forthisremnant. (different γ and both corresponding to a SNR propagating into The linearly polarized intensity is computed as: auniformmedium,i.e.withδB=δρ=0;seetable1),applying thepolar-referencedangleselection. (cid:112) I (x ,y ,ν)= Q(x ,y ,ν)2+U(x ,y ,ν)2 (14) P i i i i i i the observed morphology of this remnant can be explained and the maps of the polarization angle distribution were byconsideringthequasi-parallelcase.Schneiteretal.(2015) computed as follows: carriedoutapolarisationstudyofthisremnant,basedon3D MHD simulations, finding that only the quasi-parallel case 1 χ(x ,y )= tan−1(U(x ,y ,ν)/Q(x ,y ,ν)) (15) can succesfully reproduced the observed morpholoy of the i i 2 i i i i Stokes parameter Q (see Figure 3 of Schneiter et al. 2015). Similarly,thedistributionofthemagneticfieldorienta- For these reasons only the quasi-parallel case was consid- tion can be calculated as: eredinthiswork.Thesyntheticsynchrotronemissionmaps are obtained by integrating j(x ,y ,z ,ν) along the LoS or i i i 1 zi-axis, i.e.: χB(xi,yi)= 2tan−1(UB(xi,yi,ν)/QB(xi,yi,ν)) (16) (cid:90) where UB(xi,yi,ν) and QB(xi,yi,ν)) are calculated with I(xi,yi,ν)= j(xi,yi,zi,ν)dzi, (8) the equations (9) and (10) and replacing φ(xi,yi,ν) with LoS φ (x ,y ,ν). B i i being x and y the coordinates in the plane of the sky. i i 4 RESULTS 3.2 Stokes parameters and position angle distribution maps 4.1 Synthetic polarization maps With the purpose of comparing the numerical results with Figure 1 shows a comparison between the observed map of the observations, synthetic maps of the Stokes parameters the linearly polarized intensity of SN 1006 (upper panel), QandU werecalculatedasfollows(Clarkeetal.1989;Jun and the synthetic polarization maps obtained for runs with & Norman 1996b; Schneiter et al. 2015): γ =5/3(runR1b,middlepanel)andγ =1.3(runR2b,bot- tom panel) by using Equation (14). The observed map was (cid:90) constructedbycombiningdataobtainedwiththeATCAand Q(xi,yi,ν)= f0j(xi,yi,zi,ν)cos[2φ(xi,yi,zi)]dzi, theVLAat1.4GHz,asexplainedinReynosoetal.(2013). LoS The polarized emission image, as well as the distribution of (9) the Q and U parameters, were convolved to a 15(cid:48)(cid:48) beam. The two models used to construct the synthetic maps con- (cid:90) U(x ,y ,ν)= f j(x ,y ,z ,ν)sin[2φ(x ,y ,z )]dz , sider that the SNR is expanding into a turbulent medium i i 0 i i i i i i i LoS withperturbationsinbothdensityandmagneticfield.These (10) synthetic maps are shown in arbitrary units. Two opposite where φ(x ,y ,z ) is the position angle of the local electric i i i bright,noisyarcsareobservedinbothsyntheticmaps,hav- field in the plane of the sky and f is the degree of linear 0 ing a striking resemblance to the observations. Reducing polarization, which is a function of the spectral index α: the value of γ produces a 5% smaller expansion radius of α+1 theremnantasmeasuredonthelinearlypolarizedintensity f = . (11) 0 α+5/3 maps. MNRAS000,1–8(2002) 6 Vel´azquez et al. Gaussian mean(◦) σ(◦) g1 58.0 9.0 g2 23.0 13.0 Table 2. Two Gaussian fit of the observational position-angle distribution field perpendicular to the LoS, which is obtained through Equation (16). By using Equation (17), polar-referenced angle maps areobtainedforrunsR1bandR2b,whichareshowninFig- ure3.Thebluecolorindicatestheregionwherethemagnetic fieldisalmostparalleltotheradialdirection,whileredcolor highlights the regions where it is mostly perpendicular. Giventhatthepolar-referencedangletechniqueisause- Figure 5. Two Gaussian fit of the observed position-angle dis- tribution. ful tool for estimating the direction of the pre-shock ISM magneticfield(seeReynosoetal.2013;Schneiteretal.2015) we carried out a statistical study based on the analysis of theposition-angledistributioninthoseregionswhereχ ap- r proaches zero. The analysis was performed for both the numerical re- sults and the observations in regions with synchrotron in- tensities larger than 10% of the maximum, and where the condition χ (cid:54) 14◦ (the observational angle error) is satis- r fied. These criteria were chosen to assure a good signal to noise ratio of the sample. The results are shown in Figure 4. To obtain the observational distribution, electric vectors were computed by combining the Q and U images with the miriadtaskIMPOL,andwerelaterrotatedby π toconvert 2 themintomagneticvectors.Theobservationalcurvehastwo maxima:alargeoneat∼57◦ andasmalloneat∼25◦.Fig- ure5showsafittingoftheobservationaldistributionbytwo 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 agreementwiththemaximaofthecurvescomputedforthe Figure 6.SameasFigure4butcomparingthedistributionob- tained for observations (excluding the contribution of Gaussian runs R1a and R2a (see Figure 4). g2)andrunsR1bandR2b(withdifferentγ andforaturbulent It is reasonable to expect that the resulting ambient medium,i.e.includingbothδρandδB). magnetic field is parallel to the Galactic plane at the lati- tudeofSN1006(Bocchinoetal.2011;Reynosoetal.2013; Schneiter et al. 2015). The secondary peak of the observed Following the analysis developed by Schneiter et al. position-angledistributionwouldindicateanextramagnetic (2015), we compare the observational and synthetic maps fieldcomponent,whichisnotconsideredinoursimulations. of the Stokes parameter Q, obtained for the quasi-parallel Inordertoeliminatethissecondaryunmodelledcomponent, case using Equation 9. These maps are shown in Figure 2. wesubtractedtheGaussianfitlabeledg2fromtheobserved Agoodoverallagreementbetweenobservationsandnumer- distribution. The resulting curve will be referred to as the ical results is obtained, since the opposite bright arcs are modified observational distribution. wellreproduced.Also,notethatthesearcsareslightlymore Figure 6 compares the modified observational distribu- elongatedthantheonesobtainedforthenon-turbulentISM tionwiththoseresultingfromtherunsR1bandR2b.Clearly (Schneiter et al. 2015). thesenumericalcurvesdisplaynowwiderdistributionscom- paredwithrunsR1aandR2a,whichisaconsequenceofin- cludingperturbationsinBandρ.Themaximaofnumerical 4.2 Polar-referenced angle and a statistical study distributions remain close to 60◦. In Figure 7 we explore the effect of selectively includ- The polar-referenced angle distribution χ is given by: r ing the perturbations in B and/or ρ for γ = 1.3. The curve χ =cos−1(rˆ·ˆb ), (17) corresponding to the run R2d (δρ (cid:54)= 0 and δB= 0) does r ⊥ not show any significant widening, in contrast with those (cid:112) where rˆ = (−x,y)/ x2+y2 is the radial direction and corresponding to the runs R2b and R2c, both with δB(cid:54)=0. ˆb = (−sin(χ ),cos(χ )) is the direction of the magnetic Inordertocarryoutamorequantitativestudy,astatistical ⊥ B B MNRAS000,1–8(2002) A 3D MHD study of the SN1006: the turbulent case 7 formedonboththeobservationsandnumericalresults.This mean(◦) σ(◦) skewness kurtosis studyrevealsthattheobservationaldistributioniswideand has two maxima or components: a large one at 58◦that co- R1a 60. 3.5 −0.32 0.2 incideswiththeexpecteddirectionoftheambientmagnetic R1b 59. 7.5 0.52 5.2 field,andasmallcomponentat23◦.Thecurvescorrespond- R2a 62. 6.6 −5.80 91.0 R2b 58. 11.0 −3.0 9.7 ing to runs R1a an R2a have a single peak and narrower R2c 56. 13.0 −2.2 4.7 distributions.Thesecondarypeakoftheobservationscould R2d 58. 9.9 −3.70 14.0 beduetolocalblow-outsproducedduringtheexpansionof obs-g2 57. 22.0 −0.68 6.5 partsofthemainSNRshockfrontintolow-densitycavities, such as those explored by Yu et al. (2015). Since our sim- Table 3.Statisticsoftheposition-angledistributions ulations do not have this secondary component we filtered it our from the observations and focused our comparison withthedominantcomponent,asexplainedintheprevious section. The subsequent analysis was carried out by compar- ingtheposition-angledistributionsobtainedfromnumerical simulationswiththeobservationaldistributioncurve.Asta- tistical 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 angleof60◦,whichisparalleltotheGalacticPlaneandco- incideswiththeexpecteddirectionoftheGalacticmagnetic field around SN1006 (Bocchino et al. 2011; Reynoso et al. 2013). Insummary,thestatisticalanalysisbasedonthepolar- referenced angle technique carried out in this work made it possible to recover the orientation of the ambient mag- netic field. Based on previous observational and theoreti- cal studies of SN1006 (Bocchino et al. 2011; Reynoso et al. 2013; Schneiter et al. 2015), we have only considered the Figure 7. Comparison of the position-angle distributions ob- quasi-parallel acceleration mechanism to estimate the syn- tainedforrunsR2p,i.e.thoseconsideringγ=1.3.InmodelR2a chrotronemission.However,itisimportanttomentionthat themediumisinitializedasuniformatthebeginningofthesim- the same technique presented here can be applied with the ulation,whileinmodelR2btheSNRpropagatesintoaturbulent quasi-perpendicularcase.Finally,agoodagreementbetween medium(i.e.,δB,δρ(cid:54)=0.InmodelsR2candR2dthemediumis observations and numerical results is obtained if the simu- “partially” turbulent, with δB (cid:54)= 0,δρ = 0 and δρ (cid:54)= 0,δB = 0 lations include magnetic field perturbations. respectively. analysis was performed, comparing the distributions of the ACKNOWLEDGMENTS position-angle obtained from observations (subtracting the contributionofGaussiang2)andsimulations.Theresultsof The authors thank the referee for reading this manuscript thisstatisticalanalysisaresummarizedinTable3.Compar- and for her/his useful comments. PFV, AE, and FDC ingthevaluesobtainedforthestandarddeviation,skewness acknowledge finantial support from DGAPA-PAPIIT andkurtosis,wenotethatthemodelsachievingabetterfit (UNAM) grants IG100516, IN109715, IA103315, and with the observations are those that consider a perturbed IN117917. EMS, EMR, and DOG are members of the Car- magnetic field, and a lower γ (runs R2b and R2c). rera del Investigador Cient´ıfico of CONICET, Argentina. MVS and AM-B are fellows of CONICET (Argentina) and CONACyT (M´exico), respectively. EMS thanks Paulina Vela´zquezforherhospitalityduringhisvisittoMexicoCity. 5 DISCUSSION AND CONCLUSIONS WethankEnriquePalacios-Boneta(co´mputo-ICN)forman- tainingtheLinuxserverswhereoursimulationswerecarried AradiopolarizationstudyofSN1006wasperformedbased out. on 3D MHD simulations, considering the expansion of the remnant into a turbulent interstellar medium (including both magnetic field and density perturbations). The inclusion of perturbations in magnetic field and REFERENCES densitydonotchangethemainconclusionofSchneiteretal. BalsaraD.,BenjaminR.A.,CoxD.P.,2001,ApJ.,563,800 (2015), namely that the quasi-parallel case is the accelera- BandieraR.,PetrukO.,2016,MNRAS,459,178 tion mechanism that better explains the observedmorphol- BocchinoF.,OrlandoS.,MiceliM.,PetrukO.,2011,A&A,531, ogy in the distribution of the Stokes parameter Q. A129 Basedonthepolar-referencedanglemethod,astudyof Cassam-Chena¨ı G., Hughes J. P., Reynoso E. M., Badenes C., thepositionangledistributionofthemagneticfieldwasper- MoffettD.,2008,ApJ,680,1180 MNRAS000,1–8(2002) 8 Vel´azquez et al. C´ecereM.,Vel´azquezP.F.,AraudoA.T.,DeColleF.,Esquivel A., Carrasco-Gonz´alez C., Rodr´ıguez L. F., 2016, ApJ, 816, 64 ClarkeD.A.,BurnsJ.O.,NormanM.L.,1989,ApJ,342,700 DeColleF.,RagaA.C.,2006,A&A,449,1061 DeColleF.,RagaA.C.,EsquivelA.,2008,ApJ.,689,302 DeColleF.,GranotJ.,L´opez-C´amaraD.,Ramirez-RuizE.,2012, ApJ.,746,122 FangJ.,ZhangL.,2012,MNRAS,424,2811 FangJ.,YuH.,ZhangL.,2014,MNRAS,445,2484 FulbrightM.S.,ReynoldsS.P.,1990,ApJ,357,591 GiacaloneJ.,JokipiiJ.R.,2007,ApJL,663,L41 GuoF.,LiS.,LiH.,GiacaloneJ.,JokipiiJ.R.,LiD.,2012,ApJ., 747,98 JokipiiJ.R.,1987,ApJ,pp842–846 JokipiiJ.R.,LercheI.,SchommerR.A.,1969,ApJL,157,L119 JunB.-I.,JonesT.W.,1999,ApJ,511,774 JunB.-I.,NormanM.L.,1996a,ApJ,465,800 JunB.-I.,NormanM.L.,1996b,ApJ,472,245 KestevenM.J.,CaswellJ.L.,1987,A&A,183,118 LeeL.C.,JokipiiJ.R.,1976,ApJ,206,735 ReynosoE.M.,HughesJ.P.,MoffettD.A.,2013,AJ,145,104 Schneiter E. M., Vel´azquez P. F., Reynoso E. M., de Colle F., 2010,MNRAS,408,430 SchneiterE.M.,Vel´azquezP.F.,ReynosoE.M.,EsquivelA.,De ColleF.,2015,MNRAS,449,88 West J. L., Safi-Harb S., Ferrand G., 2016a, preprint, (arXiv:1609.02887) West J. L., Safi-Harb S., Jaffe T., Kothes R., Landecker T. L., FosterT.,2016b,A&˜A,587,A148 YuH.,FangJ.,ZhangP.F.,ZhangL.,2015,A&A,579,A35 ThispaperhasbeentypesetfromaTEX/LATEXfilepreparedby theauthor. MNRAS000,1–8(2002)