ASYMMETRIES IN THE EXPANSION AND EMISSION FROM YOUNG SUPERNOVA REMNANTS

A spatially inhomogeneous mass distribution in the ISM is probably the most obvious cause for an asymmetric expansion of SNRs. For example, the presence of atomic or molecular clouds around the SN site might have a strong effect due to the deceleration of the shock wave when impinging on a denser medium. Also, the decreasing density gradient outward from the Galactic plane may play an important role in shaping an evolved SNR. The degree of asymmetry that can be attained depends on the initial energy of the explosion and the steepness of the density gradient. Examples of this scenario have already been reported in the literature, such as the SNRs 3C 400.2 (Schneiter et al. 2006) and W50 (Zavala et al. 2008).

Let us consider the case where the progenitor star moves with a high velocity with respect to its CSM. If this star is in its RG phase, which is likely to be the case of a Type II SN explosion, there is a dense and slow stellar wind blowing toward its CSM and generating a bubble. As a result of the interaction between this moving isotropic stellar wind and the surrounding medium, a bow shock will develop, giving rise to a cometary structure. In this way, an anisotropy in the density distribution around the progenitor is progressively being built, and when the progenitor explodes, the SN shock wave will propagate through this anisotropic background. Strong departures in both morphology and emission would be observed when the SNR shock reaches the border of the bow shock. This scenario has been found to be appropriate to describe the main features of Kepler's SNR (Bandiera 1987; Borkowski et al. 1992; Velázquez et al. 2006).

In this scenario, the progenitor is a WD that belongs to a binary system. If both stars are sufficiently close to each other, mass can be transferred from the companion (a main sequence, sub-giant, or RG star) to the WD (e.g., Marietta et al. 2000). The total amount of mass transferred can be large enough to break the internal equilibrium of the WD star, triggering its explosion as a Type Ia SN. Assuming that all the mass of the progenitor is ejected in the explosion, the gravitational binding would disappear, in which case the ejecta would receive a "kick" that will separate it progressively from the companion. This "kick" has a velocity (with respect to the mutual center of mass) v_k = M₂v_orb/(M₁ + M₂), where M₁ and M₂ are the progenitor and companion masses, respectively, and v_orb is the orbital velocity, given by

$$\begin{equation} v_{\rm orb}=1.16\times 10^8 \sqrt{\frac{M_1+M_2}{a_{10}}}\ \mathrm{cm\ s^{-1}}, \end{equation} \tag{ 1 }$$

where a₁₀ is the orbital radius in units of 10¹⁰ cm and M₁ and M₂ are given in units of M_☉. In Table 1, we show the parameters assumed for a list of possible companion stars and the different values of the "kick" velocity that will be used in our simulations. In all our runs, we assume that the mass of the progenitor M₁ is 1.4 M_☉. The first column indicates the evolutionary stage of the companion star. The second one gives the mass of the companion, the third one lists the orbital radius, the fourth represents the orbital velocity of the reduced mass of the binary system, and the last one is the "kick" velocity.

Another possible consequence of the presence of a companion star is the so-called mass-loading effect that might take place when the expanding shock encounters the companion and strips a fraction of its atmospheric material. We describe the role of this effect on the asymmetry of the SNR in Section 2.5.

If the mass distribution of the progenitor's interior is itself anisotropic, the subsequent explosion will be anisotropic as well. This could be the case, for instance, for fast rotating stars.

This asymmetry is caused by the impact of the SNR shock wave on the companion star. After this collision, the SNR shock sweeps up material (mass loading) from the upper atmosphere of the companion star. The material is embedded in the SN ejecta, generating an asymmetric density distribution, which would translate into a differential expansion of the remnant. The collision between the shock wave of an SN and its companion star has been studied numerically by Marietta et al. (2000). More recent numerical studies have been performed by Kasen et al. (2004) and Pakmor et al (2008). All these studies focus on the first stages after the collision between the SN shock and the companion star (from a few days up to a few months), where physical processes other than the mass loading are also relevant. Marietta et al. (2000) find that if the companion is a main-sequence star it can lose about 15% of its mass, while an RG companion can lose up to 98% of its envelope. Since the amount of mass stripped from the companion can be an appreciable fraction of the total mass of the progenitor's mass, we believe that this can be an adequate scenario for Tycho's SNR.

