Radio polarization maps of shell-type SNRs II. Sedov models with evolution of turbulent magnetic field

TL;DR

This paper presents a detailed model for polarized radio emission from shell-type supernova remnants during their Sedov phase, incorporating turbulent magnetic fields, particle acceleration, and Faraday rotation to produce realistic polarization maps.

Contribution

It introduces a comprehensive 3D modeling approach for SNR polarization maps, including magnetic turbulence, particle acceleration, and Faraday effects, advancing the simulation of observational data.

Findings

Generated detailed polarization maps for SNRs during Sedov evolution.

Analyzed effects of magnetic turbulence and Faraday rotation on polarization.

Provided insights into the magnetic field structure within SNRs.

Abstract

Polarized radio emission has been mapped with great detail in several Galactic supernova remnants (SNRs), but has not yet been exploited to the extent it deserves. We have developed a method to model maps of the Stokes parameters for shell-like SNRs during their Sedov evolution phase. At first, 3-dimensional structure of a SNR has been computed, by modeling the distribution of the magnetohydrodynamic parameters and of the accelerated particles. The generation and dissipation of the turbulent component of magnetic field everywhere in SNR are also considered taking into account its interaction with accelerated particles. Then, in order to model the emission, we have used a generalization of the classical synchrotron theory, valid for the case in which the magnetic field has ordered and disordered components. Finally, 2-dimensional projected maps have been derived, for different orientations of SNR and of interstellar magnetic field with respect to the observer. An important effect to consider is the Faraday rotation of the polarization planes inside the SNR interior. In this paper we present details of the model, and describe general properties of the images.

Tables

Table 1. Table 1: Sets of parameters for the numerical simulations. The parameter β n subscript 𝛽 n \beta_{\mathrm{n}} is defined by Eq. ( 113 ). Set A is the reference model; while the other models differ from the case A only for the values shown in the table. ‘iso’ means isotropic injection (i.e. independent of the shock obliquity), while ‘qpar’ means quasiparallel injection (i.e. mostly at the parallel shocks).

Model	${\beta }_n$	injection	${\Omega }_{o\parallel }$	${\left(\frac{\delta B}{B}\right)}_{s\parallel }$
A	1	iso	0.4	0.3
B	0.2
C	0.5
D	2
E	5
F		qpar
G			10
H				1
I			10	1

Equations

$$u(r,t)=u_{\mathrm{s}}(t)\bar{u}(\bar{r}),\bar{r}=r/R(t)$$

$${\sigma }_{\mathrm{B}}={\left({\cos }^2{\Theta }_{\mathrm{o}}+{\sigma }^2{\sin }^2{\Theta }_{\mathrm{o}}\right)}^{1/2},$$

