Hydrodynamic simulations of the interaction of supernova shock waves with a clumpy environment: the case of the RX J0852.0-4622 (Vela Jr.) supernova remnant
Observations in all electromagnetic bands show that many supernova remnants (SNRs) have a very aspheric shape. This can be the result of asymmetries in the supernova explosion and of a clumpy circum-stellar medium. We study the generation of inhomogeneities and the mixing of elements due to these two sources in multi-dimensional hydrodynamic simulations of the propagation of a supernova blast wave into a cloudy environment. We model a specific SNR, VelaJr (RX J0852.0-4622). By comparing our results with recent observations, we can constrain the properties of the explosion. We find that a very energetic explosion of several occurring roughly about 800 years ago as well as a supernova with an energy closer to the canonic value of a few 1000 years ago are consistent with the shape and emission of the SNR.
keywords:ISM: supernova remnants
The emission of photons across all wavelengths from many supernova remnants (SNRs) reveals a complex geometry of these objects. Observations indicate global asymmetries (Manchester, 1987; Kesteven & Caswell, 1987; Gaensler, 1998) as well as a small-scale clumpy structure of many remnants (e.g., Park et al., 2004; Fesen et al., 2006; Iyudin et al., 2007). The observed asymmetries of SNRs can result from asymmetries imprinted in the matter already during the supernova explosion, or they can be created when the shock wave propagates into an inhomogeneous interstellar medium (ISM), see, for example, Bisnovatyi-Kogan & Silich (1995) for theoretical considerations on asymmetric appearances of SNRs, and Gaensler (1998); Wang et al. (2002) on the observational interpretation of some asymmetrical SNRs.
The exact nature of the explosion mechanism is not fully understood, but the presence of non-spherical hydrodynamic flows is a common feature of currently discussed explosion mechanisms. Such flows may result, e.g., from convection and the standing accretion shock instability (SASI), both of which play an important role for explosions powered by the heating of matter in the stellar core by neutrinos (see, e.g., Janka et al., 2007), or from global asymmetries caused by rapid rotation of the core, in particular in the presence of strong magnetic fields (Bisnovatyi-Kogan & Moiseenko, 2008). During its non-steady propagation through the onion-shell like structure of the stellar mantle and envelope, the supernova blast wave is subject to Rayleigh-Taylor instabilities leading to clumpy ejecta and an efficient large-scale mixing of the star’s chemical composition (Kifonidis et al., 2003, 2006; Hammer et al., 2010; Wongwathanarat et al., 2010).
The interstellar medium surrounding the explosion site represents a second source of inhomogeneities in SNRs. A non-uniform density of the environment can distort the shock wave even when the explosion itself is perfectly spherically symmetric (Bisnovatyi-Kogan et al., 1990). Various factors are responsible for shaping the spatial structure of the ISM. Most massive stars loose mass through an intense stellar wind that may possess a clumpy structure, or e.g., for rapidly rotating stars, a large-scale anisotropy. Furthermore, massive stars typically are born and evolve in regions of active star formation, where the winds of many main-sequence stars and the ejecta of many previous SNe collide forming dense clouds in the ISM.
In the following we will describe a series of three-dimensional hydrodynamic simulations of the propagation of SN shock waves into an inhomogeneous environment. Apart from studying the effects that modify the global appearance of SNRs in general, our models focus on one particular remnant, Vela Jr. (SNR RX J0852.0-4622). Observations of the environment of this remnant, which have been made in different electromagnetic bands (Dubner et al., 1992; Zinchenko et al., 1995; Dubner et al., 1998; Yamaguchi et al., 1999; Moriguchi et al., 2001; Testori et al., 2006; Pettersson, 2008), reveal the presence of a number of large clouds of atomic hydrogen (H I) close to the current location of the blast wave. They also show an enhanced emission, where the blast wave hits these clouds. While the observations provide valuable insights in the properties of the ISM, they do not allow for an unambiguous determination of the properties of the explosion. Because the mass and energy of the ejecta are not well determined by the observations, it is unknown whether we are dealing with the remnant of an ordinary supernova or a more energetic and possibly also more asymmetric hypernova event (Maeda & Nomoto, 2003).
Age estimates for Vela Jr. range from yrs (Iyudin et al., 2005) up to a few yrs (Katsuda et al., 2009). The question of the age of the SNR is closely linked to the determination of its distance, for which a similar range of values has been obtained. The distance may be as large as 750 pc, as suggested by Katsuda et al. (2008) based on an estimate of the expansion rate of the SNR. Aschenbach et al. (1999), analysing X-ray and -ray data, argued for a much lower distance of around 200 pc. Iyudin et al. (2005) present arguments for a distance at the lower end of the range of values, among these a possible spatial correlation between and the Vela Jr. SNR and other objects in the Vela region (Redman et al., 2002; Pozzo et al., 2000; Woermann et al., 2001). Kim et al. (2012) based on UV measurements of Vela SNR that include also Vela Jr., region E of Kim et al. (2012) give arguments for a smaller distance to this remnant. If these arguments hold, it is possible to use the distance of the Vela SNR as a good estimate for the distance of Vela Jr., placing it at a distance between 300 and 400 pc. We treat the distance and age as unknown parameters and perform simulations for high and low distances.
To test the compatibility of various parametrized explosion models with the observational evidence, we have set up a model for the ISM in the Vela Jr. region placing SN ejecta of a given mass and energy at its centre. We then simulate with a multi-dimensional hydrodynamic code the propagation of the shock wave into the ISM and the expansion of the ejecta.
Our ISM model consists of a spatial distribution of overdense (relative to the ambient medium) clouds similar to that of the clouds inferred from observations of the Vela Jr. region. Concerning the latter, we note that the X-ray image of Vela Jr. might be quite misleading. Although appearing as a proper ellipse in projection (see Figs. 1 in Aschenbach et al., 1999; Iyudin et al., 2007), the X-ray bright parts of the Vela Jr. rim, particular the NW part, are most probably caused by a strong interaction of the forward shock with a H I cloud located in that region (Iyudin et al., 2007, see Fig. 3 in), The properties of the H I cloud can be inferred from those given in Dubner et al. (1998), and its geometry. Taking into account projection effects, this provides an estimate of for the cloud mass, and a range of for the gas density. Together with the ejecta density of (Aschenbach et al., 1999) these estimates constrain the parameters used in our simulations of the Vela Jr. environment.
In our study we varied the explosion parameters within the range compatible with the observations. To this end we used masses and energies for the ejecta given by the anisotropic SN models of Maeda & Nomoto (2003), to test the possibility of a hypernova-like explosion. For comparison, we also simulated a ”standard” supernova explosion with conventional values for the mass and energy of the ejecta.
Only a limited number of simulations of the global structure of SNRs have been performed until present. Earlier simulations assuming spherical symmetry (Itoh & Fabian, 1984; Itoh & Masai, 1989) provided approximate relations between shock radius and SN parameters. However, they did not consider an inhomogeneous ISM, because this requires multi-dimensional models. Three-dimensional simulations similar to ours have been performed by Shimizu et al. (2012), who particularly examined the effect of a stellar wind on the ISM and the break-out of the shock from this wind. Our approach differs from their more general one by considering the specific cloud properties as observed in the Vela Jr. remnant, and by the methods we employ to construct light curves and emissivity maps from the hydrodynamic simulation data. We compute the emission across the entire electromagnetic band assuming a simple cooling function (Tucker & Koren, 1971) in a post-processing step and use tracer particles to determine the mixing of chemical elements in the ejecta.
Important aspects of the evolution of a SNR include the instabilities and turbulence developing in the flow past dense clouds, and the heating and subsequent evaporation of clouds hit by the shock wave. While only few global simulations of the dynamics of SNRs exist, these effects have been studied in great detail using high-resolution simulations of smaller sub-regions. However, these local simulations do not discuss the complex gobal geometry of the SNR. We refer the reader, in particular, to the hydrodynamic models of Cowie et al. (1981); McKee & Cowie (1977); Klein et al. (1994); Nakamura et al. (2006); Pittard et al. (2010); Pittard (2011). The possible influence of magnetic fields, which we neglect in our study, has been considered in boundaries of local models, e.g., by Mac Low et al. (1994); Inoue et al. (2009, 2010, 2012) and in global models, albeit only with two-dimensional models, by Fang & Zhang (2012).
2 Methods and initial data
We used the high-resolution finite-volume code AENUS (Obergaulinger, 2008) to perform a set of three-dimensional hydrodynamic simulations. The code employs high-order reconstruction methods (the simulations presented here were performed with a -order method), monotonicity-preserving (Suresh & Huynh, 1997), and approximate HLL Riemann solvers in the multi-stage (MUSTA) framework (Toro & Titarev, 2006). We assume that the matter is described well by a simple ideal-gas equation of state, and take into account neither ionisation nor heat conduction. The simulations were performed on a cubic Cartesian domain of size with open boundaries in all directions. The coordinate system is oriented in such a way that the -axis points towards Earth. The computational domain is covered by a numerical grid of zones.
We approximate the clumpy environment of the Vela Jr. supernova explosion site by considering two populations of dense clouds, which are embedded into an ambient ISM that has a constant gas number density and temperature . For most models used in our study we assume and .
One population of clouds is based on observations of the Vela Jr. remnant (Dubner et al., 1992; Zinchenko et al., 1995; Dubner et al., 1998; Yamaguchi et al., 1999; Moriguchi et al., 2001; Testori et al., 2006; Pettersson, 2008) and consists of four large clouds situated at a distance of roughly 10 pc from the explosion site. The positions and properties of these clouds are given in Tab. 1. We label the four large clouds by their location with respect to the explosion site Fig. 1: cloud NW, the largest one, consists of three dense cores contained in a larger, but less dense spheroidal structure. The two clouds S and SE are identical spheres, and cloud N, the densest one, is of spherical shape, too. To convert the observed positions of the clouds on the sky and their observed angular diameters into physical lengths, we had to assume a distance to the SNR, which we set equal to 300 pc (see discussion about the distance to Vela Jr. in Iyudin et al., 2007).
The second population of clouds consists of a large number of smaller spherical members, presumably below the resolution limit of the observations, at random positions in the sky, whose number density . The cloud radii are randomly distributed between 0.02 pc and 0.5 pc. This population of clouds was inferred from optical observations of Vela SNR environment, see, for example Pakhomov et al. (2011) and references therein.
All clouds are assumed to be in pressure equilibrium with the surrounding matter, i.e, . The large spherical clouds (Tab. 1) have a density falling off as , where is the distance from the cloud centre. For the spheroidal cloud NW the density is constant on spheroidal surfaces falling off as , where is given by
Here are the spherical coordinates of a point inside the cloud w.r.t. the centre of the cloud, and give the orientation of the (oblate) spheroid, i.e., the direction of its semi-major axis , and is the spheroid’s semi-minor axis.
|S25A||6.0||6.7||reference model, progenitor of 25|
|S40||10.2||10.9||progenitor of 40|
|S25Bs||6.0||1.0||weak spherical explosion|
|S25A4c||6.0||6.7||only the four main clouds|
|S25AC||6.0||6.7||1200 random clouds of|
|S25AD||6.0||6.7||SNR placed at a distance of|
We do not take into account the effects of gravity in our simulations, since gravity does not affect the evolution of the SNR, because the gravitational binding energy of the gas comprises only a small fraction of its total energy. Moreover, the gravitational time scale, i.e. the free-fall time scale of the unshocked ISM, is longer than the time scales we are interested in.
We initiate the explosion by placing gas of a total mass and energy at the centre of the computational grid, and assume that the ejecta have the shape of a prolate spheroid of semi-major axis and semi-minor axis , or of a sphere of radius . All three lengths are in the range of order 1 to 2 pc, i.e., the initial explosion site is covered by 10 to 20 grid zones per spatial dimension, only. The imposed explosions are inspired by models 25A, 25B, and 40A of Maeda & Nomoto (2003). corresponding to progenitors of 25 and 40 solar masses, respectively. The ejecta masses and the explosion energies of our models are given in Tab. 2. As our grid resolution is too coarse to allow for a direct mapping of the Maeda & Nomoto (2003) models, we had to rely on ad-hoc profiles for the hydrodynamic variables of the initial ejecta. For simplicity, we assume that the total explosion energy is composed of 5% thermal energy and 95% kinetic energy, the density and pressure of the ejecta are uniform, and the initial velocity is purely radial and varies linearly with spheroidal distance (see (2)).
In the simulations we would also like to follow the propagation of small clumps of different elements ejected in the explosion through the ISM. However, because the sizes of these clumps are far below the affordable grid resolution, we model these clumps as test particles that are passively advected with the flow, i.e. the particle positions are determined by solving the equation of motion,
where and are the time-dependent position of a particle and the flow velocity at that position, respectively. In the initial models, we place a number of particles at random positions near the explosion site and assign some elemental distribution to each of it. By following their advection, we are able to determine the spatial distribution of elements in the SNR. We distinguish two populations of tracer particles, which we call the outflow component and the spherical component. The outflow component is composed of all particles initially situated in a cone of opening angle around the semi-major axis of the prolate explosion spheroid. The remaining particles make up the spherical component. In our simulations we follow the evolution of the elemental abundances of oxygen, titanium, and nickel.
The collision rates between ions and electrons are very low in the tenuous gas inside the Vela Jr. SNR, typical ionisation ages being around (Aschenbach et al., 1999). Therefore, no equilibrium between ions and electrons can be established within the evolutionary times of interest here, i.e., the two kinds of particles are characterised by two different temperatures, and . When a fluid element is passed by the shock, the ions are heated, but thermal energy is transferred to the electrons only on very long time scales by infrequent collisions. The delayed heating of the electrons has important consequences for the emission of photons by the gas, which is not determined by the temperature of the ions but rather by the lower electron temperature. Therefore, we consider electrons and ions separately in our simulations to avoid overestimating the luminosity of the SNR.
In our treatment of the non-equilibrium ionisation (NEI), we use an expression for the time scale of equilibration between electrons and ions of mass and charge given by Spitzer (1962),
Here, and are the mass and charge of the electron, respectively (charges are measured in units of the elementary charge, ); is the number density of the ions, and is the Coulomb logarithm. Using this time scale, the equation for the internal-energy density of the electron gas, , reads
where the electron pressure is derived from the internal energy by
with , and the
equilibration source term .
Since the most detailed treatment, i.e., evolving equations governing the advection, ionisation and recombination of a large set of elements in all relevant ionisation states, is beyond the scope of this article, we used a simplified approach replacing the set by only one representative ion. We set and , i.e., appropriate for a mixture of hydrogen and helium. As a further simplification, we use a constant Coulomb logarithm, . While such an approach is inferior to a detailed NEI simulation, it allows us to capture the essential NEI effects with an acceptable uncertainty when compared to the uncertainties due to the unknown parameters of the explosion and the ISM.
Our hydrodynamic simulations provide the time-dependent three-dimensional distributions of ejecta parameters, like .e.g, density or temperature. To compare these results with the observed emission of SNR, one has to couple the hydrodynamic evolution of the matter to the radiative transfer of photons emitted, absorbed, and scattered by the gas. Fortunately, the properties of the gas allow for several simplifications to this otherwise very complex task. First and foremost, we note that the gas is optically thin to the radiation we are interested in, i.e., X-ray and UV photons. Moreover, radiation carries away only a small fraction of the total energy of the expanding ejecta. This allows us to determine the emission by post-processing the hydrodynamic simulation results.
Given the hydrodynamic nature of our models, we compute the spectral emissivity of thermal bremsstrahlung according to a simplified version of the expression given by Tucker & Koren (1971), where
Here, is the (electron) temperature in units of . Our implementation assumes only one species of ions; furthermore, we approximate the Gaunt factor by unity. By integrating (5) over wavelength between the wavelengths corresponding to the minimum and maximum photon energies and , we obtain the emission in various energy bands. We will focus our analysis on X-ray and UV bands. We will present light curves of our models in a number of intervals distributed logarithmically over photon energy (Tab. 3). An accurate treatment of other (non-thermal) processes requires a more sophisticate modelling approach. Therefore, we defer such a more complete analysis of the photon emission to a later time.
3 Simulations and results
In the following we will describe the most important results from our simulations. We will discuss the evolution of a reference model in more detail before turning to a series of models in which we vary subsequently different initial parameters.
3.1 Observational facts
We briefly summarise the main observational findings which we later compare with our numerical results:
The geometry of the remnant, as observed in, e.g., X-rays and UV (see Fig. 2 for an example of X-ray observations by Aschenbach et al., 1999) appears to be almost circular. At the level of the brightness chosen in these observations, the axis ratio is close to unity, but varies from waveband to waveband, and with the sensitivity thresholds of the imaging instruments. In addition, several large spots show enhanced emission, particularly near the outer border of the circular X-ray-emitting region.
An estimate of the expansion rate of the remnant of was obtained by Katsuda et al. (2008) for NW X-ray bright rim of the remnant. Because the rim might be bright as a result of the strong interaction of the SNR blast wave with the HI cloud in that region (Iyudin et al., 2007, Fig. 3), this measurement is related to the cloud shock velocity of this bright spot. Katsuda et al. (2008) agree with such an interpretation.
Estimates of the energy emitted in X-rays at four bright spots (at the centre and three along the rim of the SNR) have been obtained by Slane et al. (2001). To convert the fluxes to luminosities, we have to assume a distance to the SNR. Taking the total, absorption corrected flux of the high-temperature X-ray emission component of Vela Jr. as (Aschenbach, 1998) for , we derive the luminosity of remnant in the 0.1-2.4 keV energy range as . The bright spots near the outskirts of the SNR have a luminosity of a few / each (Aschenbach et al., 1999). The emission in extreme the UV is of a similar or larger magnitude (Kim et al., 2012; Nishikida et al., 2006). The observed spectra indicate that the emission of the SNR is dominated by non-thermal processes, most likely synchrotron radiation. Thus, estimates of the luminosity provide an upper bound of the thermal contribution only. Further evidence for non-thermal processes is provided by observations in the TeV range (Aharonian et al., 2007).
These findings show a rather complex structure of the remnant. However, a rather high degree of uncertainty remains, and therefore, they do not unambiguously point to a particular mechanism for creating the inhomogeneous appearance of the SNR. We test one possible, viz. the interaction of the SN blast wave with an inhomogeneous medium containing large clouds with a large overdensity w.r.t. the mean density of the ISM.
A complication arises in the comparison of our simulation results with observational data due to a potential contamination by emission from the nearby (in projection on the celestial sphere) SNR Vela. Furthermore, as noted above, the observed emission is non-thermal to a high degree, i.e., only upper limits can be given for the thermal contribution, and thermal line emission is particularly strong in the soft X-ray band. Our simplified analysis is based on thermal bremsstrahlung only. Thus, a comparison of simulation results and observations is somewhat restricted. Based on his Eqn. (52), Vink (2012) estimates that non-thermal synchrotron emission requires relatively high shock velocities of . Applying this criterion to the hydrodynamic evolution of our models we can assess the importance of non-thermal emission. A more thorough treatment is hindered by lack of detailed information about the magnetisation of the ISM around the Vela Jr. supernova remnant, which can be significantly higher behind the shock than in the ambient medium due to efficient field amplification at the shock surface (Zirakashvili & Ptuskin, 2012). Including a magnetic field into our simulations, while in principle possible, would increase the parameter space considerably, in particular if we take into account that estimates of the field strengths vary widely. Kishishita et al. (2013) estimate for the field in the NW regions of the SNR and for the ISM field, while Berezhko & Völk (2010) find a better agreement with observations if they assume a downstream field of , corresponding to an ambient field of . This uncertainty limits our ability to, e.g., setup the initial condition for self-consistent MHD simulations.
Because of observational uncertainties, but even more because of the unknown initial conditions and the simplifications in our models, we do not expect that one particular model perfectly matches the observations. However, we try to identify simulation results which resemble the observations closest.
3.2 General features
We chose a simulation with an explosion based on the SN model 25A by Maeda & Nomoto (2003) as the first model we discuss. The progenitor of this explosion is a star of 25 solar masses exploding in a prolate geometry, releasing a mass of and an energy of . We approximate the aspherical geometry of the explosion very roughly by a prolate spheroid with an axis ratio of . We follow the evolution of the simulated remnant till an age of .
During the first , the explosion roughly retains its spheroidal shape although the axis ratio decreases somewhat, i.e., the ejecta become slightly more spherical. This evolution, a posteriori, justifies that we assume in our simulations only slightly spheroidal ejecta, because the possibly more pronounced initial asphericity at the onset of the explosion has already decreased during the time it took the shock wave to expand to a size of , when we start our simulations.
The post-shock state consists of tenuous hot gas having a density around and an electron temperature around . The gas is a strong emitter of mainly X-ray photons. The total luminosities in the X-ray and UV bands at are and , respectively. While these values only slightly exceed the luminosities estimated from observations, the emission decreases quickly to a level lower than that derived from observations. We note here that most of the EUV luminosity of the remnant is produced by the nonradiative shocks inside the remnant (Kim et al., 2012; Nishikida et al., 2006; Nichols & Slavin, 2004).
The spheroidal shape of the ejecta is clearly destroyed when the shock wave hits the closest clouds (see density maps in Fig. 3). Though these are small, they are able to divert the ejecta and shield the regions immediately behind them to a certain degree. On impact, the shock wave heats the surface of the clouds, and the heated rather dense gas starts to emit X-ray radiation. This can be seen in the emissivity maps in Fig. 3.
Upon passage of the shock wave, the ionic component of the gas is heated to very high temperatures, but the gas starts to radiate only after the electrons have been heated up to roughly the ion temperature by collisions. Since the frequency of collisions increases with density, equilibration occurs fastest in the densest parts of the SNR. The central regions of the SNR possess a very low density, and though the density is considerably higher right behind the shock wave, it is still low compared to that of the clouds. Consequently, once hit by the shock, the combination of high (electron) temperatures and densities makes the clouds prominent sources of radiation. Thus, the radiation emitted by the SNR depends very strongly on the properties of the clouds.
The radiative energy loss is small, i.e., the gas remains hot and therefore radiates for a long time. Clouds can remain strong emitters of radiation throughout the entire evolution unless they are dispersed by the impact. We find a notable disruption only for clouds of very small radius or low density. Clouds above a radius of remain intact, except for some stripping of surface material and some deformations. The clouds which retain their position, size and shape, act merely as obstacles for the blast wave. This behavior holds even more for the few massive clouds NW, N, S, and SE, respectively. Neither does the shock wave disrupt any of the clouds nor is it able to modify the substructure of the cloud NW. Because of their high densities, the four large clouds are the strongest emitters of radiation in the SNR.
As the shock wave expands, the density and temperature of the post-shock gas decrease, i.e., the luminosity decreases. We show the light curves in several bands in UV and X-ray for model S25A in Fig. 4. At all times, the emission is dominated by photon energies around (green and blue solid lines). For this model, the X-ray emission between 1 and 7 keV (dash-dot-dot-dotted black line) exceeds the emission in the extreme UV band between 80 and 120 eV (dashed black line) by a factor throughout most of the evolution. The X-ray emission exhibits a plateau at a total luminosity before it decreases gradually, while the UV luminosity begins to increase at around , when the blast wave starts to impact the main clouds (the beginning of the interaction is marked by a vertical line in the figure). The X-ray emission peaks at , possibly facilitating detection of the SNR in X-rays. At the shock-cloud interaction is in full swing, and the clouds are being hit by the ejecta (cf. the positions of the tracer particles in the top right panel of Fig. 3). At this point, the flow geometry is already quite complex, leading to mixing in the gas. In both bands, several single bright spots can be identified that coincide with the clouds. In particular, we find a strong brightening of the NW limb of the SNR (see the middle and bottom panels of Fig. 3). These bright spots are the most prominent elements of a network of filamentary emission, similarly to what is observed, e.g., by Dubner et al. (1998). Comparing the emissivity maps at and , we find that the UV and the X-ray emission start to fade after the shock has ceased to heat the surfaces of the clouds; in particular the inter-cloud medium is emitting only very weakly in the UV at . A tentative comparison of the luminosities obtained in the simulations with those estimated on current observations shows that the emission of our models is roughly compatible with the observed values, though it might be on the high end of the observationally allowed range. Note, that we do not consider any influence of the foreground absorption on the luminosities derived by simulations.
The complex flow patterns developing due to the interaction between the shock and the clouds facilitate mixing of elements within the ejecta. This is reflected in both the spatial distribution of the elements given by the tracer particles, and in the velocity spectra of different elements. We show the evolution of the velocity spectrum of the oxygen ejected in the explosion in Fig. 5. The total oxygen mass amounts to , of which are in the outflow component. During early epochs the latter material has a significantly higher velocity () than the oxygen in the spherical component (). During the expansion of the blast wave, and most notably when the explosion interacts with the clouds, the ejecta are gradually slowing down. This affects the outflow component more than the spherical component. By , the velocity difference between the two components has reduced from an initial factor of to a factor slightly above unity (cf. the two distributions in the left panels of Fig. 5). As the particle distributions in the top panels of Fig. 3 show, the spatial distribution of particles belonging to the outflow and the spherical component have similar features. Both populations are characterised by fingers protruding between the heavier clouds. The evolution of the distributions of nickel and titanium, which we place in the initial model only in the outflow component, is very similar to that of the outflow component of oxygen. We note that the braking and broadening of the Ti distribution is more pronounced than that of the Ni distribution. The fluid elements containing Ni have the highest velocities, i.e., they are initially at the front of the outflow. After the generation of fingers in the ejecta material, these fluid elements remain at the largest radii and highest velocities. The mixing is most pronounced for the elements at intermediate radii.
3.3 Variation of the model parameters
3.3.1 Explosion properties
Below we assess the dependence of our results on the properties of the SN explosion, i.e., ejecta mass, explosion energy, and the geometry of the explosion.
Higher explosion energy
We compare the results described above to a simulation using a model based on a progenitor of (model S40). According to the simulations of Maeda & Nomoto (2003), the explosion of this star yields an ejecta mass of and an explosion energy of . For these parameters, the expansion of the shock wave and the mixing of the tracer particles proceeds very much like in the reference model, hitting the large four clouds at the same time. This similarity can be attributed to the fact that the ratio between explosion energy and mass is the same in both models. The light curve of this model (Fig. 6)is qualitatively similar to that of the reference model, but quantitative differences exist. We can identify the same phases of a very rapid decline from the initial level, followed by a secondary rebrightening when the shock wave interacts with the clouds, and a more gradual fading at late times. Throughout all stages, the total luminosity is higher than in the reference model by a factor of order , but the hardness ratio between the X-ray and the UV emission is similar, with the X-ray exceeding the UV by far.
Less energetic explosions
As we have seen, a higher explosion mass does not change the results very much if the explosion energy is increased accordingly to keep their ratio constant. Models S25B and S25BB ejecting a mass of but with an energy of and , respectively, only exhibits a qualitatively similar evolution, although at a slower pace. In model S25B, the blast wave hits the four large clouds at a time of (first interactions with clouds SE and N happen already at ), i.e., about after the reference model with higher explosion energy per unit mass. This would correspond to an explosion date around the century. These events are delayed by another when we decrease the explosion energy in model S25BB. While the dynamics is similar, the light curves of the model (Fig. 8) shows interesting differences to the reference model. The luminosity is, in general, lower and softer, i.e., stronger in the UV than in the X-ray band. The interaction between shock and clouds does not increase the emission. Instead, it gives rise to a long plateau phase of constant emission after , during which the UV and X-ray luminosities are around . Thus, the high total luminosity and in particular the hardness of the emission appears to be a feature distinguishing particularly strong explosions from SNe of a lower energy.
A consequence of the lower explosion energy is the lower expansion velocity of the fluid elements (see Fig. 8), which is reflected in the velocity spectrum of the tracer particles. Furthermore, we find less mixing of particles from the outflow population into the spherical one, mainly because the spherical population of particles remains more concentrated in the central parts of the SNR, with less particles escaping between the clouds close to the centre of the explosion.
We compare the X-ray emission at the outer rim at the time when the shock hits the NW cloud of this model to the reference model in Fig. 9. The overall luminosity is higher for the model with the stronger explosion, but the most striking difference is the lower amount of substructure for lower explosion energies. All other parameters equal, a higher explosion energy leads to more efficient heating of the gas and thus more emission, not only from the surfaces of the clouds, but also from the tenuous gas in the interior regions. Besides the differences in appearance, higher explosion energies also cause a larger shock velocity. To reach the NW cloud at the stages shown in Fig. 9, the high-energy-explosion model requires about 600 years, while the shock wave of the lower-explosion-energy one takes years. We estimate the relative expansion rates, , as the logarithmic time derivative of the radius of the X-ray front along the X-ray-bright ridge around , finding and for high and low explosion energies, respectively. For model S25BB, the value is the same within the measuring uncertainties.
In the bottom panel of Fig. 9, we compare the estimated rate at which the emission front expands into the direction of the cloud NW as a function of its distance from the centre of the explosion. For each time step, we detected the radius, of the emission front and then computed the expansion rate, . During the entire evolution, the emission front slows down for bot models. During the interaction of the shock wave with cloud NW, S25A has a much higher expansion rate than S25B. In both models, the rate at which the shock wave expands into cloud NW is lower by about a third w.r.t. free expansion in a direction without any clouds.
Comparing these values to derived by Katsuda et al. (2008), our simulations appear to favour an explosion energy closer to the canonical than to hypernova-like energies. We point out that this is only a shock velocity in the cloud, but not an estimate of the expansion rate of the SNR shock wave, which can be much higher.
Because taking into account just thermal bremsstrahlung, our analysis is limited to a part of the radiation emitted by the SNR only. Although a full treatment including all radiation processes is beyond the scope of this work, we can estimate the impact of our restriction to thermal bremsstrahlung on our results. To this end we employ a criterion for the most important missing process, non-thermal synchrotron radiation, given by Vink (2012). X-ray synchrotron emission requires a shock velocity exceeding 2000 km/s. This condition is fulfilled for the propagation of the shock in the inter-cloud regions, but is violated for the shock of model S25B hitting cloud NW, where . In contrast, in S25A in that cloud exhibiting considerable fluctuations in time and space. Consequently, our models would emit appreciable amounts of synchrotron radiation in the NW cloud only for a high explosion energy. Hence, the detection of non-thermal emission from that region is consistent with a high-energy supernova, albeit no conclusive evidence for it. Going beyond the qualitative statement that our models should produce X-ray synchrotron emission is difficult. Without the ability to produce detailed synchrotron maps from our numerical simulations, we can only speculate about the spatial distribution of the emission. We might expect a qualitative similarity between the synchrotron emissivity and the thermal emissivity maps presented in Fig. 9, i.e., extended filamentary structures behind the shock. However, we can only resolve this issue with a more elaborate modelling approach. Most properties of synchrotron radiation depend on the magnetic field in the emitting gas. Since our simulations do not include a magnetic field, we can only rely on estimates derived from observations, e.g., by Kishishita et al. (2013). However, the properties of the synchrotron radiation beyond the field-independent cut-off energy are determined by additional parameters which we cannot extract from our simulations. Therefore, we do not attempt a more accurate analysis here.
To explore the effects of the geometry of the explosion, we simulated models S25As and S25Bs, explosions with a mass of and energies of and , but exploding spherically. In the reference model, we deliberately point the explosion spheroid towards cloud NW. The combination of this orientation and the axis ratio of the explosion yields a roughly equal time at which the blast wave collides with all four major clouds. In the spherical case, on the other hand, the different distances of these clouds from the explosion center lead to larger differences in the times of impact: for high explosion energy, the clouds S, SE, and N are hit at , wheras it takes another for the shock to arrive at cloud NW. In the case of low explosion energies, these events take place in the same order, but roughly later.
Of course, the most obvious observation signature of a spherical explosion would be a spherical appearance of the SNR. However, in an inhomogeneous ISM, the shape of the SNR may become distorted during its expansion. Furthermore, as we have seen above, even a bipolar explosion can lead to a moderate axis ratio of the remnant after several hundreds of years. Additionally, parts of an SNR might not be easily accessible observationally, rendering a clear identification of the shape difficult. Thus, additional observational information might be helpful to obtain the properties of the SN that created the remnant. We consider in the following the most distinct properties of the light curves of these explosions.
We display the light curves of models S25Bs and S25As in Fig. 10. Since a large fraction of the total luminosity is emitted by the large clouds, the differences in the dynamics w.r.t. the anisotropic models manifest themselves in modifications of the photon emission. We observe again the common phases of a fast decline early on, a plateau or rebrightening during the interaction with the clouds, and a fading afterwards. Similarly to models S25A and S25B discussed above, the spectra are harder for higher explosion energies, though even for the high-energy explosion the emission is not as strongly dominated by the X-ray emission as for the anisotropic explosions. Note also that the high-energy model S25As emits stronger in UV than in X-rays between and , i.e., when the shock hits the clouds. The low-energy model S25Bs has a very soft emission at late times, since the X-ray luminosity decreases rather rapidly after .
3.3.2 Properties of ISM and clouds
The position and the distribution of sizes and densities of small clouds is poorly constrained from observations. Therefore, we can only hope to model the real state of the ISM in a very approximate manner by inserting clouds of random sizes at random points into the ISM. In two models with and , we did not include any clouds at all (S25Anc) or only the four large ones for which we have observational indications (S25A4c). Consequently, there is no mixing until when the shock wave hits cloud N.
In another series of models (models S25Anc, S25A4c, S25AC, S25Aw, S25Ad0.1, S25Ad1.0, and S25Ad10), rather than reducing the number of clouds, we explored the consequences of the presence of more and denser clouds in the vicinity of the explosion, while keeping the explosion energy and mass the same as in the reference model S25A.
Comparing the results of these models among each other and to those of the model S25A, we can conclude that the expansion velocity of the blast wave is not affected strongly by the absence of a multitude of small clouds. The post-shock flow, on the other hand, shows much less features without the random population of clouds. This is reflected in the distribution of the particles tracing the mixing of elements in velocity space. We find much less dispersion of both the spherical and the outflow components without the additional clouds. At , e.g., the distributions of titanium and nickel in velocity space exhibit a rather pronounced maximum at . A certain broadening of the distributions occurs after the interaction with the four clouds took place. Afterwards, the distribution retains its shape until the end of simulations at . The impact of the ejecta onto clouds clouds leads to both up-scattering and down-scattering of particles in the velocity space. While some parts of the flow are decelerated in interactions with the clouds, other parts are accelerated when the flow is squeezed through the less dense, narrow channels existing between the clouds. On average, though, the energy loss exceeds the energy gain, and the total kinetic energy is decreasing by about 10 % when the random clouds are included.
In model S25AC, we placed 1200 small clouds with a density randomly in the simulation domain. Though this is most likely more than can be expected in the vicinity of the actual SNR, we used this model to explore the general effect of many small clouds. Most prominently, the expansion of the explosion is slowed down by the large number of interactions the shock wave suffers. The luminosity is very high during the entire evolution. Furthermore, we find rather soft spectra, the UV band exceeding the X-ray emission after and decaying slower than in the latter band (see Fig. 11, bottom left panel). The frequent interactions between shock and clouds broaden the velocity spectra of Ti and Ni considerably, and they lead to rather slow ejecta (see Fig. 11, bottom left panel).
We find a clear correlation between the photon emission and the cloud population (see Fig. 11, left panels). Without any clouds (top panel), the emission has a high hardness ratio between X-ray and UV. The most striking difference to the light curves of other models is probably the very slow decline of the UV emission after . The X-ray luminosity, on the other hand, decreases exponentially during this time.
Varying the density and temperature of the ISM
Since the density and temperature of the ISM surrounding the ejecta and the clouds are not well constrained by observations, we have varied these two parameters in our simulations as well. While the temperature of the ISM appears to have only a minor influence on the dynamics of the gas, increasing the gas density leads to remarkable differences from the reference model. In model S25Aw, we start the explosion in a medium of a temperature of instead of the standard value of . This choice does neither affect the time at which the SN shock wave hits the four main clouds, nor the emission from the gas, nor the velocity spectra of the different elements significantly. Thus, though we assumed a very cold ISM in all our simulations, our results are valid as well for a warm ISM.
A denser ISM, on the other hand, delays the expansion of the shock wave. We have simulated a model (S25Ad0.1 with a moderately increased ISM density of rather than . While the hydrodynamics of the expansion and the mixing are hardly affected by this change, we find a stronger emission in all bands (see middle panel of Fig. 12).
A further increase of the ambient density (models S25A-d1.0 and S25A-d10 with and , respectively) leads to a modified dynamics. The blast wave, sweeping up ambient matter much faster in this model, is decelerated more efficiently and reaches the four main clouds at a much later time: cloud NW at for model S25Ad1.0 and for model S25Ad1.0, respectively. Although the evolution is considerably slower than in the reference model, the mixing is similar to that of the reference model. We find that the tendency of stronger emission for increasing ISM density continues here: the luminosity of model S25A-d1.0 at the time of contact with cloud NW is , and the UV band shows a flat light curve exceeding during the entire time of the simulation. Because of the enhanced density in the interior regions of the SNR, the luminosity is dominated far less by the cloud than in the reference model, showing instead a comparably featureless volume-filling emission.
3.3.3 Distance of the SNR
The distance to Vela Jr. is a matter of intense debate, with estimates between and over . The models described above were performed assuming a low value for the distance, . Here, we show the results of a simulation where we placed the reference model at a distance of (model S25AD).
The distance enters our setup and the interpretation of our results in two ways:
Given the celestial coordinates of the clouds, their distance from the centre of our simulation grid are proportional to the assumed distance. Thus, if we place an explosion at a larger distance without changing its parameters like mass and energy, the time required until the blast wave interacts with the clouds, producing the current observational appearance, is considerably longer.
A larger distance means that the apparent brightness of the remnant in UV and X-rays corresponds to a higher absolute luminosity.
The interaction with the clouds starts with cloud SE at ; it takes about 400 yr more to reach cloud NW. For such a model, current appearance of Vela Jr. would be compatible with the supernova explosion around the year 0 AD. Overall this model shows the same dynamics as the reference model, and therefore the observational appearance at the time of the shock-cloud interaction is rather similar to the reference model (see Fig. 13, top left panel). However, after a much longer expansion, the energy and mass of the explosion are distributed over a much larger volume, and therefore, the absolute values of the emissivity is less than for the same model assuming a lower distance. Apart from this scaling, the emission properties are comparable. In particular, the model with a high explosion energy, exhibits a rather hard spectrum during the entire evolution we considered here (cf. the UV and X-ray lightcurves in the top right panel of Fig. 13). Additionally, as can be seen in the bottom panels of Fig. 13, the velocity spectra are very similar to our reference model, but the long time of expansion before the shock wave encounters the clouds lead to a somewhat less pronounced broadening of the population of fluid elements in velocity space.
4 Summary and conclusions
We have performed a series of hydrodynamic simulations of the expansion of a supernova explosion into an anisotropic ambient interstellar medium. In choosing the properties of the ISM, we had the clouds observed in the supernova remnant Vela Jr. in mind. The ISM in this SNR is dominated by four large clouds at roughly equal (projected) distance from the centre of the remnant. In addition to these clouds, we distributed a number of smaller, less dense clouds randomly in the vicinity of the explosion site.
The age of the SNR, and the mass and energy of the ejecta are only poorly constrained by observations. The goal of our simulations was to explore the parameter space of these properties of the explosion and to compare the simulation results to the UV and X-ray emission observed in the SNR (Aschenbach, 1998; Aschenbach et al., 1999; Tsunemi et al., 2000; Slane et al., 2001; Nichols & Slavin, 2004; Iyudin et al., 2005, 2007; Bamba et al., 2005; Nishikida et al., 2006; Pannuti et al., 2010; Iyudin et al., 2010; Pakhomov et al., 2012; Kim et al., 2012). Our main focus was on very energetic bipolar explosions that belong to the class of hypernovae. We initiated the explosion by placing the ejecta into a small spheroid (of roughly 1 pc size) at the centre of the computational domain. The explosion energy was mostly (95%) assumed to be in the form of kinetic energy. We chose the properties of the explosion (mass, energy, and axes ratio of the spheroid) according to the simulations by Maeda & Nomoto (2003), who computed the late phases of the hypernova explosions of main sequence stars of 25 and 40 . These authors provide masses and explosion energies as well as abundances of different elements in the ejecta.
We based the comparison of our simulation results with the observations on the following diagnostics:
The morphology of the ejecta and its interactions with the ISM is a direct result of the hydrodynamics of the explosion. In particular, we focus on the time scales of the expansion of the blast wave and on the development of small-scale structures.
At , we place a large number of test particles in the ejecta and assign different chemical compositions to them based on their initial positions. The particles are advected passively with the flow. The positions of these particles at later times allow us to gain a detailed description of the mixing of different elements in the ejecta caused by the decay of the flow into smaller structures. We construct time-dependent velocity spectra of oxygen, titanium and nickel that can be compared to observations.
We compute the thermal bremsstrahlung emission from the gas in different electromagnetic bands and construct mock light curves. Since the collision time of electrons and ions can be long in the tenuous gas forming the interior of the SNR compared to its evolutionary time scales, we have to account for the non-equilibrium ionisation (NEI) of the gas. Adding an equation for the energy of the electron gas to the system of hydrodynamic equations and modelling the coupling of the two species with a simplified ansatz for the collision time (Spitzer, 1962), we compute the electron temperature, which is, in general, different from the temperature of the ions. Based on these variables, we compute the photon emission due to bremsstrahlung in different bands at UV and X-ray energies (Tucker & Koren, 1971). We also look for correlations between the model parameters and the light curves and spectral properties of the models.
The main results of our simulations can be summarised as follows:
The supernova ejecta expanding into the ISM hits the four large clouds after a few hundred years. For our reference model, a hypernova-like bipolar explosion with an ejecta mass of and an energy of , this happens at roughly . Current observations of the SNR Vela Jr. indicate that this moment of first interaction between the blast wave and the clouds has occurred already. Though it is difficult to estimate the exact time that has passed since, we can conclude that the current appearance of the SNR would be consistent with a very energetic supernova explosion in the century, with an uncertainty of a few hundred years. The uncertainty of this result is rather large, as can be seen when varying the properties of the explosion and of the SNR environment. A lower explosion energy closer to the canonical value of or a higher density of the inter-cloud ISM delay the impact of the blast wave on the clouds significantly, and we would estimate a time of explosion about 2000 years earlier, consistent with the speculations by Katsuda et al. (2009). Besides the geometry of the SNR, we can estimate the rate at which the X-ray bright front marking the interaction of the blast wave and the clouds expands. The results of our models with a canonical explosion energy are in agreement with the value derived from observations by Katsuda et al. (2008), while a higher explosion energy gives rise to a faster expansion.
The emission of the SNR models is the result of a series of processes. Gas hit by the blast wave is heated, but at first the rise of the temperatures affects only the ions. Afterwards, collisions transfer energy to the electron gas, gradually establishing equilibrium between the two fluids. Since the collision frequency is highest in dense regions of the SNR, in particular the surfaces of the clouds, the electron temperature increases rather inhomogeneously, and, thus, the bremsstrahlung emission from the remnant is correlated very clearly to the presence and properties of the clouds. Typically, the early light curves of the models exhibit a steep decline in all energy bands as the shock wave expands rapidly, leading to cooling of the post-shock matter. Once the shock wave starts to interact with the clouds, the decline is much slower; some models even enter a plateau phase of essentially constant luminosity. During this phase, which corresponds most likely to the current observational appearance of RX J0852.0-4622, we estimate luminosities in the UV and X-ray bands between a few and a few , i.e., of the same order of magnitude as observations suggest. As noted above, the emission reflect the inhomogeneities of the ISM, yielding patchy emissivity maps in which the clouds are connected by filamentary structures, similar to what can be seen in observations. After the passage of the clouds by the shock wave, the effects of the interaction tend to become less prominent compared to the expansion, and the luminosities decrease again, though slower than in the initial stages. Models with high explosion energies per unit mass of the ejecta tend to produce an emission at the upper bound of luminosities given above; the same tendency can be observed if we increase the ISM density or add a large number of small clouds in the vicinity of the explosion. The light curves show pronounced differences across the range of photon energies. As the SNR expands and the mean temperatures drop, the emission gradually shifts to lower frequencies. High-energy bands fade away very quickly, while the bands we are most interested in, UV and low-energy X-ray, remain strong during the entire evolution. In most models, the X-ray band of 1 - 7 keV is more luminous than the UV band of 80-120 eV by a factor of a few. Exceptions from this behaviour can be found for spherical and low-energy explosions, as well as a very dense ISM with many clouds.
Using (passively advected) tracer particles we follow the evolution of three elements. Based on the results by Maeda & Nomoto (2003), we assume that titanium and nickel are present only in the bipolar outflow, while oxygen is distributed across all angles of the explosion. During the expansion of the ejecta, all elements are redistributed to lower velocities. Their velocity spectra broaden considerably when the flow of the post-shock matter past the clouds decays by instabilities into small scale structures, and particles are trapped inside the SNR at low velocities. Based on the evolution of the velocity spectra, we can distinguish a spherical explosion from a bipolar one: the latter exhibits a broader spectrum and shows nickel at much higher velocities.
The main caveat for a comparison of our models with observations is our limited modelling of the radiation processes. The observations of the Vela Jr. SNR show large contributions due to non-thermal processes and thermal line emission, which we currently do not include. We tried to assess the importance of this shortcoming by applying a criterion for X-ray synchrotron emission (Vink, 2012). We do indeed expect that our models would be strong sources of non-thermal synchrotron radiation. However, in regions where the shock hits the major clouds, X-ray synchrotron emission seems to be suppressed due to the low shock velocity unless the explosion was very energetic.
The main difficulty in identifying a possible progenitor star is that most of the observational properties of the SNR are affected by the characteristics of the ISM, in particular by the distribution of clumps and clouds, and by foreground absorption. The initial conditions of the explosion are washed out to a high degree during the interaction of the shock wave with the ISM.
Based on our results, we conclude that the current appearance of the Vela Jr. SNR is compatible with most of our simulations. Though the differences between our models do not allow us to determine the progenitor of the SN with certainty, we can at least point towards the most likely combinations of explosion properties and ages of the SNR. The assumption of a hypernova around the century (with a considerable uncertainty that may be up to a few hundred years ) and a less energetic explosion a few thousand years ago agree well with the observations. Hence, observations constraining either the explosion energy or the age of the remnant would constrain the other of the two parameters. Confronting results of our simulations with additional observations of the Vela Jr. SNR and of its environment in radio, optical, UV, and X- and gamma-rays can provide valuable information on the SNR progenitor type. Given the current set of models, observations identifying the thermal contribution of the emission would be most helpful, while the inclusion of non-thermal processes and line emission is the most important requirement for future modelling. It would also help us to better constrain the expansion rate and the age of the SNR, and its distance.
AFI and GFS were partially supported through the Grant of RF “11.G34.31.0076”. AFI acknowledges discussions and joint work with N. Chugai and Yu. V. Pakhomov related to the optical properties of Vela Jr., and thanks V. Burwitz, K. Dennerl and F. Haberl for their contribution to the X-ray imaging of RX J0852.0-4622, and especially Bernd Aschenbach for his unwaivering support of the whole effort to prove the young age of the remnant. MO acknowledges discussions with M.A. Aloy, D. Patnaude, R. Fesen, and D. Milisavlejic as well as support from the European Research Council (grant CAMAP-259276), and from the Spanish Ministerio de Ciencia e Innovación (grant AYA2010-21097-C03-01 Astrofísica Relativista Computacional). We are grateful for technical support by the system administrators of the Universitat de València, in particular C. Aloy. We thank the anonymous referee for his/her valuable comments.
- Aharonian et al. (2007) Aharonian F. et al., 2007, ApJ, 661, 236
- Aschenbach (1998) Aschenbach B., 1998, Nature, 396, 141
- Aschenbach et al. (1999) Aschenbach B., Iyudin A. F., Schönfelder V., 1999, A&A, 350, 997
- Bamba et al. (2005) Bamba A., Yamazaki R., Hiraga J. S., 2005, ApJ, 632, 294
- Berezhko & Völk (2010) Berezhko E. G., Völk H. J., 2010, A&A, 511, A34
- Bisnovatyi-Kogan et al. (1990) Bisnovatyi-Kogan G. S., Lozinskaia T. A., Silich S. A., 1990, Ap&SS, 166, 277
- Bisnovatyi-Kogan & Moiseenko (2008) Bisnovatyi-Kogan G. S., Moiseenko S. G., 2008, Progress of Theoretical Physics Supplement, 172, 145
- Bisnovatyi-Kogan & Silich (1995) Bisnovatyi-Kogan G. S., Silich S. A., 1995, Reviews of Modern Physics, 67, 661
- Cowie et al. (1981) Cowie L. L., McKee C. F., Ostriker J. P., 1981, ApJ, 247, 908
- Dubner et al. (1992) Dubner G., Giacani E., Cappa de Nicolau C., Reynoso E., 1992, A&AS, 96, 505
- Dubner et al. (1998) Dubner G. M., Green A. J., Goss W. M., Bock D. C.-J., Giacani E., 1998, AJ, 116, 813
- Fang & Zhang (2012) Fang J., Zhang L., 2012, MNRAS, 424, 2811
- Fesen et al. (2006) Fesen R. A., Hammell M. C., Morse J., Chevalier R. A., Borkowski K. J., Dopita M. A., Gerardy C. L., Lawrence S. S., Raymond J. C., van den Bergh S., 2006, ApJ, 645, 283
- Gaensler (1998) Gaensler B. M., 1998, ApJ, 493, 781
- Hammer et al. (2010) Hammer N. J., Janka H.-T., Müller E., 2010, ApJ, 714, 1371
- Inoue et al. (2009) Inoue T., Yamazaki R., Inutsuka S.-i., 2009, ApJ, 695, 825
- Inoue et al. (2010) Inoue T., Yamazaki R., Inutsuka S.-i., 2010, ApJ, 723, L108
- Inoue et al. (2012) Inoue T., Yamazaki R., Inutsuka S.-i., Fukui Y., 2012, ApJ, 744, 71
- Itoh & Fabian (1984) Itoh H., Fabian A. C., 1984, MNRAS, 208, 645
- Itoh & Masai (1989) Itoh H., Masai K., 1989, MNRAS, 236, 885
- Iyudin et al. (2005) Iyudin A. F., Aschenbach B., Becker W., Dennerl K., Haberl F., 2005, A&A, 429, 225
- Iyudin et al. (2007) Iyudin A. F., Aschenbach V., Burwitz V., Dennerl K., Freyberg M., Haberl F., Filipovic M., 2007, in ESA Special Publication Vol. 622 of ESA Special Publication, Multiwavelength Appearance of Vela Jr.: Is it up to Expectations?. p. 91
- Iyudin et al. (2010) Iyudin A. F., Pakhomov Y. V., Chugai N. N., Greiner J., Axelsson M., Larsson S., Ryabchikova T. A., 2010, A&A, 519, A86
- Janka et al. (2007) Janka H.-T., Langanke K., Marek A., Martínez-Pinedo G., Müller B., 2007, Phys. Rep., 442, 38
- Katsuda et al. (2008) Katsuda S., Tsunemi H., Mori K., 2008, ApJ, 678, L35
- Katsuda et al. (2009) Katsuda S., Tsunemi H., Mori K., 2009, Advances in Space Research, 43, 895
- Kesteven & Caswell (1987) Kesteven M. J., Caswell J. L., 1987, A&A, 183, 118
- Kifonidis et al. (2003) Kifonidis K., Plewa T., Janka H.-T., Müller E., 2003, A&A, 408, 621
- Kifonidis et al. (2006) Kifonidis K., Plewa T., Scheck L., Janka H., Müller E., 2006, A&A, 453, 661
- Kim et al. (2012) Kim I.-J., Seon K.-I., Min K.-W., Han W., Edelstein J., 2012, ApJ, 761, 135
- Kishishita et al. (2013) Kishishita T., Hiraga J., Uchiyama Y., 2013, A&A, 551, A132
- Klein et al. (1994) Klein R. I., McKee C. F., Colella P., 1994, ApJ, 420, 213
- Mac Low et al. (1994) Mac Low M.-M., McKee C. F., Klein R. I., Stone J. M., Norman M. L., 1994, ApJ, 433, 757
- Maeda & Nomoto (2003) Maeda K., Nomoto K., 2003, ApJ, 598, 1163
- Manchester (1987) Manchester R. N., 1987, A&A, 171, 205
- McKee & Cowie (1977) McKee C. F., Cowie L. L., 1977, ApJ, 215, 213
- Moriguchi et al. (2001) Moriguchi Y., Yamaguchi N., Onishi T., Mizuno A., Fukui Y., 2001, PASJ, 53, 1025
- Nakamura et al. (2006) Nakamura F., McKee C. F., Klein R. I., Fisher R. T., 2006, ApJS, 164, 477
- Nichols & Slavin (2004) Nichols J. S., Slavin J. D., 2004, ApJ, 610, 285
- Nishikida et al. (2006) Nishikida K., Edelstein J., Korpela E. J., Sankrit R., Feuerstein W. M., Min K. W., Shinn J.-H., Lee D.-H., Yuk I.-S., Jin H., Seon K.-I., 2006, ApJ, 644, L171
- Obergaulinger (2008) Obergaulinger M., 2008, PhD thesis, Technische Universität München
- Pakhomov et al. (2012) Pakhomov Y. V., Chugai N. N., Iyudin A. F., 2012, MNRAS, 424, 3145
- Pannuti et al. (2010) Pannuti T. G., Allen G. E., Filipović M. D., De Horta A., Stupar M., Agrawal R., 2010, ApJ, 721, 1492
- Park et al. (2004) Park S., Zhekov S. A., Burrows D. N., Garmire G. P., McCray R., 2004, ApJ, 610, 275
- Pettersson (2008) Pettersson B., 2008, Young Stars and Dust Clouds in Puppis and Vela. ASP Monograph Publications, p. 43
- Pittard (2011) Pittard J. M., 2011, MNRAS, 411, L41
- Pittard et al. (2010) Pittard J. M., Hartquist T. W., Falle S. A. E. G., 2010, MNRAS, 405, 821
- Pozzo et al. (2000) Pozzo M., Jeffries R. D., Naylor T., Totten E. J., Harmer S., Kenyon M., 2000, MNRAS, 313, L23
- Redman et al. (2002) Redman M. P., Meaburn J., Bryce M., Harman D. J., O’Brien T. J., 2002, MNRAS, 336, 1093
- Shimizu et al. (2012) Shimizu T., Masai K., Koyama K., 2012, PASJ, 64, 24
- Slane et al. (2001) Slane P., Hughes J. P., Edgar R. J., Plucinsky P. P., Miyata E., Tsunemi H., Aschenbach B., 2001, ApJ, 548, 814
- Spitzer (1962) Spitzer L., 1962, Physics of Fully Ionized Gases. New York: Interscience
- Suresh & Huynh (1997) Suresh A., Huynh H., 1997, \jcop, 136, 83
- Testori et al. (2006) Testori J. C., Arnal E. M., Morras R., Bajaja E., Pöppel W. G. L., Reich P., 2006, A&A, 458, 163
- Toro & Titarev (2006) Toro E. F., Titarev V. A., 2006, \jcop, 216, 403
- Tsunemi et al. (2000) Tsunemi H., Miyata E., Aschenbach B., Hiraga J., Akutsu D., 2000, PASJ, 52, 887
- Tucker & Koren (1971) Tucker W. H., Koren M., 1971, ApJ, 168, 283
- Vink (2012) Vink J., 2012, A&A Rev., 20, 49
- Wang et al. (2002) Wang L., Baade D., Höflich P., Wheeler J. C., 2002, The Messenger, 109, 47
- Woermann et al. (2001) Woermann B., Gaylard M. J., Otrupcek R., 2001, MNRAS, 325, 1213
- Wongwathanarat et al. (2010) Wongwathanarat A., Janka H., Müller E., 2010, ApJ, 725, L106
- Yamaguchi et al. (1999) Yamaguchi N., Mizuno N., Moriguchi Y., Yonekura Y., Mizuno A., Fukui Y., 1999, PASJ, 51, 765
- Zinchenko et al. (1995) Zinchenko I., Mattila K., Toriseva M., 1995, A&AS, 111, 95
- Zirakashvili & Ptuskin (2012) Zirakashvili V. N., Ptuskin V. S., 2012, Astroparticle Physics, 39, 12