We are interested in studying the evolution of Type Ia SNe, such as Tycho's SNR. The initial evolution of this kind of remnant is usually described considering that by the time of the SN explosion, the inner 4/7 of the total ejected mass has a constant density ρ_c up to a radius r_c, while the remaining outer 3/7 of the ejected mass follows a power-law density ρ ∼ r⁻⁷ (Colgate & McKee 1969; Chevalier 1982; Jun & Norman 1996).

Under this assumption, Truelove & McKee (1999) obtained approximate shock trajectories for the free expansion and adiabatic phases, which are given by

$$\begin{equation} R^*_b= 1.06\ t^{*4/7} \end{equation} \tag{ 2 }$$

for the free expansion phase and

$$\begin{equation} R^*_b= (1.42\ t^*\ - 0.321)^{2/5} \end{equation} \tag{ 3 }$$

for the adiabatic or Sedov–Taylor phase. In Equations (2) and (3), R*_b = R_b/R_ch is a dimensionless SNR radius and t* = t/t_ch is a dimensionless time, defined by characteristic length and timescales, respectively, given by (Truelove & McKee 1999)

$$\begin{equation} R_{\rm ch}=3.07 \biggl (\frac{M_*}{n_0}\biggr)^{1/3}\ \mathrm{pc} \end{equation} \tag{ 4 }$$

and

$$\begin{equation} t_{\rm ch}=423\ \frac{M^{5/6}_*}{n^{1/3}_0\ E^{1/2}_{51}}\ \mathrm{yr}, \end{equation} \tag{ 5 }$$

where M_⋆ is the ejected mass in units of M_☉, n₀ is the unperturbed ISM particle density in cm^-3, and E₅₁ is the initial energy explosion in units of 10⁵¹ erg.

The transition from the free expansion to the Sedov–Taylor phase occurs when the amount of swept up mass is larger than the ejected mass. Truelove & McKee (1999) estimated this transition to occur at a time

$$\begin{equation} t_{\rm ST}=0.732\ t_{\rm ch}\ \mathrm{yr} \end{equation} \tag{ 6 }$$

and a radius

$$\begin{equation} R_{\rm ST}=0.881\ R_{\rm ch}\ \mathrm{pc}. \end{equation} \tag{ 7 }$$

Depending on whether the remnant is in the free expansion or in the Sedov–Taylor phase, it is possible to obtain an estimate for E₅₁ combining Equation (2) or Equation (3) with Equations (4) and (5) for historical SNRs because their ages t_f are known and their radii R_f can be inferred. To this end, a value for the ISM particle density n₀ needs to be assumed, which will be one of the inputs of our numerical simulations. It is straightforward to show that E₅₁ is given by

$$\begin{equation} E_{51}= 2.9\times 10^4 (n_0\ M_{\star })^{1/2} \ \frac{R^{7/2}_f}{t^2_f} D^2_{3} \end{equation} \tag{ 8 }$$

if the remnant is in the first phase of evolution or by

$$\begin{equation} E_{51}= 8.9\times 10^4 \frac{M^{5/3}_{\star }}{t^2_f\ n^{2/3}_0} \Bigg [\!6.05\times 10^{-2}\ R^{5/2}_f\ \Bigg (\frac{n_0}{M_{\star }}\Bigg)^{5/6}\!\!\!\!\! {+}\, 0.321\!\Bigg]^2 D^2_3 \end{equation} \tag{ 9 }$$