$$\left\{\begin{array}[]{lcl}{\cal I}&=&{\cal I}^{\prime}\\ {\cal Q}&=&{\cal Q}^{\prime}\cos 2\chi-{\cal U}^{\prime}\sin 2\chi={\cal Q}^{\prime}\cos 2\chi\\ {\cal U}&=&{\cal Q}^{\prime}\sin 2\chi+{\cal U}^{\prime}\cos 2\chi={\cal Q}^{\prime}\sin 2\chi\\ {\cal V}&=&{\cal V}^{\prime}=0\end{array}\right.$$

$$\left\{\begin{array}[]{lcl}\cos 2\chi&=&\displaystyle\frac{B_{\mathrm{y}}^{2}-B_{\mathrm{x}}^{2}}{B_{\mathrm{x}}^{2}+B_{\mathrm{y}}^{2}},\\ \\ \sin 2\chi&=&-\displaystyle\frac{2B_{\mathrm{x}}B_{\mathrm{y}}}{B_{\mathrm{x}}^{2}+B_{\mathrm{y}}^{2}}.\end{array}\right.$$

$$\beta (l,z_1)=RM(l,z_1){\lambda }^2$$

$$RM(l,z_1)=\frac{e^3}{2\pi m_{\mathrm{e}}^2{c}^4}{\int }_l^{z_1}{n}_{\mathrm{e}}{B}_{\mathrm{z}}dz.$$

$$\left\{\begin{array}[]{lcl}I&=&\int_{z_{2}}^{z_{1}}{\cal I}^{\prime}(l)\,dl\\ \\ Q&=&\int_{z_{2}}^{z_{1}}{\cal Q}^{\prime}(l)\cos 2\left(\chi(l)+\beta(l,z_{1})\right)dl\\ \\ U&=&\int_{z_{2}}^{z_{1}}{\cal Q}^{\prime}(l)\sin 2\left(\chi(l)+\beta(l,z_{1})\right)dl\end{array}\right.$$

$$\Pi =\frac{\sqrt{Q^2+U^2}}{I}$$

$$\Psi =\frac{1}{2}\arctan \left(\frac{U}{Q}\right).$$

$${\mathcal{B}}_{\mathrm{x}}$$

$${\mathcal{B}}_{\mathrm{y}}$$

$$⟨\Pi ⟩=\frac{1}{N}\sum _i^N{\Pi }_i.$$

$${\mathcal{I}}_{\mathrm{single}}^{\prime }(\nu ,p)=\frac{\sqrt{3}{e}^3}{4\pi m{c}^2}{B}_{\mathrm{xy}}F\left(\frac{\nu }{{\nu }_{\mathrm{c}}}\right),$$

$${\mathcal{Q}}_{\mathrm{single}}^{\prime }(\nu ,p)=\frac{\sqrt{3}{e}^3}{4\pi m{c}^2}{B}_{\mathrm{xy}}G\left(\frac{\nu }{{\nu }_{\mathrm{c}}}\right)$$

$${\mathcal{I}}^{\prime }$$

$${\mathcal{Q}}^{\prime }$$

$$\mathcal{C}=\frac{\sqrt{3}{e}^3}{4\pi m{c}^2}{\left(\frac{\nu }{c_1}\right)}^{-(s-1)/2}\Gamma \left(\frac{s}{4}+\frac{7}{12}\right)\Gamma \left(\frac{s}{4}-\frac{1}{12}\right),$$

$${\Pi }_{\max }=\frac{s+1}{s+7/3}.$$

$${\mathcal{I}}^{\prime }$$

$${\mathcal{Q}}^{\prime }$$

$$\Pi$$

$$\frac{\partial E_{\mathrm{w}}}{\partial t}+\frac{\partial F_{\mathrm{w}}}{\partial r}=u\frac{\partial P_{\mathrm{w}}}{\partial r}+{\sigma }_{\mathrm{w}}{P}_{\mathrm{w}}-{\Gamma }_{\mathrm{w}}{P}_{\mathrm{w}}$$

$$F_{\mathrm{w}}=\frac{\delta B^2}{4\pi }u+\frac{\delta B^2}{8\pi }u.$$

$$v_{\mathrm{A}}=\frac{B}{\sqrt{4\pi \rho }}.$$

$$\frac{\partial P_{\mathrm{w}}}{\partial t}+u\frac{\partial P_{\mathrm{w}}}{\partial r}+P_{\mathrm{w}}\frac{3}{2}\frac{\partial u}{\partial r}=\frac{1}{2}\left({\sigma }_{\mathrm{w}}{P}_{\mathrm{w}}-{\Gamma }_{\mathrm{w}}{P}_{\mathrm{w}}\right).$$

$${\sigma }_{\mathrm{w}}(k)=\frac{4\pi }{3}\frac{v_{\mathrm{A}}\left|\hat{b}\cdot \hat{n}\right|}{kW(k)}{p}^{4}v\frac{\partial f}{\partial r}$$

$${\sigma }_{\mathrm{w}}{P}_{\mathrm{w}}=\int {\sigma }_{\mathrm{w}}(k)W(k)dk$$

$$P_{\mathrm{c}}=\frac{4\pi }{3}\int p^{4}vf\frac{dp}{p},$$

$${\sigma }_{\mathrm{w}}{P}_{\mathrm{w}}=v_{\mathrm{A}}\cos \Theta \frac{\partial P_{\mathrm{c}}}{\partial r}.$$

$${\Gamma }_{\mathrm{NL}}(k)W(k)=c_{\mathrm{K}}{v}_{\mathrm{A}}k\frac{E_{\mathrm{w}}(>k)}{B^2/4\pi }W(k)$$

$$\int {\Gamma }_{\mathrm{NL}}(k)W(k)dk=\frac{c_{\mathrm{K}}{v}_{\mathrm{A}}}{B^2/4\pi }\int k{E}_{\mathrm{w}}(>k)W(k)dk.$$

$$1\equiv E_{\mathrm{w}}{\left(\int W(k)dk\right)}^{-1},$$

Full text

††pubyear: 2016††pagerange: Radio polarization maps of shell-type SNRs II. Sedov models with evolution of turbulent magnetic field–A

1 Introduction

Supernova remnants (SNRs) are very important objects for High-Energy Astrophysics. They provide experimental information, in particular, on the processes related to magnetic turbulence, cosmic rays (CRs), their interactions and influence on the fluid dynamics and the shock properties. Interstellar magnetic field (ISMF, of order of few ${\,\rm\mu G}$) considerably affects the propagation of all CRs with energies $\leq 3\times 10^{15}{\,\rm eV}$ deviating them from their original directions. Therefore, the only possibility to study these cosmic accelerators individually is to consider different kinds of emission resulted from interactions of accelerated particles with magnetic fields, with photons, or with other particles.

SNRs are observed in all electromagnetic domains, from radio waves to TeV gamma-rays. There is a wealth of information produced by many experiments but only a fraction of it has been exploited.

Fluxes and spectra, both integrated and spatially resolved, are widely used. In some cases, fine details of the SNR surface brightness distribution have been deeply investigated. For instance, the radial thickness of the X-ray SNR rim is associated to the strength of the magnetic field (MF) behind the shock front (Vink & Laming, 2003) if the dominant factor of the fading process is the density decrease and the energy downgrading (and/or diffusion) of the emitting particles (Berezhko et al., 2003; Berezhko & Völk, 2004; Bamba et al., 2005; Völk et al., 2005; Warren et al., 2005). Alternatively, it may indicate the length-scale of the MF damping downstream (Pohl et al., 2005, 2015; Ressler et al., 2014; Tran et al., 2015). Additional constraints on models come from the energy dependence of the rim thickness (Ressler et al., 2014). The presence (or absence) of the X-ray emission in the precursor gives hints on the particle acceleration properties (Long et al., 2003; Morlino et al., 2010; Winkler et al., 2014). The azimuthal variation of the radio brightness depends on the obliquity dependence of the injection efficiency (Fulbright & Reynolds, 1990; Petruk et al., 2009a) while the X-ray azimuthal profile might help in determining whether the maximum particle momentum is time- or loss-limited (Petruk et al., 2011a). High spatial resolution observations may reveal the distance between the shock and the contact discontinuity (e.g. Miceli et al., 2009). A considerably reduced distance can be effect of particle acceleration with an efficient back-reaction on the flow (Warren et al., 2005) or otherwise of protrusions of ejecta beyond the forward shock (Rakowski et al., 2011). Similar effects may be due to the ejecta clumping (Orlando et al., 2012). Rapidly varying spots in the SNR shell might be a sign of the presence of either a highly amplified MF and ongoing acceleration (e.g. Uchiyama et al., 2007) or of the plasma waves and steady electron distributions (Bykov et al., 2008, 2009). The ordered stripes observed in Tycho SNR not far from the shell (Eriksen et al., 2011) are a challenge for models of plasma microphysics (Bykov et al., 2011; Malkov et al., 2012; Laming, 2015).

In contrast to the widespread analysis of spectra and local spatial characteristics of the nonthermal emission from SNRs in different spectral ranges, the maps of the surface brightness and of the polarized synchrotron emission – which are known for many SNRs – are lacking an adequate analysis. Such situation should be changed by developing new analysis methods. In general, there are two ways to deal with SNR images: to analyze observed maps with minimal assumptions or to model maps starting from basic theoretical principles.

Using the properties of the emission processes and the observed maps in different bands it is possible, for example, to separate the thermal and nonthermal X-ray images out of the mixed observed one (Miceli et al., 2009), to predict a gamma-ray image of SNR (Petruk et al., 2009c) or to determine the MF strength in the limbs of SNR (Petruk et al., 2012).

As to the latter strategy to explore SNR images, a method to model the synchrotron radio and X-ray images of Sedov SNRs (that is the shell-like adiabatic remnant of a spherical explosion in a uniform medium) was developed by Fulbright & Reynolds (1990) and Reynolds (1998), who also analyzed the main features of the these maps. The method was adapted to the gamma-ray SNR images due to the inverse-Compton emission (Petruk et al., 2009b, 2011b) and to those due to the hadronic interactions with protons inside SNR (Beshley & Petruk, 2012). The method was generalized to the case of a SNR evolving in an ISM with nonuniform distributions of density and/or MF: the resulting asymmetries on the radio maps have been studied by Orlando et al. (2007) and in X-rays and $\gamma$-rays by Orlando et al. (2011).

The SNRs morphologies shown by their maps have the potential for understanding the plasma microphysics (Reynolds, 2004) and determining the properties of the SNR and its environment. In particular, they allow for determination of the 3-dimensional orientation of ISMF and even its gradient around SN 1006 (Petruk et al., 2009a; Bocchino et al., 2011) or, by analyzing a sample of SNRs, to reveal properties of the Galactic MF (West et al., 2016). By modeling the surface brightness maps, it was shown that a certain class of SNRs, namely the thermal X-ray composites (or mixed-morphology SNRs, i.e. remnants with a thermal X-ray peak within the radio shell), might be explained as a projection effect of SNR evolving in essentially nonuniform environment (Petruk, 2001).

The work of Bandiera & Petruk (2016) (hereafter Paper I) and the present paper continue a series of articles devoted to development of methods to model SNR images. Namely, we are interested in maps of the surface distribution of the polarized radio emission that could also be important to test theories. It is known that random MF generally reduces the polarization of the synchrotron emitting sources (e.g. Stroman & Pohl, 2009), but it is able also to produce stochastic small-scale patchy structures, like spots and filaments, with high polarization fraction in synchrotron X-rays (Bykov et al., 2008, 2009). Recently, Schneiter et al. (2015) have considered a task to simulate polarization images of SNR as a whole under two simplifications: the authors neglected the turbulent component of MF and assumed the (variable) Faraday effect in the SNR interior to be minor comparing to the (uniform) Faraday rotation of the polarization planes outside. In our model, we consider both these components and demonstrate that they might essentially modify polarization images comparing to cases with vanishing turbulent MF and internal Faraday rotation.

With Paper I, we have essentially reached two main achievements. First, we have developed a method to evaluate the synchrotron polarized emission in the case in which the MF is the composition of an uniform plus a random component, and we have computed exact analytic formulae for all Stokes parameters in some particular cases, the most useful probably being that of a power-law particle energy distribution and a Gaussian isotropic MF random component: the resulting formulae show strong similarities with those for the classical case (i.e. of a homogeneous MF), and are in fact a generalization of them. Then, using a thin-shell approximation, we have simulated the observed maps of the total radio intensity and of the Stokes parameters for a wide number of cases. In spite of the simplification used, we have been able to tackle in detail a number of effects, like for instance the effectiveness of the internal Faraday rotation in distorting the polarization pattern as well as on affecting the observed polarization fraction; and the partial depolarization also due to projection effects and to the random MF component. Another advantage of that approach is that it is computationally rather light.

In the present paper we present instead a more detailed model; namely, we consider the distribution of polarized synchrotron emission not only in the thin shell but everywhere inside a shell-like SNR. In this way more reliable maps can be obtained, which may be used in particular to test various effects as well as to verify how accurate the simplified models of Paper I can be, under typical conditions.

2 Polarization in adiabatic SNR

In the present paper, aimed at modeling the radio emission from shell-type SNRs, we have improved the treatment presented in Paper I, by implementing a detailed and self-consistent treatment of SNRs in the adiabatic regime, valid for objects old enough to have swept up from the ambient medium a mass much larger than the original mass of the supernova ejecta, but that are still evolving adiabatically, i.e. the radiative losses and the back reaction of accelerated particles are inefficient (in fact, the SNR statistics is consistent with their radio emission to be strongly reduced about the end of the adiabatic phase; Bandiera & Petruk, 2010). The approach we present here (Sects. 2.2, 2.3, 2.4, 3) is general, in the sense that it may be used with the numerical 3D MHD simulations of an adiabatic remnant of an asymmetric supernova explosion in a medium with nonuniform distributions of density and/or magnetic field. However, in order to simplify the analysis, to make conclusions clearer and to have a basic reference model, in the present paper, we apply our model to the so-called Sedov SNR, i.e. the remnant of a spherical point-like explosion in the uniform medium (Sect. 2.1). Polarization maps of SNRs from fully numerical 3-D MHD simulations with asymmetric ejecta and/or expansion in a nonuniform environment are in preparation.

2.1 Magnetohydrodynamics and relativistic electrons in Sedov SNRs

Before calculations of the polarization maps, one has to simulate the 3D MHD structure of an SNR.

In the simulations reported in the present paper, we assume a strong unmodified shock (with $\gamma=5/3$), adiabatically expanding in a uniform medium, after a spherical point-like explosion. This problem admits an analytic solution (Sedov, 1959), which is self-similar: namely, the spatial and temporal variation of parameters inside the SNR, for example the flow velocity $u(r,t)$, can be written in the form

[TABLE]

where the index ‘s’ corresponds to the immediately post-shock position and the profile $\bar{u}(\bar{r})$ does not explicitly depend on time. In other words, it is sufficient to calculate once the normalized profile $\bar{u}(\bar{r})$, to obtain (just by scaling) the profile at any given time. The scaling properties of the Sedov solution reduce the number of free parameters in simulations and reveal how the supernova and the environment characteristics combine.111The Sedov solution describes the remnant of a point-like explosion and therefore neglects the structure of ejecta which may give some effect in the polarization maps at the early stages of SNR evolution. Accurate approximations of the Sedov solution have been proposed for calculations, either in the Eulerian coordinate $r$ (Cox & Franko, 1981) or in the Lagrangian coordinate $a$ (Sect. 4.5 in Petruk, 2000).222The Lagrangian approach is useful if one needs to trace the evolution of parameters inside a given fluid parcel. $r$ is the common spatial coordinate, while $a$ is like a number attached to a fluid element. It is defined as the radial coordinate of the element before the shock crosses this element and remains unchanged later; i.e. $a\equiv R(t_{\mathrm{i}})$ where $t_{\mathrm{i}}$ is the time when the fluid element $a$ was shocked. There is a possibility to convert $a$ to $r$ if one knows the hydrodynamical structure downstream. Such conversion is a way to find a spatial distribution of parameters.

The model assumes that the MF is a composition of an ordered ($B$) and a disordered ($\delta B$) component. The distribution of $B$ inside the SNR may be calculated following Chevalier (1974): the MF in each point may be treated as the sum of a radial and a tangential component, which evolve in different ways (for detailed references and expressions see Sect. 2.2 in Petruk et al., 2016). In a Sedov SNR, like for the other hydrodynamic parameters, also the spatial profiles of these two quantities might be expressed with self-similar functions of the normalized coordinates $\bar{r}$ or $\bar{a}=a/R$. The MF experiences a compression at the shock, with the compression factor $\sigma_{\mathrm{B}}={B_{\mathrm{s}}}/{B_{\mathrm{o}}}$ which depends on the shock obliquity angle $\Theta_{\mathrm{o}}$ (the angle between the ambient MF and the shock normal) according to the law:

[TABLE]

where the index ‘o’ refers to the pre-shock position, and $\sigma=\rho_{\mathrm{s}}/\rho_{\mathrm{o}}$ is the (density) shock compression factor. Consistently with our assumptions, we will use $\sigma=4$ in our calculations. Therefore, $\sigma_{\mathrm{B}}$ runs from $1$, for a parallel shock ($\Theta_{\mathrm{o}}=0^{\mathrm{o}}$), to $4$, for a perpendicular shock ($\Theta_{\mathrm{o}}=90^{\mathrm{o}}$). Instead, the model for evolution of the disordered MF component will be presented in Sect. 3.

In order to synthesize the SNR synchrotron maps one has to know also the spectral and spatial distributions of relativistic electrons. We assume the power-law spectrum of emitting particles with the index $s=2$ and the injection efficiency (defined as the fraction of accelerated particles) independent of the shock obliquity (unless otherwise stated). The evolution of the relativistic electron population downstream of the shock is also self-similar; we model it following Reynolds (1998); details are presented in Petruk et al. (2011b, Appendix A).

2.2 Polarized emission

Once the distributions of MF and emitting particles have been evaluated for all locations within the SNR, we can compute the local Stokes parameters, ${\cal I}$, ${\cal Q}$, and ${\cal U}$, which contain full information on the (linearly) polarized synchrotron emission per unit path; while ${\cal V}$, the parameter associated to circular polarization, vanishes. We can then integrate them along the line of sight, to get the observed Stokes parameters.

For a general orientation $x$-$y$ of the projected axes ($z$ is the line-of-sight coordinate with the positive direction toward the observer), the standard formulae give ${\cal Q}=\langle E_{\mathrm{x}}E^{\mathrm{*}}_{\mathrm{x}}\rangle-\langle E_{\mathrm{y}}E^{\mathrm{*}}_{\mathrm{y}}\rangle$ and ${\cal U}=\langle E_{\mathrm{x}}E^{\mathrm{*}}_{\mathrm{y}}\rangle+\langle E_{\mathrm{y}}E^{\mathrm{*}}_{\mathrm{x}}\rangle$. If we choose the special orientation ($x^{\prime}$-$y^{\prime}$), such that the $E^{\prime}$ vector is directed towards the unit vector $\hat{x}^{\prime}$ (namely the $B^{\prime}$ vector is directed towards the unit vector $\hat{y}^{\prime}$), we readily derive ${\cal Q}^{\prime}\geq 0$ and ${\cal U}^{\prime}=0$.

Let $P_{\perp}$ and $P_{\|}$ be, respectively, the synchrotron power per unit frequency $\nu$ polarized perpendicular and parallel to the projected MF direction, then ${\cal I}^{\prime}=(P_{\perp}+P_{\|})/4\pi$ and ${\cal Q}^{\prime}=(P_{\perp}-P_{\|})/4\pi$. For an arbitrary orientation, the Stokes parameters can be derived using the standard rotation laws:

[TABLE]

where $\chi$ is an angle between the $x$ axis and the local polarization measured counter-clockwise, from the observer point of view. This is the same as an angle between the $y$ and the local vector $\mathbf{B}_{\mathrm{xy}}$ (a component of the ordered MF $\mathbf{B}$ in the plane orthogonal to the line of sight333In the present paper, $\chi$ is the rotation angle of the MF vector while in the Paper I it was the rotation angle of the coordinate system. The sign in this definition is opposite to that in Paper I, and now gives an orientation consistent with the angle $\beta$, for the Faraday rotation. It also applies to (7), cf. Eqs. (10) and (69) in Paper I (see Fig. 1). On the other hand, this change of definition does not affect the results.).

In terms of the components $B_{\mathrm{x}}$ and $B_{\mathrm{y}}$ of the ordered projected MF, we get $\cos\chi=B_{\mathrm{y}}/B_{\mathrm{xy}}$, $\sin\chi=-B_{\mathrm{x}}/B_{\mathrm{xy}}$, where $B_{\mathrm{xy}}^{2}=B_{\mathrm{x}}^{2}+B_{\mathrm{y}}^{2}$ (see Appendix A). Therefore,

[TABLE]

Before integrating the Stokes local values along the line of sight, to obtain observed Stokes parameters, we must account also for the rotation of the polarization plane inside the SNR due to Faraday effect (called ‘internal Faraday rotation’ throughout this paper):

[TABLE]

where $\beta$ is the angle of rotation of the polarization plane during propagation from the emission site $l$ inside the SNR to $z_{1}$, namely where the line of sight crosses the edge of the SNR near to the observer (see Fig. 1), measured in the same direction as $\chi$ (i.e. counter-clockwise from the observer point of view); $\lambda$ is the wavelength of emission, and $RM(l,z_{1})$ is the rotation measure, defined as:

[TABLE]

The foreground Faraday rotation (i.e. from the point $z_{1}$ to the observer location) is not considered in our simulations because our goal is to investigate the propagation effects internal to the source (SNR) rather than those in the interstellar medium.

We detect the composition of the emission from each point $l$ along the line of sight, propagated until the point $z_{1}$. Therefore, the projected Stokes parameters should be calculated as

[TABLE]

where $z_{2}$ is the farther edge of the SNR (Fig. 1).

The angle-sum trigonometric identities allow us to use Eqs. (4) and (5) here.

Once we have the projected Stokes parameters calculated for each point of the SNR projection, we may obtain the map of polarization fraction

[TABLE]

and of the angle of the ‘observed’ rotation of the polarization plane

[TABLE]

The angle $\Psi$ gives us the orientation of the observed ‘projected’ MF ${\cal B}$ relative to the projection of the ambient MF $\mathbf{B}_{\mathrm{o}}$:

[TABLE]

In the discussion of results (Sect. 5), we consider also an average polarization fraction for a SNR image $\langle\Pi\rangle$ which is an arithmetic average of $\Pi$ from each (out of $N$) point over the SNR projection:444We would like to make clear that this average polarization fraction is not how observers typically calculate the polarization fraction for an entire object. The standard practice is the total polarized power divided by the total flux density (our Eq. 12 does not weight regions by total flux). Instead, Eq. (12) gives an average of the local fractions over the image and it is used solely to characterize the dominant values of $\Pi$ in the image, not the object as a whole.

[TABLE]

It should also be noted that when gradients in Stokes parameters occur, unless they are fully resolved, they will entail an effective ‘beam depolarization’. That is, the polarized fractions will be resolution-dependent and can be expected to rise with increasing $N$. In our simulations, we use a high resolution, therefore our local polarization fractions may reach the maximum theoretical value. It is important to account for this effect when comparing the simulations with observations.

2.3 Stokes parameters. Completely ordered MF

In the classic synchrotron emission theory, the Stokes parameter ${\cal I}^{\prime}$ for particles with the single momentum $p$ is

[TABLE]

and the parameter ${\cal Q}^{\prime}$ is

[TABLE]

where $B_{\mathrm{xy}}$ is the component of MF in the projection plane (i.e. perpendicular to the line of sight), $F$ and $G$ are known functions, $\nu_{\mathrm{c}}$ the critical frequency.

The Stokes parameters for electrons distributed with the power law in momentum $N(p)=Kp^{-s}$ are given by convolution over $p$:

[TABLE]

where

[TABLE]

and $c_{1}=3e/(2\pi m^{3}c^{5})$. In the case of the ordered MF, the polarization fraction is the maximum possible one:

[TABLE]

Results of simulations of the polarized emission from the Sedov SNR with the only ordered MF are presented in Sect. 5.1.

2.4 Stokes parameters. Ordered plus disordered MF

If a model considers the only ordered MF then it results in the high fraction of the polarized emission. Observations reveal rather small $\Pi$ in SNRs, around $15\%$. Therefore, a model of polarized synchrotron emission should also include the disordered component of MF inside SNR.

At this point we face the two problems.

First, the classic theory of synchrotron emission is developed for MF which is uniform on the lengthscale much greater than the gyroradius of the emitting electron, in other words, it may be used with the ordered MF only. Thus, the theory of synchrotron emission should be generalized to cases when the random MF $\delta B$ is not negligible comparing to the ordered MF $B$. We have developed such generalization in Paper I. Namely, if the particles with the power-law momentum distribution emit in MF which has ordered component of the strength $B$ and the random component represented by the spherical Gaussian555 A possible realization of a turbulent MF by a spherical Gaussian, which preserves the zero divergence of MF, by an ensembles of Alfvén waves, is given by Bykov et al. (2008). with the standard deviation $\delta B$ then the Stokes parameters are (for details see Sect.2.2 in the Paper I)

[TABLE]

where ${{}_{1}F_{1}}(a,b,z)$ is the Kummer confluent hypergeometric function, ${\cal I}^{\prime}_{\mathrm{o}}$ and ${\cal Q}^{\prime}_{\mathrm{o}}$ are respective Stokes parameters for the case of the ordered MF (i.e. given by Eqs. 15 and 16). The polarization fraction is

[TABLE]

The terms in braces $\{...\}\rightarrow 1$ for $\delta B/B\rightarrow 0$ and the results of the classic synchrotron emission theory are recovered for a vanishing random component of MF. Fig. 2 demonstrates the dependence of ${\cal I}^{\prime}$, ${\cal Q}^{\prime}$ and $\Pi$ on the strength of the disordered MF. Note, that ${\cal I}^{\prime}$ increases with $\delta B/B$; this is of importance for fitting of the observed synchrotron spectra. In particular, if one assumes $\delta B/B\sim 1$ than the flux is twice the flux given by the classic synchrotron theory. Similar effect, namely, that the turbulent MF is able to enhance the synchrotron emission, is noticed also by Bykov et al. (2008).

The second problem is that we have to know the value of $\delta B$ everywhere inside SNR. An approach to model the evolution of the disordered component of MF on the shock and in the SNR interior is one of the main goals of the present paper. It is developed in Sect. 3 where the equation for the evolution of waves is considered. Solutions of this equation are analyzed in Sect. 4. In our approach, the wave properties are related to properties of CRs accelerated on the forward shock and evolved downstream.

Once $B$ and $\delta B$ are known in each point of SNR, the projected Stokes parameters may be calculated by Eq. (7) with ${\cal I}^{\prime}$ and ${\cal Q}^{\prime}$ given by Eqs. (19) and (20). Note, that in our approach to the synchrotron emission theory (where the random MF is represented by the spherical Gaussian) neither $\chi$ (spatial structure of MF) nor $\beta$ (internal Faraday rotation) are affected by the disordered MF.

Results of simulations of the polarized emission from the Sedov SNR including disordered MF are presented in Sect. 5.2.

3 Modelling disordered component of MF in SNR

One of the basic ingredients needed to calculate detailed synchrotron polarization maps of SNRs is the spatial distribution of the MF disordered component, everywhere inside the SNR. In this section we present a description of the evolution of this MF component: this treatment is rather general, in the sense that it may be applied to a shock moving in a medium with nonuniform distribution of density and/or MF; however, for the calculation of the numerical models presented in this paper, we will apply it only to cases of a SNR Sedov expansion. The present treatment is based on the following main assumptions: i. the disordered MF is essentially the effect of the composition of Alfvén waves and we are limited therefore to small values of the ratio $\delta B/B$; ii. these waves, as well as the CRs, behave like “test-particles”, i.e. they are developed atop of a MHD structure, which is unaffected by them; iii. MHD instabilities are not taken into consideration once we consider the Sedov SNR; iv. waves and CRs interact with each other through resonant interactions only; v. the wave damping is described by a non-linear Landau approach.

More sophisticated descriptions of the random MF component are present in the literature (e.g. Bell, 2004; Balsara & Kim, 2005; Giacalone & Jokipii, 2007; Mizuno et al., 2011; Guo et al., 2012; Fraschetti, 2013, and others) and should be considered in the future studies. In the present paper, we adopt instead a rather standard and ‘simple’ prescription for the wave evolution. This allows us to reveal the basic trends caused by the disordered MF, as well as the complexity which arises when taking into consideration the random MF. Also the choice of modelling a Sedov SNR in the present paper is in order to provide some basic ‘reference’ results, for a typical evolutionary phase.

We leave to future simulations the task of revealing effects of other, more complex, elements (like MHD instabilities, nonuniform ISM or ISMF, non-adiabatic or CR-modified shocks, other types of waves) in the modelling of polarized synchrotron emission from SNRs.

3.1 Equation for the downstream evolution of $\delta B$

The equation that describes the evolution of waves downstream of the shock is (McKenzie & Völk, 1982)

[TABLE]

where $u$ is the flow velocity in respect to the shock reference frame, $E_{\mathrm{w}}=\delta B^{2}/4\pi$ is the wave (magnetic and kinetic) energy density, $P_{\mathrm{w}}=\delta B^{2}/8\pi$ is the pressure of waves, $\sigma_{\mathrm{w}}P_{\mathrm{w}}$ and $\Gamma_{\mathrm{w}}P_{\mathrm{w}}$ are the growth and damping rates for waves respectively. Note that the quantities in this equation are not functions of the wave number, but are the integrated ones. Let us assume that the energy flux $F_{\mathrm{w}}$ is dominated by Alfvén waves, so that

[TABLE]

It accounts for the purely magnetic (the first term to the right) and kinetic (the second term) contributions (McKenzie & Völk, 1982), $v_{\mathrm{A}}$ is the Alfvén velocity

[TABLE]

Therefore, behind the shock, $E_{\mathrm{w}}=2P_{\mathrm{w}}$, $F_{\mathrm{w}}=3uP_{\mathrm{w}}$ and equation for $P_{\mathrm{w}}$ becomes

[TABLE]

Alfvén waves interact resonantly with Larmor radius of accelerated protons: $k\propto 1/p$ where $k$ is the wavenumber. The growth rate is derived by Skilling (1975b):

[TABLE]

where $W(k)$ is the spectrum of the wave energy density, $\hat{b}$ is the unit vector in the direction of magnetic field (Skilling, 1975a) and $\hat{n}$ is the unit vector outward along the radius. These unit vectors are defined in such way in order their dot product to be $\hat{b}\cdot\hat{n}=\cos\Theta$ where $\Theta$ is a ‘local‘ obliquity angle (between magnetic field and the radial direction). The absolute value of the dot product is taken because the growth rate is positive while the dot product may be negative depending on the local orientations of these vectors. The growth term in Eq. (22) is

[TABLE]

We therefore integrate (26) using the the property $dk/k=-dp/p$ and the definition of the CR pressure

[TABLE]

and derive that

[TABLE]

The damping term consists of the linear and non-linear contributions (e.g. Ptuskin & Zirakashvili, 2003): $\Gamma_{\mathrm{w}}=\Gamma_{\mathrm{L}}+\Gamma_{\mathrm{NL}}$. The first one is due to collisions of protons with neutral hydrogen of the density $n_{\mathrm{H}}$ (e.g. Kulsrud & Cesarsky, 1971). It is taken $\Gamma_{\mathrm{L}}\propto n_{\mathrm{H}}=0$ assuming that no neutral survive after passage across the shock front. There is no commonly accepted approach to the description of the non-linear damping (see e.g. Ptuskin & Zirakashvili, 2003, and references therein). We take the expression for the rate of non-linear Landau damping from Ptuskin & Zirakashvili (2003), namely from their Eq. (12):

[TABLE]

where $c_{\mathrm{K}}=(2C_{\mathrm{K}})^{-3/2}$, $C_{\mathrm{K}}\approx 3.6$ being the quantity introduced by Ptuskin & Zirakashvili (2003), and $E_{\mathrm{w}}(>k)$ the energy density of waves with wave numbers larger than $k$. Note, that there is no explicit obliquity dependence in this term because the waves interact with themselves and this process does not depend on the shock obliquity. Eq. (30) can be integrated over $k$, giving:

[TABLE]

Multiplying this expression by unity

[TABLE]

and noting that the average $\left\langle kE_{\mathrm{w}}(>k)\right\rangle$ is dominated by the lowest wavenumber in the spectrum ($k_{\mathrm{min}}$): $\left\langle kE_{\mathrm{w}}(>k)\right\rangle\simeq k_{\mathrm{min}}E_{\mathrm{w}}$ where $E_{\mathrm{w}}(>k_{\mathrm{min}})=E_{\mathrm{w}}$, we have the damping term in Eq. (22) of the form:

[TABLE]

Thus, Eq. (22) is now

[TABLE]

In order to follow more easily the evolution of waves in a given fluid element, we use the Lagrangian coordinate $a$ and the relations

[TABLE]

where ‘E’ and ‘L’ denote Eulerian and Lagrangian derivatives respectively, and the latter relation directly follows from mass conservation ($\rho r^{2}dr=\rho_{\mathrm{o}}a^{2}da$.) In terms of $a$, Eq. (34) has the form

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

Eq. (36) is the Riccati differential equation and may be solved numerically. Solutions of Eq. (36) are considered in detail in Sect. 4. In the rest of this section, we describe how to model different properties of waves which are needed to solve this equation.

The solution $P_{\mathrm{w}}(a,t)$, with known $r(a,t)$, $\rho(a,t)$, $u(a,t)$, $B(a,t)$, $P_{\mathrm{c}}(a,t)$ describes the wave pressure. The evolution of $\delta B$ downstream is given by

[TABLE]

This dependence on $a$ may be converted to the spatial dependence on $r$ with the use of the relation $r(a,t)$ which is known explicitly from the hydrodynamic solution.

3.2 Shock compression factor for waves

Let us define the ‘shock compression factor’ for the wave pressure as the ratio of pressures immediately downstream $P_{\mathrm{ws}}$ and upstream $P_{\mathrm{wo}}$:

[TABLE]

This ratio may be calculated from the stationary form ($\partial P_{\mathrm{w}}/\partial t=0$) of Eq. (25) written as

[TABLE]

By integrating this equation across the shock, we get:

[TABLE]

where $\tau_{\mathrm{\sigma}}=\sigma_{\mathrm{w}}^{-1}$ and $\tau_{\mathrm{\Gamma}}=\Gamma_{\mathrm{w}}^{-1}$ are the time-scales for the wave growth and damping. In order to come to this equation, we estimated the functions $\sigma_{\mathrm{w}}$ and $\Gamma_{\mathrm{w}}$ with their average: $(\sigma_{\mathrm{ws}}+\sigma_{\mathrm{wo}})/2$ and $(\Gamma_{\mathrm{ws}}+\Gamma_{\mathrm{wo}})/2$.

If the time-scale $\Delta t$ for the shock to cross the fluid element is much smaller than all $\tau_{\mathrm{\sigma s}}$, $\tau_{\mathrm{\sigma o}}$, $\tau_{\mathrm{\Gamma s}}$ and $\tau_{\mathrm{\Gamma o}}$ then the right side in (43) vanishes and then

[TABLE]

where we have used the relation $u_{\mathrm{o}}=\sigma u_{\mathrm{s}}$. An important result is that, since $\sigma$ is constant for strong shocks, the factor $\sigma_{\mathrm{Pw}}$ neither depends on time nor on the shock obliquity. The shock jump in $\delta B$ is then simply

[TABLE]

3.3 Model for the wave pressure on the shock $P_{\mathrm{ws}}$

The energy density of waves is assumed to reach the saturation level before the shock front. Its value $P_{\mathrm{wo}}$ is given by the stationary form (Amato & Blasi, 2006) of Eq. (25) written in the upstream region:

[TABLE]

Let us introduce two parameters

[TABLE]

[TABLE]

Note from Eqs. (29) and (33) that $\xi$ and $\lambda$ are independent of $P_{\mathrm{w}}$. Eq. (46) becomes, after multiplication by $\lambda\xi/(uP_{\mathrm{co}})$,

[TABLE]

The solution of (49) may be derived by separation of variables, assuming $\xi=\xi_{\mathrm{o}}$ and $\lambda=\lambda_{\mathrm{o}}$ constant on the scale of the precursor, i.e. from the shock to distance $x_{\mathrm{o}}$ upstream of the shock. Thus, immediately upstream

[TABLE]

where $\Omega_{\mathrm{o}}=x_{\mathrm{o}}/\lambda_{\mathrm{o}}$. The CR and wave pressures are set to zero beyond the distance $x_{\mathrm{o}}$ because our model does not consider seed CRs and waves (note that $0\leq P_{\mathrm{wo}}/(P_{\mathrm{co}}\xi_{\mathrm{o}})\leq 1$).

In order to see the meaning of the parameters $\xi_{\mathrm{o}}$ and $\Omega_{\mathrm{o}}$, we use Eqs. (29) and (33), written for the upstream region and use some simplifications. For instance, the minimum wavevector $\mathbf{k}_{\mathrm{min}}$ (whose modulus appears in Eq. 33) is set to $\mathbf{k}_{\mathrm{o}}\simeq\hat{b}_{\mathrm{o}}r_{\mathrm{Lo}}^{-1}$ where $r_{\mathrm{Lo}}=p_{\mathrm{max}}c/(eB_{\mathrm{o}})$ is the Larmor radius of protons with maximum momentum in the field $B_{\mathrm{o}}$. The value of the CR pressure gradient in the shock precursor may be estimated as $\partial P_{\mathrm{c}}/\partial r\simeq P_{\mathrm{co}}/x_{\mathrm{o}}$ where $x_{\mathrm{o}}$ is the distance the same protons are able to diffuse off the shock (the size of the CR precursor); it is directed opposite to $\hat{n}_{\mathrm{o}}$. The unit vectors $\hat{b}_{\mathrm{o}}$ and $\hat{n}_{\mathrm{o}}$ are in the directions of the upstream MF and of the normal to the shock respectively; their dot product is $\hat{b}_{\mathrm{o}}\cdot\hat{n}_{\mathrm{o}}=\cos\Theta_{\mathrm{o}}$. The distance $x_{\mathrm{o}}\simeq D(p_{\mathrm{max}})/u_{\mathrm{o}}$ (Amato & Blasi, 2005) where $D(p_{\mathrm{max}})=r_{\mathrm{Lo}}c/3$ is the Bohm diffusion coefficient. So,

[TABLE]

[TABLE]

where $2c_{\mathrm{K}}\simeq 0.1$, and $M_{\mathrm{Ao}}={u_{\mathrm{o}}}/{v_{\mathrm{Ao}}}$ is the Alfvén Mach number of the shock.

Now, the post-shock pressure of waves is given by

[TABLE]

3.3.1 Limit cases of the wave evolution

Let us introduce the variable $\alpha$ as the ratio of the wave growth and damping rates

[TABLE]

Following the definition (47), immediately before the shock

[TABLE]

Comparing this relation with the solution (50), we see that the ratio between the two rates is related to the ratio of the exponents:

[TABLE]

The plot $\alpha_{\mathrm{o}}(\Omega_{\mathrm{o}})$ is shown on Fig. 3. The two limiting regimes of the wave evolution are evident from this plot: the limit A of the small damping rate at the shock ($\alpha_{\mathrm{o}}=(2/\Omega_{\mathrm{o}})^{2}$ is larger than unity, $\Omega_{\mathrm{o}}\ll 2$) and the limit B of the same rates at the shock ($\alpha_{\mathrm{o}}$ is unity, $\Omega_{\mathrm{o}}\gg 2$). The maximum damping rate may not be larger than the growth one, as expected (we do not assume a seed turbulence).

Fig. 3 demonstrates that the transition between the two limiting cases is in a rather narrow range. Namely, the first one may be used already for $\Omega<1$ and the second one for $\Omega>4$.

It is apparent from the figure that the ratio of the growth to damping rates increases with obliquity, for the same $\xi_{\mathrm{o}}$ and $M_{\mathrm{Ao}}$, because $\Omega_{\mathrm{o}}$ decreases with $\Theta_{\mathrm{o}}$, Eq. (52).

Eq. (53) gives the possibility to calculate the wave pressure at the present time, on the shock The wave pressure on the shock evolves as it follows from Eqs. (53) and (56) and the expression for this evolution may be found as a ratio:

[TABLE]

where the index ‘i’ refers to the time $t_{\mathrm{i}}$, i.e. $P_{\mathrm{wsi}}$ and $P_{\mathrm{csi}}$ are the the wave and CR pressures in the element $a$ when it was shocked. This formula is general; its limits (which are useful for interpretation of results) are considered in the next two subsections.

3.3.2 Limit case A: small wave-damping at the shock

If the damping rate in the shock precursor is negligible compared to the growth rate (limit A), then Eq. (46) becomes

[TABLE]

with the obvious solution

[TABLE]

The same expression may be obtained from the more general formula (53). Namely, if $\Gamma_{\mathrm{w}}P_{\mathrm{w}}\ll\sigma_{\mathrm{w}}P_{\mathrm{w}}$ then $\Omega_{\mathrm{o}}\rightarrow 0$. Decomposition of the exponential terms with small arguments in (53) into series and use of the first-order term and Eq. (52) result in (59).666This way to derive Eq. (59) proves also the correctness of Eq. (52).

The CR pressure is smooth across the shock: $P_{\mathrm{co}}=P_{\mathrm{cs}}$. The post-shock pressure in form of waves is therefore

[TABLE]

in this limit.

In a non-uniform medium with a density profile $\rho_{\mathrm{o}}(R)$, the velocity $V$ of the strong adiabatic shock is well approximated by $V(R)\propto\left(\rho_{\mathrm{o}}(R)R^{3}\right)^{-1/2}$ (Hnatyk & Petruk, 1999, this formula is exact in the case of a constant density). Numerically, $u_{\mathrm{o}}=V$. In the general case of an ambient medium with non-uniform distributions of density and MF, the variation of the post-shock value of the wave pressure in time follows from (60). Namely, the evolution of the wave pressure on the shock may be found as the ratio

[TABLE]

The term in braces $\{...\}=1$ for a uniform ISMF.

3.3.3 Limit case B: equal rates at the shock

If $\Gamma_{\mathrm{w}}P_{\mathrm{w}}=\sigma_{\mathrm{w}}P_{\mathrm{w}}$ at the shock (limit B) then $\Omega_{\mathrm{o}}$ is large. Unity may be neglected comparing to the exponential term both in the nominator and denominator in (53), i.e. $\alpha_{\mathrm{o}}=1$ and the expression for pressure of waves comes simply from (55):

[TABLE]

Therefore, for a general case of non-uniform ISM and ISMF,

[TABLE]

The term in braces $\{...\}=1$ for ISM with uniform distributions of density and MF. Equation (63) is derived as a ratio $P_{\mathrm{wsi}}/P_{\mathrm{ws}}$ directly from (62) using expression (51) and the approximation $V(R)\propto\left(\rho_{\mathrm{o}}(R)R^{3}\right)^{-1/2}$.

3.3.4 Obliquity dependence of $P_{\mathrm{ws}}$

In order to simulate polarization images, we need also to know how does the wave pressure vary with obliquity. We will use Eq. (53) for this purpose, and have therefore to find the obliquity dependences for all the quantities there.

The compression factor $\sigma_{\mathrm{Pw}}$ is constant, Eq.(44).

The obliquity dependence of $\xi_{\mathrm{o}}$ comes from Eq. (51). In order to find it, we write the ratio of this parameter for a parallel shock (marked with the additional index $\|$) and for the shock with a generic obliquity angle (without additional index). Then we have, for a general case of non-uniform ISM and ISMF, that

[TABLE]

The term in braces is equal to unity for uniform ambient density and MF. The obliquity variation of $\Omega_{\mathrm{o}}$ is given by definition (52) and (64):

[TABLE]

In deriving Eqs. (64)-(65), the approximate formula $V(R)\propto\left(\rho_{\mathrm{o}}(R)R^{3}\right)^{-1/2}$ was used.

Now, we know how $\xi_{\mathrm{o}}$ and $\Omega_{\mathrm{o}}$ vary with obliquity. Thus, Eq. (53) gives us the ratio

[TABLE]

which represents the obliquity dependence of the wave pressure. This reduces, for the Sedov shocks, to

[TABLE]

in the limit A and to

[TABLE]

in the limit B.

3.4 Evolution of the CR pressure $P_{\mathrm{c}}$

For a given momentum distribution $f(p)$ of CRs, the CR pressure is

[TABLE]

In order to simplify the calculations of the evolution of $\delta B$, we approximate $P_{\mathrm{c}}$ assuming an isotropic test-particle CR distribution $N(p)=4\pi p^{2}f(p)$ which is represented by a power law $N(p)=Kp^{-s}$, with $s$ close to $2$, in the range between the minimum and maximum momenta of CRs $p_{\mathrm{min}}$ and $p_{\mathrm{max}}$. The CR pressure in this case is

[TABLE]

Since an approximate treatment of $P_{\mathrm{c}}(r)$ is sufficient to our purposes, in this subsection, we fix $s=2$, a choice that minimizes the dependence of the CR pressure on $p_{\mathrm{max}}$, while describing reasonably well also the cases with $s$ slightly different from $2$.

Then, in the regime $p_{\mathrm{min}}\ll mc\ll p_{\mathrm{max}}$, being $v\simeq c$ for $p\geq mc$, we finally derive

[TABLE]

In all reasonable cases the logarithmic term can be approximated with $\ln\left({p_{\mathrm{max,s}}}/{mc}\right)$, where $p_{\mathrm{max,s}}\geq 10^{3}mc$.

The normalization $K_{\mathrm{s}}$, related to the injection efficiency $\eta$ by

[TABLE]

evaluates

[TABLE]

At the shock ($a=R$), the expression for the CR pressure, Eq. (71) with (73), coincides with Eq. (30) in Amato & Blasi (2009).

The minimum momentum may be estimated by its relation to the thermal momentum. We assume that the distribution of the thermal particles downstream may be approximated by the Maxwellian

[TABLE]

where $y=p/p_{\mathrm{th}}$ with $p_{\mathrm{th}}=\sqrt{2mkT_{\mathrm{s}}}$ is the thermal momentum and

[TABLE]

is the shock temperature, $\mu_{\mathrm{o}}$ the mean particle mass in units of the proton mass. Further, we assume that all thermal particles with momenta $p>p_{\mathrm{min}}$ are to be accelerated. Thus, the injection efficiency is

[TABLE]

In general, in order to find the minimum momentum $p_{\mathrm{min,s}}=y_{\mathrm{min}}p_{\mathrm{th}}$, one has to solve this equation to find $y_{\mathrm{min}}$ for a given value of $\eta$. This formula gives, in particular, $y_{\mathrm{min}}=3.07$ for $\eta=3\times 10^{-4}$ and $y_{\mathrm{min}}=3.91$ for $\eta=10^{-6}$. Since $y_{\mathrm{min}}$ varies quite slowly for a wide range of $\eta$, we set in our simulations $p_{\mathrm{min,s}}=3.4p_{\mathrm{th}}(T_{\mathrm{s}})$.

Taking all these pieces together,

[TABLE]

where all parameters are referred to the time $t$ and some of them depend on the shock obliquity.

Let us parametrize the time dependence of the injection efficiency as $\eta(t)\propto V(t)^{-(b+1)}$ with a parameter $b$ (i.e. $K_{\mathrm{s}}\propto V^{-b}$ because $p_{\mathrm{min}}\propto V$). Then the ratio of CR pressures needed in (57) is

[TABLE]

The way to calculate $K(a)/K_{\mathrm{s}}$ is described in Orlando et al. (2007, 2011) for general case of nonuniform ISM and ISMF and in Petruk et al. (2011b) for a Sedov shocks in ISM with uniform density and MF. In particular, the evolution of the normalization $K$ behind the Sedov shock is self-similar: $K(a,t)/K_{\mathrm{s}}(t)=\bar{K}(\bar{a})$ where

[TABLE]

is time-independent, $\bar{a}=a/R$, $R$ the radius of the shock. The profile of the CR pressure, Eq. (77), is (to the level of approximation used for the logarithmic factor in Eq. 71) self-similar as well, i.e. it may be written in the form $P_{\mathrm{c}}(a,t)=P_{\mathrm{cs}}(t)\bar{P}_{\mathrm{c}}(\bar{a})$ where

[TABLE]

[TABLE]

We have just described the time evolution of the CR pressure, at the parallel shock. We also need the prescription for the whole SNR surface, i.e. at the oblique shocks as well. The obliquity dependence of $P_{\mathrm{cs}}$ is due to the respective dependence of the injection efficiency and variations (if any) of the ambient density and the shock radius. It is given by the ratio

[TABLE]

derived from (80) where we neglected for simplicity the obliquity dependence of the logarithmic term and used the approximate formula $V(R)\propto\left(\rho_{\mathrm{o}}(R)R^{3}\right)^{-1/2}$.

3.5 Evolution of $k_{\mathrm{min}}$

The function $q_{2}$, Eq. (39), contains $k_{\mathrm{min}}$. In our model, the downstream evolution of $k_{\mathrm{min}}$ is described as follows. Since $W(k)$ is a decreasing function, the damping is highest for waves with the smallest $k$, Eq. (30). Therefore, we assume that, in a given fluid element, at a given time, the only survived waves are those with the wavenumbers larger than $k_{\mathrm{min}}$ which is able to interact resonantly with the Larmor radius of particles of the highest energy: $k_{\mathrm{min}}\simeq 1/r_{\mathrm{L}}(p_{\mathrm{max}})=eB/(cp_{\mathrm{max}})$. Then, we consider the evolution of $B$ and $p_{\mathrm{max}}$ within this fluid element.

In order to describe $p_{\mathrm{max}}(a,t)$, we parametrize the time evolution of the maximum momentum which particles reach at the shock as $p_{\mathrm{max,s}}(t)\propto V(t)^{q}$; then

[TABLE]

Numerically, for the Sedov phase, $q=1/3$ if $p_{\mathrm{max}}$ is time-limited777In this case, $p_{\mathrm{max}}\propto V^{2}t$ (Lagage & Cesarsky, 1983) and $V\propto t^{-3/5}$ for the Sedov shock; therefore, $p_{\mathrm{max}}\propto V^{1/3}$., $q=0$ if it is escape-limited and $q=1$ if it is limited by the radiative losses (Reynolds, 1998, p. 379). Downstream, even protons with the maximum momentum $p_{\mathrm{max}}$ experience only adiabatic losses, i.e. $p_{\mathrm{max}}(a)\propto E_{\mathrm{max}}(a)\propto\rho(a)^{1/3}$ (e.g. Reynolds, 1998); the energy losses due to pion production is neglected assuming that the shock is propagating in medium with density less than $\sim 10^{4}{\,\rm cm^{-3}}$ (e.g. Beshley & Petruk, 2012). The evolution of $p_{\mathrm{max}}$ downstream is therefore given by

[TABLE]

Taking (83) and (84) together,

[TABLE]

It appears from the above formulae, that the evolution of the minimum wavenumber is self-similar for the Sedov shock:

[TABLE]

with

[TABLE]

[TABLE]

4 Solutions of equation for waves

The evolution of the wave pressure $P_{\mathrm{w}}$ within a fluid element $a$ after its passage through the shock is given by Eq. (36). In this section, we consider solutions of this equation. Sect. 4.1 deals with a particular case which allow for analytic expression. The formulae are presented for the general case and for a particular model of Sedov shocks in uniform medium and uniform MF. Numerical solutions are presented in Sect. 4.2.

4.1 Limit case A: small wave damping

4.1.1 Solution for a general case

If the damping rate is negligible at the shock (i.e. $\Omega_{\mathrm{o}}<1$) and everywhere downstream then the damping term in Eq. (36) may be set $q_{2}P_{\mathrm{w}}^{2}=0$ and the solution of Eq. (36) is

[TABLE]

[TABLE]

In order to evaluate this solution numerically, one needs to know the time history of the whole SNR structure; in other words, one should have the MHD datacube for each time moment $t_{\mathrm{i}}$. It is more suitable for calculations to express the integrals in terms of $a$ rather than $t^{\prime}$. With the use of

[TABLE]

one has the spatial profile of the wave pressure

[TABLE]

[TABLE]

$P_{\mathrm{wsi}}\equiv P_{\mathrm{ws}}(t_{\mathrm{i}})$ is the value at time $t_{\mathrm{i}}$, while other quantities are referred to time $t$. In this approach, it is enough to know the present-time SNR structure only. Having the present time value of $P_{\mathrm{ws}}$ and the ratio $P_{\mathrm{wsi}}/P_{\mathrm{ws}}$ from Eq. (61) (which, like the integrals in (92) and (93), depends on the distribution of MHD parameters at the present time), we may calculate $P_{\mathrm{w}}(a,t)$ numerically.

4.1.2 Solution for a Sedov shock

Let us consider this solution for a particular model, namely, for a Sedov shock in a uniform medium and uniform MF.

For a Sedov SNR evolution in a medium with uniform density and MF Eq. (92) becomes

[TABLE]

where the profiles $P_{\mathrm{wa}}(\bar{a})$, $P_{\mathrm{wb}}(\bar{a})$ are time independent, $\bar{a}=a/R$. Namely,

[TABLE]

[TABLE]

[TABLE]

In order to derive Eq. (95), we made use of Eqs. (61), (78) and $u_{\mathrm{s}}=V/\sigma$. In derivation of (96), we made use of Eqs. (60), (44) and $V(R)/V(a)=\bar{a}^{3/2}$.

Eqs. (95)-(97) are written in a way to express the solution through the basic profiles of MHD parameters. It is clear from here that the downstream evolution of $P_{\mathrm{w}}$ is self-similar. Eq. (40) shows that $\delta B(a)$ is self-similar as well, $\delta B(a,t)=\delta B_{\mathrm{s}}(t)\cdot\overline{\delta B}(\bar{a})$, with

[TABLE]

and $\overline{\delta B}=\bar{P}_{\mathrm{w}}^{1/2}$, i.e.

[TABLE]

The profiles $\overline{\delta B}(\bar{r})$ are shown on Fig. 4 in comparison with profiles of $\bar{B}(r)$. In the considered limit A, the turbulent MF has a flat distribution in the deep interior (no wave damping) while $B$ decreases. As a results, the ratio $\delta B/B$ increases toward the center of SNR. The obliquity affects the normalized profiles $\overline{\delta B}(\bar{r})$: the normalized strength of the turbulent MF is larger for a perpendicular shock. However, in the considered limit, the postshock value $\delta B_{\mathrm{s}}\propto P_{\mathrm{ws}}^{1/2}\propto\cos\Theta_{\mathrm{o}}^{1/2}$ (Eq. 67), i.e. it zero for perpendicular shock.

4.2 General solution for Sedov shock

In a general case when both the wave growth and damping should be considered, the solution of Eq. (36) may be found numerically.

Let us consider Sedov SNR evolution in medium with uniform density and MF and rewrite Eq. (36) in order to have an equation for $\bar{P}_{\mathrm{w}}(\bar{a})$. First, we substitute Eq. (36) with the self-similar form

[TABLE]

Then we use the same procedure as before, namely, transform from the variable $t$ to the variable $a$ with Eq. (91). These steps yield the following equation for the normalized wave pressure $\bar{P}_{\mathrm{w}}(\bar{a})$:

[TABLE]

with

[TABLE]

The factor ${\Omega_{\mathrm{o}}\alpha_{\mathrm{o}}^{-1/2}}{2^{-1}}=(\Omega_{\mathrm{o}}/{2})^{2}$ in the limit of small wave damping. Therefore, $q_{2}\rightarrow 0$ in this limit because $\Omega_{\mathrm{o}}\ll 2$. In addition, ${\Omega_{\mathrm{o}}\alpha_{\mathrm{o}}^{1/2}}{2^{-1}}=1$ and we recover results of Sect. 4.1.2.

In the limit of equal rates $\alpha_{\mathrm{o}}^{1/2}=1$. It should also be $\bar{q}_{2}\bar{P}_{\mathrm{w}}^{2}=\bar{q}_{0}$ at the shock. This yields

[TABLE]

where we have used the relation $\cos\Theta_{\mathrm{o}}=\sigma_{\mathrm{B}}\cos\Theta_{\mathrm{s}}$. Though the expression (105) comes from a limit, we shall use it for any $\alpha_{\mathrm{o}}$. Therefore, we adopt

[TABLE]

With (81), (79) and, for Sedov solutions (Hnatyk & Petruk, 1996, 1999)

[TABLE]

we have

[TABLE]

We use $s=2$ and $\gamma=5/3$ for plots and numerical simulations.

In order to have the spatial profile $\bar{P}_{\mathrm{w}}(\bar{a})$, we should start from the initial value of $P_{\mathrm{w}}(a,t_{\mathrm{i}})$ in the fluid element $a$ at time $t_{\mathrm{i}}$ when it was shocked, i.e. from $P_{\mathrm{wsi}}(a)$; its normalized version is $\bar{P}_{\mathrm{wsi}}(\bar{a}_{\mathrm{ini}})$. Then the Eq. (101) has to be used to evolve the wave pressure in this fluid element up to the time $t$, i.e. to find $P_{\mathrm{w}}(a,t)$; its normalized version is $\bar{P}_{\mathrm{w}}(1)$. Note that $\bar{P}_{\mathrm{w}}(1)$ is different for different initial values $\bar{P}_{\mathrm{wsi}}(\bar{a})$. In other words, the solution $\bar{P}_{\mathrm{w}}(\bar{a})$, $\bar{a}\in[\bar{a}_{\mathrm{ini}},1]$, of Eq. (101) with the initial value $\bar{P}_{\mathrm{wsi}}(\bar{a}_{\mathrm{ini}})$ does not represents the spatial distribution of the wave pressure for different $\bar{a}$ but the time evolution of $P_{\mathrm{w}}$ within the element $a$ written in terms of time in an indirect way. Therefore, in order to have the spatial variation $\bar{P}_{\mathrm{w}}(\bar{a})$, one need to solve this equation many times with the initial conditions for different $\bar{a}_{\mathrm{ini}}$. Different values of $\bar{P}_{\mathrm{w}}(1)$ derived from each solution correspond to the wave pressure in respective fluid elements $\bar{a}_{\mathrm{ini}}=a/R(t)$ at time $t$ and therefore represent such a spatial variation, $\bar{P}_{\mathrm{w}}(\bar{a}_{\mathrm{ini}})$.

This could look like a complicate procedure. However, the direct solution of the equation (36) is not simpler from the point of view of numerical realization. Namely, in our approach, we do not need to keep the time history of spatial distributions of all the parameters to make temporal integration. We just take one, the present-time spatial distribution of parameters and do the spatial integration instead of the temporal one (cf. Sect. 4.1).

The initial condition, as derived from Eqs. (57), (51) and (78), is:

[TABLE]

where we substituted $\Omega_{\mathrm{oi}}=\Omega_{\mathrm{o}}\bar{a}^{3(b+3)/4}$ (as from Eqs. 52, 51, 78) and assumed a uniform ambient medium density and MF. In the limit of small damping, this results in

[TABLE]

(cf. Eq. 61); while in the limit of equal rates we obtain (cf. Eq. 63)

[TABLE]

Fig. 5 shows the radial distribution of the wave pressure downstream of the Sedov shock for some values of $\Omega_{\mathrm{o}}$. The larger the $\Omega_{\mathrm{o}}$ the closer the damping rate to the rate of the wave growth at the shock (Fig. 3). However, for larger $\Omega_{\mathrm{o}}$, the values of $\bar{P}_{\mathrm{wsi}}$ were higher at previous times (Fig. 6). Combination of these features results in the larger wave pressure $P_{\mathrm{w}}$ in the interior of SNR for larger $\Omega_{\mathrm{o}}$ (Fig. 5).

How to know which regime of the wave behavior (limit A or B) is more reliable for a certain parameter set? One can note from Eqs. (52), (51) and (80) that

[TABLE]

This means that a typical SNR shock ($\eta^{1/2}c/V\sim 0.1\div 10$) could be either in the small damping limit or in the limit of the equal damping and growth rates. The quasi-perpendicular shocks are always in the regime of the small damping because of the dependence of $\Omega_{\mathrm{o}}$ on the cosine of the obliquity angle.

The solution for the limit of the small damping, Eq. (94), is shown on Fig. 5 by the dot-dashed line. The radial profiles for small $\Omega_{\mathrm{o}}$ values are close to that approximate solution. The differences are evident deep downstream. The solution in the limit A is obtained under assumption that the condition $q_{2}P_{\mathrm{w}}^{2}\ll q_{0}$ is valid everywhere inside the SNR. However, in the solution of the more general equation (101), the ratio $q_{2}P_{\mathrm{w}}^{2}/q_{0}$ increases downstream with distance from the shock. Therefore, the condition $q_{2}P_{\mathrm{w}}/q_{0}\ll 1$ could not be valid at radii lower than some distance from the center. Nevertheless, the differences are not large and the approximate solution may be used if one needs to reduce the computational cost.

At a first glance, it seems that a similar approximate solution might be obtained for the limiting case B. In fact, this case assumes equal rates at the shock. If the growth and damping would be almost the same everywhere downstream then the equation would be much simpler: $d\bar{P}_{\mathrm{w}}/d\bar{a}+\bar{q}_{1}\bar{P}_{\mathrm{w}}=0$. However, the ratio of the damping to growth rates ($\alpha^{-1}$) increases downstream quite quickly with decrease of $\bar{a}$ (less waves survived) preventing us from such a simplification.

Calculations show that profiles $\bar{P}_{\mathrm{w}}(\bar{r})$ are almost insensitive to the value of $q$. Smaller values of $b$ results in higher $\bar{P}_{\mathrm{w}}$ in the SNR interior because the saturated level $P_{\mathrm{wsi}}\propto\bar{a}^{3b/2}$ ($\bar{a}\leq 1$) was larger at previous times for smaller values of the parameter $b$. In order to reduce the parameter space, $b$ and $q$ are set to zero in our calculations.

5 Polarization maps for Sedov SNRs

In this section we present the results of numerical simulations for SNRs evolving in the Sedov regime (adiabatic expansion in the medium with the uniform distribution of density and magnetic field). For these simulations we have used the assumptions and formulae as described in the above sections. In all our images, the projected ambient MF is along the horizontal $x$ axis. Thanks to the semi-analytic nature of our simulations (Sedov solutions are analytic) we are not confined by the fixed grid of the numerical simulations with limited resolution in three dimensions. We can calculate all the necessary parameters in any point inside SNR and we may have as large resolution as we need. In practice, all our images have $100\times 100$ pixels in the projection plane and very large resolution (determined automatically by the integration algorithm) along the line of sight.

5.1 Completely ordered magnetic field

For the sake of comparison with what will follow, let us present first the results of numerical calculations for the case of a fully ordered MF.

Fig. 7 shows maps of the Stokes parameter $I$ (surface brightness), $Q$, $U$, the polarization fraction $\Pi$, the magnetic polarization angle (angle of ‘projected’ MF with the $x$ axis) $\Psi$ and the associated vector map of the ‘MF orientations’ inferred from the projected Stokes parameters, in the absence of the internal Faraday rotation ($\beta=0$), for the aspect angle (between the ambient MF direction and the line of sight) $\phi_{\mathrm{o}}=90^{\mathrm{o}}$. The most prominent features are the quite symmetric patterns for $Q$ and $U$, the smoothly ordered ‘projected’ MF and the high level of polarization: it is close to the maximum theoretical value, $\simeq 0.69$, over the whole projection.

Fig. 10 of Paper I shows the same set of images but simulated in the thin-layer approximation. The maps from the present three-dimensional model (Fig. 7) resemble those images rather well. Some minor differences can be noticed, especially closer to the SNR edge: they are essentially consequence of the fact that we now model in details the internal structure of SNR which was neglected in the thin-shell approximation adopted in Paper I. In particular, the most apparent differences are the four ’islands‘ on the map of the polarization fraction, which were not present in the same image in Paper I. They are real, not just a numerical effect, and the reason is the following: in the case of the thin layer approximation, all the material is in an infinitely-thin layer around the SNR shock and there are only two emitting ’points‘ along each line of sight with different orientations of the MF vectors (in the case of an aspect angle of $90^{\mathrm{o}}$ the orientations of the MF vectors are aligned on the two points). In contrast, the present model deals with different MF orientations in each point in the SNR interior. The depolarization on the $\Pi$ image is a bit higher where MF vectors vary to larger extent along the line of sight. In other words, it is an effect of the MF jiggling along the line of sight, and $\Pi$ thus is a ’measure’ of the MF disorder along this line.

Fig. 8 shows the same maps, calculated including the internal Faraday rotation. The polarization fraction is still high on Fig. 8 but the lowest value is now $\Pi=0.522$ that is smaller than on Fig. 7. However, the polarization fraction is still high over the most of the SNR projection, as it is demonstrated by the average value of $\Pi$, $\langle\Pi\rangle=0.686$.

Also in this case, the patterns shown in Paper I (see Fig. 14 there) are reasonably reproduced. The most apparent differences are again in the $\Pi$ map. The regions with the lowest polarization fraction are not lying on the projected boundary, but slightly inside. In addition, this map is not symmetric any more with respect to the horizontal axis: also in this case, because (differently from the thin layer case) the level of the MF jiggling along the line of sight is now no longer symmetric. The internal Faraday effect causes the asymmetry between the upper and lower half of the $Q$ and $U$ maps. This translates in an asymmetry in the map of the polarization fraction. In fact, this effect is not present in the corresponding figures in Paper I (Fig.  14 there) which were simulated in the limit of the infinitely thin layer.

The pattern of this asymmetry depends on the orientation of the ambient MF. In our simulations, it is directed towards the right. Therefore, in the nearer half-sphere of SNR (i.e. for $z\geq 1$, see Fig. 1), the line-of-sight component of the ordered MF ($B_{\mathrm{z}}$) is preferentially positive (toward the observer) in the left half and preferentially negative in the right half of image. The orientation of $B_{\mathrm{z}}$ in the nearer half-sphere of SNR is dominant for the internal Faraday rotation (the emission from the rearer part also passes through this SNR part). Therefore, the rotation is preferentially counterclockwise to the left (from the observer point of view) and clockwise to the right (Fig. 9). The polarization planes are rotated by the same angle in the points symmetric with respect to the horizontal axis. However, the plane-of-the-sky components of MF have different orientation in the upper and lower half of image. As a result, the upper-lower asymmetry appears.

The degree of the asymmetry depends on the ’strength’ of the internal Faraday effect. Eqs. (5) and (6) show that, due to the scaling property of the Sedov solutions, the angle of internal rotation $\beta$ in each point of the projection is proportional to the same product

[TABLE]

where, $B_{\mathrm{o\mu}}$ is the strength of the ambient MF in ${\,\rm\mu G}$, $R_{\mathrm{pc}}$ is the radius of SNR in pc, $\lambda_{\mathrm{m}}$ is the wavelength in m. The parameter is, for example, $\beta_{\mathrm{n}}=2.7$ for $n_{\mathrm{o}}=1{\,\rm cm^{-3}}$, $B_{\mathrm{o}}=3{\,\rm\mu G}$, $R=10{\,\rm pc}$, $\lambda=30{\,\rm cm}$. Thus, the parameter $\beta_{\mathrm{n}}$ determines how strong is the Faraday effect in the SNR interior.

5.2 Models with the partially disordered magnetic field

In this subsection, we present maps of polarized emission in models which consider both the ordered $B$ and disordered $\delta B$ MF components. The parameters for a number of considered models are summarized in Table 1. The basic set of parameters is that for the model labelled as A. Other models have one, or at most two, parameters changed with respect to model A and only these parameters are displayed in the table.

Fig. 10 shows the same set of images as on Fig. 7 and Fig. 8 but now including the Faraday rotation and the disordered magnetic field. We can see that the addition of the disordered MF (with $(\delta B/B)_{\mathrm{s\|}}$ less than unity) does not alter considerably the images of the polarized emission.

However, there are some differences visible. In particular, the surface brightness is somehow higher in the model which includes the partially turbulent MF: the area of the region with brightness between 2.5 and 3 units is larger on Fig. 10a comparing to Fig. 8a. This is because the Stokes parameter $I$ is higher if consider the turbulent MF, Eq. (19); this effect should be more prominent if $\delta B/B$ is larger than unity (Fig. 2). Another feature is visible on the same maps, in the shape of the grayscale border marked by the value 3: on Fig. 10a it is like $\}$ while it is like $)$ on Fig. 8a. The reason of this difference is that the ratio $\delta B/B$ is smaller in the perpendicular shock compared to the parallel one: the shock compression factor for $B$ increases from 1 to 4 respectively while $\delta B$ depends on the cosine of the obliquity angle, Eq. 67; thus $\delta B\rightarrow 0$ toward the perpendicular shock (which is projected as a vertical diameter) and the brightness decreases to the value it has without the turbulent MF.

More prominent differences are in the map of the polarization fraction (Fig. 10d, to be compared with Fig. 8d). As expected, the disordered MF depolarizes emission: the area of highly polarized regions is considerably smaller; these regions are located where the perpendicular shock is projected with $\delta B\rightarrow 0$. The values of $\Pi_{\mathrm{min}}=0.488$ and $\langle\Pi\rangle=0.602$ shows how strong is this effect (to be compared to the maximum possible value $0.692$). The polarization fraction is smaller for the same model if the angle $\phi_{\mathrm{o}}$ between ISMF and the line of sight is smaller. In particular, $\Pi_{\mathrm{min}}=0.007$ and $\langle\Pi\rangle=0.428$ for $\phi_{\mathrm{o}}=60^{\mathrm{o}}$.

As we already noticed in Sect. 5.1, the internal Faraday effect changes the measured polarization directions and the overall effect of FR on SNR images depends on a parameter $\beta_{\mathrm{n}}$, Eq. (113). The role of $\beta_{\mathrm{n}}$ in determining the pattern of the ‘observed MF directions’ is apparent from Fig. 11 where these patterns are calculated for models B-E in Table 1 (see also model A on Fig. 10e,f).888The image of $I$ has the same pattern (though with different amplitude) for the models A-E because $\beta_{\mathrm{n}}$ does not affect the contrasts in the map of this Stokes parameter. The polarization images are sensitive to $\beta_{\mathrm{n}}\geq 0.5$. The larger the parameter $\beta_{\mathrm{n}}$ the stronger the effect of the internal Faraday rotation.

Dickel & Milne (1976) have suggested that a change of the ‘inferred MF’ directions from the radial to the mostly tangential one could be the evolutional trend from young to old SNRs. Since $R$ is increasing with the SNR age, we checked if such temporal trend is visible with increase of $\beta_{\mathrm{n}}$ in our model. In fact, there is no clear trend from preferentially radial to tangential polarization directions in our model (Fig. 11). Since we consider here the Sedov phase only, our results agree with the expectation of Dickel & Milne (1976) that the change in magnetic-field orientation is ‘possibly corresponding to the first and third evolutionary phases of an SNR discussed by Woltjer (1970).’ In fact, i) we do not consider ejecta (Sedov model is just the point explosion) and thus the MHD instabilities which might result in the radial ‘fingers’ in young SNRs (e.g. Orlando et al., 2012) are not accounted for; ii) the radiative losses of the shocked plasma (efficient after the Sedov stage) increase the tangential component of MF in old SNRs (Petruk et al., 2016). Future simulations have to be extended to the SNR ages before and after the adiabatic phase in order to reconsider the idea of Dickel & Milne (1976) which seems to be confirmed by recent observations (Dubner & Giacani, 2015).

Another observed property – decrease of the polarization at lower frequencies – is restored in our simulations. In fact, the polarization is smaller for larger $\beta_{\mathrm{n}}$ (Fig. 11, look at values of $\Pi_{\mathrm{min}}$ and $\langle\Pi\rangle$) and $\beta_{\mathrm{n}}$ inversely depends on the frequency $\beta_{\mathrm{n}}\propto\lambda^{2}\propto\nu^{-2}$.

Till now, we have considered the SNR projections for the aspect angle $\phi_{\mathrm{o}}=90^{\mathrm{o}}$, i.e. when the uniform ISMF lies in the projection plane. The polarization is smaller for other orientations of the ambient MF versus the observer (Fig. 12). We see again that the main factor determining the polarization patterns is the internal Faraday rotation. Our simulations demonstrate that the SNR images do not generally reflect the actual orientation of MF inside SNR because they gives the polarization directions of emission affected by the internal Faraday rotation. The larger the role of the internal Faraday rotation (i.e. the higher $\beta_{\mathrm{n}}$), the less similar is the observed polarization pattern to the actual orientation of MF in SNR. Therefore, observed patterns may not be used to extract the MF orientations in SNR until it is proved that the internal Faraday rotation is negligible (i.e. that $\beta_{\mathrm{n}}<0.5$; for example, observing at a high enough frequency, the Faraday rotation becomes asymptotically negligible, see e.g. the moments approach in Paper I). This is clearly seen by comparing lower panels on Fig. 12 and Fig. 13; in the latter one, the images are calculated for the value $\beta_{\mathrm{n}}=0$ and the maps trace therefore the orientation of MF.