if the SNR is in the Sedov–Taylor phase. In Equations (8) and (9) t_f is the SNR's age given in years and R_f is the radius given in pc. The distance to Tycho's SNR is ill constrained and still in debate, with values ranging from 1.5 to 4.5 kpc (e.g., Chevalier et al. 1980; Albinson et al. 1986; Schwarz et al. 1980, 1995; Smith et al. 1991; Ruiz-Lapuente et al. 2004). As a working hypothesis, we shall adopt an intermediate value of 3 kpc throughout this paper, which is coincident with the distance suggested by Ruiz-Lapuente et al. (2004). In Equations (8) and (9), D₃ = D(kpc)/3 kpc is a correction factor if a different distance to the source is assumed. In Table 2, we use different values of ISM density (Column 1) and total ejected mass (WD and companion's envelope, Column 2) to compute the radius of the shell, R_ST, at the beginning of the Sedov phase (Column 3), the transition time to the Sedov phase (Column 4), and the initial energy of the explosion (Column 5). In all cases, the final age is 440 years and the final radius 3.5 pc, which is estimated assuming a distance of 3 kpc.

Table 2. Estimate of Initial Energy Explosion and Initial Radius for the Sedov–Taylor Phase for Different Input Parameters

n0(cm−3)	M⋆(M☉)	RST(pc)	tST(yr)	$E_{51} (\mathrm 10^{51} \,\rm erg)$
0.5	1.4	3.8	516	1.0
0.5	1.7	4.1	579	1.1
0.5	2.0	4.3	634	1.2
1.0	1.4	3.0	346	1.4
1.0	1.7	3.2	380	1.6
1.0	2.0	3.4	423	1.7

Download table as:  ASCIITypeset image

The effects of deceleration of the forward SNR blast wave due to the acceleration of cosmic rays have not been considered when deriving Equations (8) and (9). This effect can be observationally studied through the relative location between the forward shock and the contact discontinuity. Warren et al. (2005) carried out this analysis for the case of Tycho's SNR and reported evidence in support of this process. Similar results were reported by Miceli et al. (2009) for the case of SN 1006. Recently, Ferrand et al. (2010) have shown that this mechanism can produce a difference of the order of 14 arcsec in the whole expansion of a remnant such as Tycho. This process appears to have a global effect on the SNR dynamics but it cannot explain the observed asymmetries. Besides, it is not clear how much of the initial energy is spent in accelerating particles at the main shock. For these reasons we have not taken into account this effect in our study.

In this section, we compare two basic sources of asymmetry: a kick velocity, v_k, and a gradient in the ISM density of the form n_ISM(x) = n₀ exp(-x/H), where x is a Cartesian coordinate with its origin at the remnant center and H is a density scale length. We can think of v_k as a perturbation to the expansion velocity profile that translates into an energy perturbation E₁, i.e., $E\simeq 0.5\ m\ (v_0\ \hat{r}\ + v_k \hat{x})^2\sim E_0+E_1$, where E₀ = 0.5 m v²₀ and E₁ = m v₀ v_k cos θ (θ is the angle between the expansion velocity and the "kick" direction). The relative importance between these two effects can be evaluated taking into account Equations (2) and (3). The perturbation caused by both effects can be written as (see also Equation (5))

$$\begin{eqnarray} \displaystyle E^{1/2} n^{1/3}_{\rm ISM} &&\approx E^{1/2}_0 n^{1/3}_0 (1 + \epsilon _k)(1-\epsilon _n)\nonumber\\ &&\approx E^{1/2}_0 n^{1/3}_0 (1 + \epsilon _k - \epsilon _n)\,, \end{eqnarray} \tag{ 10 }$$

where _k = v_k cos θ/v₀ and _n = x/3H. Therefore, replacing Equation (10) into Equation (2), we obtain

$$\begin{equation} \displaystyle (R^*_b)^{7/4} \approx 2.6\times 10^{-3} \frac{E^{1/2}_{51}\ n^{1/3}_0\ t}{M^{5/6}_*}\ \left[1 + \left(\epsilon _k-\epsilon _n\right)\right]. \end{equation} \tag{ 11 }$$

Based on Equation (11), it is possible to assess which of the two effects is quantitatively more important: the kick (characterized by _k) or the density gradient (described by _n). As a result of either _k or _n, R_S(t) in any given direction could depart from the symmetric case. Figure 2 shows the relative importance between _n and _k. The solid straight line represents the condition _k = _n, corresponding to comparable effects between the kick and the ISM density gradient. In region II, the asymmetries are dominated by the kick velocity, while region I shows the zone where anisotropies due to a density gradient are more important. Density gradients can arise, for example, in the presence of local molecular clouds or in a turbulent ISM. The gray area in Figure 2 displays the region of kick velocities smaller than 500 km s⁻¹, which is a reasonable assumption for an SNR like Tycho. We find that even for relatively large values of v_k, between 100 km s^-1 and 500 km s^-1, and for reasonably mild density gradients characterized by H between 5 pc and 50 pc, the most important cause of asymmetry is the density gradient.

Zoom In Zoom Out Reset image size

To allow for the mass-loading effect described in Section 2.5, our simplified model assumes that shortly after the explosion, the mass distribution is given by two axisymmetric cone-shaped regions, each one with a homogeneous density distribution, as sketched in Figure 3. We calculated the total energy and mass for such a density distribution, while maintaining the center of mass at rest. Zone I contains only material from the progenitor star, while Zone II also includes the extra mass provided by the outer envelope of the companion, after the main SNR shock wave collides with it. We have the following set of equations for mass, energy, and momentum:

$$\begin{equation} \displaystyle \hbox{\fontsize{8.5}{7.5}\selectfont ${ \left\lbrace \begin{array}{@{}lllll} M & = & M_I + M_{II} & = & \frac{1}{2} M_0 (1 + \cos \theta) + \frac{1}{2} M_0 (1 - \cos \theta) \epsilon ^{-1} \\ \\ E & = & E_I + E_{II} & = & \frac{1}{2} E_0 (1 + \cos \theta) + \frac{1}{2} E_0 (1 - \cos \theta) \epsilon ^{-1}\\ \\ \rho _I v_I & = & \rho _{II} v_{II} & \Rightarrow & v_I = \epsilon v_{II} \end{array}\right.}$}, \end{equation} \tag{ 12 }$$

where we define the total energy as $E_0 = \case{2}{5}\pi r_{\rm in}^3 \rho _I v_I^2$, the mass as $M_0 = \case{4}{3} \pi r_{\rm in}^3 \rho _I$, and the density contrast as $\epsilon = \case{\rho _{II}}{\rho _I}$. The mass excess, $m_{\rm exc} = \frac{1}{2} M_0 (1-\cos \theta) (\epsilon -1)$, is included in M_II where we have assumed ⩾ 1. Additional parameters are r_in, the initial radius considered, and the aperture of the cone θ₀. For = 1, the problem reduces to the symmetric case with E = E₀ and M = M₀. We perform simulations with three different masses: 1.4 M_☉ for the symmetric explosion case, 1.7 M_☉ and 2 M_☉ for the models that include a binary system with a massive companion.

Zoom In Zoom Out Reset image size

The hydrodynamic numerical simulations shown in this work were carried out with the two-dimensional axisymmetric and a full three-dimensional version of the adaptive mesh Yguazú-a code, which is described in detail by Raga et al. (2000).

In this code, the gasdynamic equations are integrated with a second-order accurate (in space and time) implementation of the "flux-vector splitting" algorithm of van Leer (1982). A system of rate equations for atomic/ionic species is also integrated together with the gas dynamic equations. With these chemical equations, it is possible to compute a non-equilibrium cooling function (with a parameterized cooling rate applied for the high-temperature regime). Rate equations for H i ii, He i iii, C ii iv, and O i iv are considered. The reaction and cooling rates which have been included are described in more detail in Raga et al. (2002) and Raga et al. (2007).

For the two-dimensional axisymmetric simulations, a five-level binary adaptive grid with a maximum resolution of 3 × 10¹⁶ cm was used, in a 10 pc × 5 pc (axial × radial) computational domain. The three-dimensional simulations are also done in a five-level binary adaptive grid but with a maximum resolution of 6 × 10¹⁶ cm. A cubic computational domain was considered with a physical size of 10 pc along each Cartesian direction.

As mentioned above, our main goal is to model an SNR with the characteristics of Tycho's SNR. The average radius of this remnant is ∼3.5 pc if a distance of 3 kpc is considered (Section 3.1), and its age is 440 years. If a number density n₀ of 0.5 cm^-3 is assumed, from Table 2 and Equation (8), we derive a value of E₀ = 10⁵¹ erg for the initial SN energy explosion. This value is in agreement with previous estimates given by Badenes et al. (2006). Employing detailed one-dimensional simulations, Badenes et al. (2006) showed that delayed detonation models with kinetic energies around 10⁵¹ erg can reproduce the fundamental properties of the X-ray emission from the shocked SN ejecta, constraining the explosion of SN 1572 to be a normal Type Ia. In addition they were also able to reproduce some morphological characteristics, such as the positions of the reverse shock and contact discontinuity. This has later been confirmed by the light echo spectrum presented in Krause et al. (2008). Using these values for n₀ and E₀ in Equations (6) and (7), we find that the remnant is close to finish its free expansion phase (see Table 2).

The SNR was initialized by applying radial velocities to a sphere of radius 6.2 × 10¹⁷ cm (21 cells in the finest grid) for the axisymmetric simulations, centered along the symmetry axis at x = 5 pc and r = 0, where a reflective boundary condition is imposed. For the fully three-dimensional simulations, the energy is placed into a sphere of radius 10¹⁸ cm (equivalent to 16 pixels at the maximum resolution), centered on a cubic Cartesian domain.

For the "kick" models (where the ejecta acquires an extra velocity due to the orbital motion of the progenitor when it explodes as a Type Ia SN), we integrate the hydrodynamic equations in the reference frame of the progenitor. This condition is equivalent to assuming that the ISM, with n₀ = 0.5 cm^-3 and a temperature of 10³ K, moves with a velocity v_ISM = −v_k. According to Table 2, we choose "kick" velocities of 0, 200, and 500  km s^-1. Furthermore, we assume that an initial mass for the ejecta (M_⋆) of 1.4 M_☉ is uniformly distributed within the initial remnant radius, i.e., the initial remnant density is constant. The velocity was modeled with a linear profile, with a maximum velocity at the initial SNR radius given by $v_{\rm max}=\ssty\sqrt{\frac{10}{3} E_k/M_{\star }}$. The initial kinetic energy, E_k, was set as 95% of E₀, while the remaining 5% was imposed as thermal energy.⁶

To model an asymmetric explosion due to a mass excess from the companion, we initialize the internal structure of the SN as shown in Figure 3 and Equation (12), for both our two-dimensional and three-dimensional simulations considering mass excesses of 0.3 M_☉ and 0.6 M_☉. For the three-dimensional simulations, in order to break the axisymmetry, we impose a fractal density distribution with a spectral index of −11/3, which is consistent with the turbulent ISM (see Esquivel & Lazarian 2005; Esquivel et al. 2003; Ossenkopf et al. 2006). The fractal density was imposed on the entire computational domain and has an amplitude of fluctuations on the order of 20% of the mean value.

To compare the results arising from our simulations against observations, we also generated synthetic X-ray emission maps.

In the low-density regime, the emission coefficient is calculated as j_ν(n, T) = n²η(T), where n and T are the particle density and temperature distributions obtained from the numerical simulations. The function η(T) depends smoothly on T and is computed over the energy band (0.3–1.6) keV using the chianti⁷ atomic database and its associated IDL software (Dere et al. 1997, 2009), where the X-ray spectrum can be fitted with a thermal plasma model, considering both the line and free–free emission. Furthermore, the effects of the interstellar extinction (considering the extinction curve of Morrison & McCammon 1983) have been taken into account. We used the ionization equilibrium (IEQ) model developed by Mazzotta et al. (1998) (also Landini & Fossi 1991) and the element abundances given by Vancura et al. (1995).