It is visible on our images that the Faraday effect could be responsible for the spatial domains in the maps with abrupt change of the polarization directions. This is a feature observed in many SNRs (e.g. Dubner & Giacani, 2015).

If injection preferentially happens at the parallel shocks (model F) then the images of $I$ and of the polarized emission $\sqrt{Q^{2}+U^{2}}$ have maxima located around the parallel shock regions, where emitting particles reside (Fig. 14). Some increase of the polarization fraction from the SNR as a whole is prominent comparing to the isotropic injection case (compare values of $\Pi_{\mathrm{min}}$, $\langle\Pi\rangle$ and $\Pi_{\mathrm{max}}$ given in the figure captions).

How do the properties of the turbulent MF affect images of polarized emission of SNRs? Two parameters in our model regulate behavior of waves in SNR. The first one, $\Omega_{\mathrm{o\|}}$, is $\Omega$ immediately before the parallel shock. This parameter, defined by Eq. (52), determines the ratio between the rates of the wave growth and damping (see Eq. 56). The value $\Omega_{\mathrm{o\|}}=0.4$ in the basic model A brings us to the limiting case of the small wave damping (Fig. 3). In this regime the pressure of waves decreases to zero toward the regions of the perpendicular shock, Eq. (67). The role of the fluctuating field is therefore prominent mostly around the regions with the parallel shocks (i.e. around ‘poles’ which are on the left and right edges of images with $\phi_{\mathrm{o}}=90^{\mathrm{o}}$). Another limit of the wave evolution, when the growth and damping of waves have the same rates, is represented by the model G. Now, the ratio of the rates does not depend on the obliquity. Thus, the effect should lead to overall decrease of the polarization fraction (because of increase of the volume with turbulent MF) and morphological differences, mostly where the regions with the perpendicular shocks are projected. This is what we see on Fig. 15 (models A and G).

The second parameter which determines properties of the waves is the ratio $\left(\delta B/B\right)_{\mathrm{s\|}}$. We have seen that the differences between images with and without fluctuating component of MF are rather small. The reason is that the role of the disordered MF on that SNR images is not strong because of the small value $\left(\delta B/B\right)_{\mathrm{s\|}}=0.3$. Such value was taken because our model assumes that the nature of the fluctuating MF is Alfvén waves. Theory for these waves implicitly assumes $\delta B\ll B$. Nevertheless, in order to see what could be the effect of the stronger waves, we run simulations also for models H and I, with $\left(\delta B/B\right)_{\mathrm{s\|}}=1$ (Fig. 15). These simulations are done just to imagine an effect of efficient disordered MF on the polarized emission from SNRs. The main effects are the increase of the surface brightness with increase of $\delta B/B$ and the decrease of the polarization fraction, in agreement with Fig. 2.

The model in the present paper is based on a principal condition of unmodified shocks. To this extent, our model is self-consistent: we consider test-particle acceleration, Alfvén waves, no MF amplification in the precursor due to CR streaming. This is in agreement with our assumption of adiabatic SNRs which is typically applied to evolved SNRs where the CR back-reaction is not effective, as commonly agreed, and as we have shown in (Bandiera & Petruk, 2010) with our statistical analysis. We hope that in the future the effects of the nonlinear particle acceleration as well as other modes (including Bell’s one) of waves, on the polarized images of SNRs may be considered, in a self-consistent way. We would like to note that there are almost no studies (see review in the Introduction) where the overall maps of the SNR surface brightness are considered with the non-linear effects included (though there are papers dealing with such effects locally, e.g. close to the shock only). In order to study these effects in the polarization images, one has to invest much more effort than to simulate the brightness maps.999Though some aspects of the problem are studied, to our knowledge, there is neither model nor numerical simulations which consider both MF amplification and its evolution deep downstream in the whole SNR volume, for any obliquity. The reason is the different time and length scales: the amplification involves micro-physics which develops on the scale of the shock precursor while MF evolves downstream on the scales of the SNR radius. In the present paper, we consider the Alfvén waves because their theory is well developed. However, there is an approximate approach to consider some aspects of how the CR shock modification affects the polarization images. It is suggested by Ellison et al. (2004) and Ferrand et al. (2010) and is implemented, in particular, in our papers on SNR images (Petruk et al., 2011b; Orlando et al., 2011). It consists in the use of an ’effective’ adiabatic index $\gamma=4/3$ or even $\gamma=1.1$ which is intended to mimic respectively a moderate and a strong shock modification. With these indices the shock compression is higher: $\sigma=7$ and $\sigma=21$ respectively (to be compared to the classical value $\sigma=4$ for $\gamma=5/3$). The consequences of such an approach are the radially-thinner features in the brightness images (e.g. Fig. 7 in Orlando et al. (2011) and approximate analytical formulation in Petruk et al. (2011b)). In addition, the increased compression ratio decreases the electron spectral index $s$ ($s=1.5$ and $s=1.15$ respectively). This results in a small change of the Stokes parameters $Q$ and $U$ (Fig.1 of the Paper I) and $\Pi_{\mathrm{max}}$ (according to Eq. 18; it is $0.652$ and $0.617$ respectively while $\Pi=0.692$ for $\gamma=5/3$), and is completely ineffective to change the ratio $\Pi/\Pi_{\mathrm{max}}$ (Fig.2 of the Paper I). Therefore, we expect that the non-linear acceleration will create higher-order corrections to our test-particle results which are able to catch the most prominent effects.

6 Conclusions

The present paper is devoted to development of the model and to simulations of images of the polarized synchrotron emission from Sedov SNRs. The main goal is a general analysis of the polarization maps. Namely, we are interested to see how the properties of SNR, its internal structure and the turbulent component of MF affect the polarization patterns of these objects. We leave discussion on certain objects to future studies.

Some aspects of the problem were considered by Schneiter et al. (2015), in limits of the negligible internal Faraday rotation and of the zero turbulent component of MF. The approach we developed in Paper I and in the present paper accounts for both of these components.

In order to model polarized emission from relativistic electrons in SNRs, we have solved the two essential problems. First, the classical theory of synchrotron emission has to be generalized to describe the emissivity in the ordered plus disordered MF. A way to such generalization is presented in Paper I; it restores, in the limit $\delta B\ll B$, the known description. Second, the distribution of the disordered MF component has to be known inside the SNR interior. To this purpose, we developed a model for the evolution of the random component of MF downstream of the adiabatic shock under the assumption that it is due to the growth and damping of Alfvén waves. This type of waves is suitable because the whole task is complex in itself but the physics of the Alfvén waves is more or less known. In doing so, we consider the interactions of waves with particles accelerated by the SNR forward shock. To this extend, our treatment is non-linear: evolution of waves are connected to relativistic particles. However, we do not consider possible influence of particles and waves on the MHD structure of SNR. The evolution of the random MF component is considered to be on the magneto-hydrodynamic ‘background’ of Sedov shock.

Having traced the 3-D distribution of all the relevant physical quantities inside SNR, we can consider the detailed spatial distribution of the polarized radio emission and synthesize the projected images of polarization parameters under different SNR orientations with respect to the observer.

Our assumption of adiabatic SNRs implicitly shifts the attention towards not very young SNRs, namely to cases in which the effects of the non-linear particle acceleration should not be important. Probably it could be worth specifying this: we miss some part of generality, but we are right in the frames of the clearly described model which may be applied to many SNRs.

We have demonstrated that the random MF component reduces the polarization fraction, as expected. A new effect is that the presence of the random MF component increases the synchrotron emissivity (Fig. 2). This produces important consequences for fitting the broadband SNR spectra. If, for example, $\delta B\sim B$ then the synchrotron flux is twice the flux for the case where the random MF is neglected.

The role of the Faraday effect in the SNR interior (i.e. the different rotation of the polarization plane from each point along the line of sight inside SNR) is important in formation of the Stokes parameter maps. The uniform derotation of the observed polarization maps might give a ‘true’ image of MF only if the internal Faraday rotation is negligible, that is the exception rather than the rule. The uniform derotation corrects for the Faraday effect in ISM (between the SNR and the observer) but not for the internal rotation which vary essentially over the SNR projection. The pattern of the polarization maps is affected by the internal Faraday rotation if the parameter $\beta_{\mathrm{n}}$ (Eq. 113) is greater than $\simeq 0.5$.

Our simulations restore the observed properties: i) there are polarization domains with abrupt change of polarization direction in SNR maps; ii) the polarization fraction is small; iii) the polarization fraction decreases for larger $\lambda$.

As discussed in Sect. 5.1, the simulated maps presented here show a general agreement with those already presented in Paper I. The main differences can be noticed close to the SNR projected edge, as a consequence of the fact that we have released here the thin-shell approximation. This leads essentially to two kinds of effects: i) the presence of a shell with a finite geometrical thickness leads to a shift of some features away from the edge, giving rise to features like, for instance, some ‘islands’ now appearing in maps of the polarization fraction; ii) the fact that the direction of projected MF now can change along the line of sight, enhances the level of depolarization, and this effect is even more apparent in the cases with a low local random MF component.

In Paper I we have drawn the basis to treat the synchrotron emission in a partially ordered magnetic field, as well as to model the basic and most important geometrical aspects in a shell-type SNR. In the present paper we have added a treatment of the physics of particle acceleration and evolution, in a test-particle regime. The similarity between the results in these two papers should be seen as a proof of the validity of the simplified approach chosen for Paper I. The simpler thin-layer approach may be used in cases when one need an idea about the overall pattern of polarization images. When, however, the goal is to understand the observed patterns or to study some details in images then 3D simulations are needed. We believe that an advantage of our approach is that it could be implemented, with a moderate effort, to some real cases, with more complicate geometry.

There are observational evidences (Dickel & Milne, 1976; Dubner & Giacani, 2015) that patterns of polarization maps of the young SNRs demonstrate predominantly radial orientation of the projected MF directions while maps of the older SNRs have mostly tangential orientation. Our simulations suggest that, in order to restore this effect numerically, the model has to go beyond the ‘point-explosion’ approach and beyond the adiabatic stage, in order to treat the ejecta-ISM interactions in the young SNRs and the influence of the radiative losses of the shocked plasma on the magnetic field.

Acknowledgements

We thank Pasquale Blasi and Elena Amato for interesting discussions and suggestions, as well as the referee for useful recommendations. This work is part of a project started with funds of PRIN INAF 2010. It is partially funded by the PRIN INAF 2014 grant ‘Filling the gap between supernova explosions and their remnants through magnetohydrodynamic modeling and high performance computing’. All simulations were performed on the computational cluster in Institute for Applied Problems in Mechanics and Mathematics, partially in the frame of the project 0117U002815.

Appendix A Components of magnetic field

In this Appendix, we present the formulae to calculate the components of magnetic field.

Let us consider the reference frame with z-axis directed toward the observer. The ambient MF is in the yz plane (without loss of generality). The aspect angle $\phi_{\mathrm{o}}$ is the angle between the ambient MF $\mathbf{B}_{\mathrm{o}}$ and the line of sight. The corresponding spherical coordinates are $(r,\theta,\varphi)$.

In each point, MF $\mathbf{B}$ is a vector sum of the radial $\mathbf{B}_{\mathrm{r}}$ and the tangential $\mathbf{B}_{\mathrm{t}}$ components which may be re-written as $\mathbf{B}=\bar{B}_{\mathrm{r}}\mathbf{B}_{\mathrm{ro}}+\sigma\bar{B}_{\mathrm{t}}\mathbf{B}_{\mathrm{to}}=\bar{B}_{\mathrm{r}}\mathbf{B}_{\mathrm{ro}}+\sigma\bar{B}_{\mathrm{t}}(\mathbf{B}_{\mathrm{o}}-\mathbf{B}_{\mathrm{ro}}),$ where $\bar{B}_{\mathrm{r}}=B_{\mathrm{r}}/B_{\mathrm{rs}}$, $\bar{B}_{\mathrm{t}}=B_{\mathrm{t}}/B_{\mathrm{ts}}$, indices ’s’ and ’o’ mark the immediate post- and pre-shock values, $\sigma$ is the shock compression factor. Therefore,

[TABLE]

where $\mathbf{n}$ is the unit vector in the radial direction. The projection of $\mathbf{B}$ on the Cartesian x axis is given by the dot product

[TABLE]

where $\Theta_{\mathrm{o}}$ is the shock obliquity angle. Expressions for $B_{\mathrm{y}}$ and $B_{\mathrm{z}}$ are the same but the index ’x’ should be changed to ’y’ and ’z’ respectively. The Cartesian components of $\mathbf{n}$ are $(\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta).$ The components of $\mathbf{B}_{\mathrm{o}}$ are $(0,B_{\mathrm{o}}\sin\phi_{\mathrm{o}},B_{\mathrm{o}}\cos\phi_{\mathrm{o}}).$ Taking all these together, we have the MF components:

[TABLE]

As to the obliquity angle, it is given by the dot product $\cos\Theta_{\mathrm{o}}=(\mathbf{B}_{\mathrm{o}}\cdot\mathbf{n})/B_{\mathrm{o}}$ and is therefore

[TABLE]