We present 3D hydrodynamical simulations of a precessing jet propagating inside a supernova remnant (SNR) shell, particularly applied to the W50−SS433 system in a search for the origin of its peculiar elongated morphology. Several runs were carried out with different values for the mass-loss rate of the jet, the initial radius of the SNR, and the opening angle of the precession cone. We found that our models successfully reproduce the scale and morphology of W50 when the opening angle of the jets is set to 10° or if this angle linearly varies with time. For these models, more realistic runs were made considering that the remnant is expanding into an interstellar medium with an exponential density profile (as H i observations suggest). Taking into account all these ingredients, the large-scale morphology of the W50−SS433 system, including the asymmetry between the lobes (formed by the jet–SNR interaction), is well reproduced.
The standard supernova remnant (SNR) evolution theory, proposed by Woltjer (1970, 1972), predicts a spherical morphology for remnants expanding into a uniform interstellar medium (ISM). Observations show, however, that a spherical morphology is not the prevalent appearance of SNRs. Non-spherical morphologies can be produced by external or internal factors. Among the external factors is the existence of dense clouds or strong density gradients in the ISM where the SNR is expanding. Internally, in remnants resulting from a Type II supernova (Type II SN) explosion, a compact object (a neutron star or a pulsar for example) can remain inside the SNR, injecting enormous amounts of energy through its magnetic field and particle outflow, thus modifying the evolution of the host SNR.
Based on astronomical observations (mainly in radio continuum and X-ray), SNRs can be classified in three main morphological groups: (i) shell-type, when the remnant exhibits a shell morphology in both radio and X-ray images; (ii) plerions or Crab-like type, when the morphology of the remnant resembles the appearance of a centre-filled nebula, similar to the Crab's SNR (in this case the X-ray and radio emission mainly come from the central compact object and the SNR shell is not observed); and (iii) composite, when the SNR has a shell-type morphology plus a central component. The thermal-composite or mixed SNRs belong to this group; they are characterized by a radio shell with thermal X-ray emission in the interior.
There are remnants with more complex morphologies that can be included in group (iii); such is the case of W50 (Geldzahler, Pauls & Salter 1980; Downes, Pauls & Salter 1986; Elston & Baum 1987). The morphology of this remnant can be described as an almost spherical shell distorted by two symmetrical lobes along the east–west direction (see e.g. fig. 1a and 1b of Dubner et al. 1998). The source of relativistic jets SS433 (Abell & Margon 1979; Fabian & Rees 1979; Clark & Murdin 1978; Hjellming & Johnston 1982; Margon 1984; Hjellming & Johnston 1985) is located at the centre of W50. These jets have a speed of 0.26c and a precession period of 164 d and form a cone around the source with a half-angle of ∼20°. Analysis of the nature of the system reveals that the source of SS433 is a binary system, with an orbital period of 13 d; the jets' axis has an additional nodding motion with a period of 6 d (Newsom & Collins 1981, see also Collins & Scher 2002). The precession axis is aligned with the east–west direction of the W50 lobes. Several authors have given different explanations for this alignment. These explanations can be grouped in two possible scenarios: first, the SS433–W50 system is formed by the interaction of the jets of SS433 with the ISM (homogeneous or previously swept up by the stellar wind of SS433), without being affected by the SNR (Begelman et al. 1980; Königl 1983); and secondly, the eastern and western lobes of the system were formed by the interaction of the jets with the SNR shell (Zealey, Dopita & Malin 1980; Downes, Pauls & Salter 1986; Murata & Shibazaki 1996). Observationally, Elston & Baum (1987) found evidence supporting the latter idea. Later on, Dubner et al. (1998) also found evidence based on high-resolution radio-continuum images obtained with the Very Large Array (VLA). This was recently confirmed by the X-ray XMM–Newton observations of Brinkmann et al. (2007).
A recent work carried out by Durouchoux et al. (2001) (and references therein) examines SNR catalogues looking for other possible examples of SNRs with a similar elongated appearance. They identify three possible candidates: G309.2−00.6, G41.1−0.3 located in the 3C397 complex, and G54.1+0.3.
Following the jet–SNR interaction scenario, Velázquez & Raga (2000) carried out 2D axisymmetric numerical simulations to describe the W50−SS433 system, obtaining a rough agreement between their model and observations. Their model, however, was constrained by several hypothesis and assumptions.
In this work, we continue the study of Velázquez & Raga (2000), but remove several of their assumptions. We present 3D numerical simulations for the W50−SS433 system with the objective of providing a more realistic approach than the 2D model, following carefully the dynamics of the system as the precessing jets inject energy and momentum into the SNR cavity, collide with the SNR shock front, and change the spherical shape of the remnant. These simulations were carried out with the yguaz'u-a code (Raga, Navarro-González & Villagrán-Muniz; Raga et al. 2002).
This paper is organized as follows. Section 2 gives the initial conditions of the numerical simulations and a brief description of the yguaz'u-a code. The results are presented in Section 3. Finally, the discussion and conclusions of our work are given in Section 4.
2 INITIAL CONDITIONS AND THE yguaz'u-a CODE
As input data, we will use known values of the initial SN energy explosion, jet velocity, density, precession period and angle, etc., based on the only well-analysed observed case: the interaction between the relativistic jets from SS433 with the W50 SNR. Previous works have used 2D axisymmetric models of conical jets to describe this problem (Kochanek & Hawley 1990; Velázquez & Raga 2000), because the 164 d precession period of the jet (Margon et al. 1979; Hjellming & Johnston 1988, 1985, 1982) is considerably shorter than its dynamical evolution time, defined as tdyn=Rl/vj≃ 140 yr. Nevertheless, several assumptions on the geometry of the jet are needed in this type of model. 2D axisymmetric simulations model the symmetry axis as a reflecting boundary. This condition can produce undesirable effects such as an increase in the momentum flow along the symmetry axis (see e.g. Raga et al. 2007). Furthermore, in 2D axisymmetric simulations, it is not clear if the precession cone should be modelled by a filled conical jet or a hollow conical jet to represent the actual physical processes that take place in such a system. With a full 3D simulation, the assumption of axisymmetry is no longer needed and problems regarding the reflection condition disappear.
We have carried out numerical simulations of five models for precessing jets evolving into a SNR cavity. The main parameters for these models (the half-angle of precession α, the initial SNR radius RSNR and the jet's mass-injection rate ) appear in Table 1.
The numerical simulations were made employing the 3D version of the yguaz'u-a code (see Raga, Navarro-González & Villagrán-Muniz; Raga et al. 2002). The main characteristic of the code is that it has an adaptive grid and integrates the gas-dynamical equations with a second-order accuracy (in space and time) implementation of the ‘flux-vector splitting’ algorithm of van Leer (1982). Together with the gas-dynamical equations, a system of rate equations for atomic/ionic species is also calculated, which are then used to compute a non-equilibrium cooling function (a parametrized cooling rate being applied for the high temperature regime). They are also integrated with the gas-dynamical equations. Rate equations for H i–H ii, He i–He iii, C ii–C iv and O i–O iv are considered (for more details about reaction and cooling rates see Raga et al. 2002).
In all simulations, a five-level binary adaptive grid was used with a maximum resolution of 4.8 × 1017 cm. A Cartesian computational domain of 5123 pixels with a size of 80 pc along each axis was employed. The centre of the SNR and the source of the jets is located at the centre of the xy-plane with z= 0. Free outflow conditions were employed on each computational domain boundary.
The initial conditions for all models were adopted from different multiwavelength observational studies. For the jet, we adopt velocity vj= 7.5 × 109cm s−1, temperature Tj= 1 × 107K, radius 2.4 × 1018cm (equivalent to 5 pixels for the maximum resolution) and length 4.8 × 1018cm. For the ISM, we adopt density 0.01 cm−3 (Safi-Harb & Oegelman 1997), temperature 1 × 104K and velocity equal to zero; finally, the SNR was initialized employing the autosimilar Sedov solution (Sedov 1959) considering a remnant radius of 6.2 × 1019 cm or 20 pc (except for the model M3 that has half of this radius) and an SN explosion energy of 5 × 1050 erg.
2.1 The ISM density gradient
It is known that the system W50−SS433 is embedded in an ISM with a strong density gradient (Dubner et al. 1998), unlike the constant density ISM we have assumed so far. Thus, we made other runs of models M4 and M5 including a stratified ISM. In this case, the Sedov solution is not valid, so we divide the simulation in two parts. In the first one, the SN explosion is simulated considering that an initial energy of 5 × 1050erg (the same energy used in the previous models) is put into a sphere with a radius of 1.5 pc and a mass of 3 M⊙, uniformly distributed within this radius. Actually, the details of the SN explosion are not important if our interest is to study an evolved SNR, such as W50, because when the remnant enters into the adiabatic or Sedov evolution phase, the SNR has already lost the memory of the initial conditions. The only important parameter is the initial explosion energy E0. The SNR evolution is followed until the remnant shell reaches a radius of ≃ 21 pc, at an integration time of t= 6000 yr. In the second part, the jets are ‘turned on’ within the obtained SNR shell using the same parameters as in the previous models.
For these new runs, a computational domain of 512 × 512 × 1024 pixels was used, corresponding to a physical size of 80 × 80 × 160 pc. Both SNR and source jets were located at the centre of the computational domain. Free outflow conditions were applied on each grid boundary.
The density stratification is given by ρ(x, z) =ρ0 exp (−zp/hz), where zp=x sin β−z cos β, hz= 20 pc is a characteristic scalelength, and β= 20° is the inclination of the precession cone axis with respect to the z-axis.
As a first analysis, we study the interaction of the precessing jets with the shell of the SNR for α= 20°. With the objective of reproducing the observed morphology of W50, we explore different values for the jet's mass-injection rate (models M1 and M2) and the SNR radius (models M1 and M3). The range used for lies in the range inferred from X-ray observations (Safi-Harb & Oegelman 1997; Kotani et al. 2006). The time-delay between the SN explosion and the ‘turn-on’ of the jets is not well known; this motivated us to make a run with a shorter SNR radius to simulate a shorter time-delay (model M3).
Fig. 1 shows the temporal evolution of the density distribution for these models. The maps were obtained making a cut in the xz-plane passing through the centre of the jet's source. Kochanek & Hawley (1990) reported that hollow conical jets, expanding into a uniform ISM, experience a geometrical dilution of the momentum flux which prevents the jet's expansion. A similar situation occurs in models M1–M3. Even though we are not using the hollow conical jet approximation, we expect that the symmetric approximation will be valid for distances sufficiently close to the source. The asymmetry in the density distribution due to the jet's precession can be appreciated by zooming in on the central region; this is shown in Fig. 2 for model M1. The left-hand and right-hand panels display the xz and yz density distributions, respectively. The bottom panel is the point-to-point ratio between these two projections, and it clearly shows the asymmetry of the system, being larger for larger distances to the source. Even when the axial asymmetries are not very large, the figure shows the advantage of a 3D simulation over the usual assumption of axial symmetry.
Comparing models M1 and M2, it is clear that the increase in does not change the general structure of the cavity shell; the main difference is the larger cavity size for model M2. In the case of model M3, with RSNR= 10 pc, the resulting morphology is very elongated along the jet axis, and the SNR shell quickly loses its identity.
In these models, the jet is moving into a stratified medium (the SNR cavity) which could bend the initial jet's trajectory. However, a deflection of the jet is not observed. We conclude that a precessing jet with α= 20° moving into a remnant cavity cannot reproduce the observed shape of W50.
Because it was not possible to reproduce the morphology of W50 with models M1–M3, we ran a new simulation, model M4, with the same parameters as M1 but with α= 10°.
The left-hand panels of Fig. 3 display the temporal evolution of the density stratification for model M4. By t= 1000 yr, the SNR shell has been broken by the precessing jet. The evolution snapshots show that the jet has suffered a collimation process caused by the interaction of the jet's bow shock and the SNR shell. As a result, secondary shocks are formed which further reflect themselves on the z-axis. For t= 3000 yr, we can see that there is a rupture of the primary lobe present at t= 2000 yr. This pattern resembles the helical structure shown in the eastern lobe of W50. The result is also in agreement with the 2D axisymmetric modelling of Velázquez & Raga (2000).
3.1 The semi-aperture angle of the jets
These results point to an apparent contradiction between the observed semi-aperture angle α and our numerical models. How can they be made compatible?. Previous X-ray observations (Watson et al. 1983; Brinkmann, Aschenbach & Kawai 1996; Safi-Harb & Oegelman 1997) have shown that the angle subtended by the ‘X-ray lobes’ is smaller than expected from the precession cone of the innermost SS433 jets (40°). A possible solution is to consider that α is not constant during the jet evolution. We will use a model where α starts with a value equal to zero and then grows with time. We assume a simple model where α depends linearly on time: α=αmωt for 0 ≤t≤τc and α=αm for t > τc, where ω= 1/τc, τc is a free parameter that was chosen to be 2000 yr in order to produce lobes with the desired size, and αm= 20° is the estimated value of the opening angle today (Abell & Margon 1979).
This hypothesis is included in model M5 and the results are shown in the right-hand panels of Fig. 3. The snapshots are for the same evolution times as the left-hand panels. The shock front of the jet breaks the SNR shell as early as t= 1000 yr and for t= 3000 yr a well-formed lobe is observed. Unlike model M4, there is no secondary breaking of the lobe.
3.2 The ISM density gradient effect
The effects from the inclusion of the ISM density gradient to our model, described in Section 2.1, can be seen in Fig. 4. The left-hand panel shows the density distribution on the xz-plane for a run with the same parameters as for model M4, while the right-hand panel has the same parameters as model M5. These maps were obtained for an integration time of 3600 yr after the jet is turned on (corresponding to a total evolution time of 9600 yr). The ISM density gradient produces an extra asymmetry in the system which is clearly seen by the difference in the evolution of the lobes corresponding to the upper (western) and lower (eastern) jets. In the region where the ISM density is higher, the effect is to slow the SNR shell evolution and also to slow the penetration of the SNR shell by the jet's cocoon.
In order to have a more direct comparison with the actual observed system, we produced maps of the electronic density integrated along the line of sight for the new runs of models M4 and M5. These maps are shown in the left-hand and right-hand panels of Fig. 5, respectively. An inclination angle of 21° with respect to the plane of the sky has been taken into account. The angle was chosen according to the observed value given by Hjellming & Johnston (1988). The electronic density is a measure of the ionized zones in the system; it is expected that it will be enhanced just behind the shock fronts. For this reason, the electronic density is a good tracer of the structure observed in radio-continuum images. The physical sizes of the resulting distorted SNR shell are 56 and 103 pc, along the x- and z-directions, respectively. These physical sizes correspond to angular sizes of 1° 4′× 1° 58′ (if a distance of 3 kpc to SS433 is assumed), which are very similar to the observed dimensions (1°× 2°). The main difference between the panels of Figs 4 and 5 is that the more evolved lobe (lower) is more collimated for the case of model M4 than for M5.
Fig. 6 illustrates the comparison between simulations (left-hand panel of Fig. 5) and radio observations. The similarity with the system W50−SS433 is appealing. The eastern (left-hand panel) lobe in the simulation is, however, wider than the observed one. A possible cause for this difference is the presence of an H i cloud in the southern part of the eastern lobe (as suggested by Dubner et al. 1998), which could have contributed to an extra collimation of the lobe. Since this cloud is not included in our models, the result is a wider lobe. It is important to point out that an ideal comparison would be with synthetic radio maps generated directly from the simulations. However, our models do not include magnetic fields and thus it is not possible to produce such maps.
4 SUMMARY AND CONCLUSIONS
We have analysed the interaction between a precessing jet and a SNR shell using a 3D hydrodynamical simulation with the goal of studying the large-scale structure of the SS433–W50 system. The motivation for doing so is to put aside one of the usual assumptions made in the past regarding this kind of system, namely, using a cone (hollow or filled) as an approximation to model the jets from the source SS433. The precession period of the jets is lower than the dynamical time of the system; therefore, it would be expected that the precessing jets can be modelled by a single axially symmetric conical jet. Nevertheless, analysing our simulations, we found that there is a significant axial asymmetry, especially for material located far away from the central source (∼20–40 per cent of the ratio between the xz and yz density projections, see Fig. 2). These results confirm the advantages of using 3D over 2D axisymmetric simulations.
We found that the simple simulated model with a constant value for the aperture angle, as reported observationally for the jets in SS433: α= 20°, fails to reproduce the observed morphology of this system. This result holds in spite of exploring a large range of values for the mass-loss rate and the initial radius of the SNR (models M1–M3, see Fig. 1). This radius is connected with the time-delay between the SN explosion and the ‘turn-on’ of the jets. If this delay is too short, the jet's cocoon sweeps up the SNR shell erasing all traces of the initial spherical morphology (model M3, right-hand panel of Fig. 1).
If the value of α is arbitrarily reduced by half (model M4, see left panel of Fig. 3), then the shock front of the jets is collimated enough to break the SNR shell, forming two lobes very similar to those observed. This result is related to the observational fact that at a distance of ∼6 × 1019cm from the source centre, the angular extent of the X-ray jet is ∼18° (see, for instance, Brinkmann et al. 2007), far below the expected 40° of the original model. This suggests a re-collimation mechanism that is still unknown. It is no surprise then that the large-scale morphology of the system is reproduced well by our model when the opening angle is artificially reduced by half.
Such results led us to explore a different model where the value of α changes with time (model M5, see right-hand panel of Fig. 3), starting from zero and growing linearly until it reaches the constant value of α= 20°. This model also successfully reproduces the expected morphology. It is interesting to note that the corkscrew patterns observed in radio maps, at the tips of both eastern and western lobes of W50 (possibly produced by precession of the arcsecond-scale jet from SS433) are limited spatially by a cone which seems to be smaller than the derived ∼20° half-opening angle of the subarcsecond-scale jet. This is clearly illustrated in Fig. 1 of Moldowan et al. (2005), where a ∼20° half-opening cone is superimposed on the radio map of W50. The authors also suggest that the subarcsecond jet is playing a crucial role in determining the morphology of W50 nebula, since their Chandra observations of the inner region of the west lobe of W50 show a good correlation between the precession of the jet and the X-ray emission. Based on this evidence, we presented here a numerical exercise, in order to investigate the morphological effects of a time-varying opening angle jet on the whole structure of the nebula.
Finally, in order to give a more realistic model, two more simulations were made by using the parameters of models M4 and M5 but taking into account the observed density gradient of the ISM surrounding W50 (Dubner et al. 1998). In these runs, the Sedov solution for the SNR was not used, instead, the simulations were divided in two stages: in the first one, the SNR evolves alone; and the second part follows with the turn-on of the jets. The effect of the ISM density gradient is to produce different evolutionary patterns for the two lobes, making one of them more developed (see Fig. 4). Given the fact that in our description the magnetic field is not included, it is not possible to generate synthetic radio continuum maps; an alternative was to present maps of the electronic density integrated along the line of sight, which is a good tracer of the shock fronts. An inclination of the system with respect to the plane of the sky of ∼ 21° was considered to obtain this map (see Fig. 5).
In regard to the scale of the system, our models give physical sizes of 56 × 103 pc2, which are similar to the observed dimensions, if a distance of 3 kpc is considered for W50 (based on the H i study of Dubner et al. 1998). There are other distance estimations for this remnant in the literature. Recently, Lockman, Blundell & Goss (2007) give a distance between 5.5 and 6.5 kpc. This, however, does not change our main results, it would only be necessary to re-scale the simulated sizes and times.
In summary, the presented 3D models of the precessing jets of SS433 interacting with the shell of W50 reproduce the global characteristics (morphology and size) of this astrophysical object (see Fig. 6), giving a closer qualitative and quantitative agreement with observations. It is, however, important to say that the model is still far from being complete; it is of key importance to propose and test a model that can account for the early re-collimation of the jets. The magnetic field could be included in future works as a possible solution to this problem (this was analytically treated by Kochanek 1991, but for conical jets moving into an uniform ISM). Other details of the structure of the system, such as the reported bend in the eastern jet's propagation direction (see e.g. Brinkmann et al. 2007), can be studied by simulating a more complex ISM. Also, future works could focus on the late evolution of the jets, when they are far from the central regime.
The authors thank the referee for her/his comments, which helped us to improve the previous version of this manuscript. The authors acknowledge support from CONACyT grants 46628-F, and DGAPA-UNAM grant IN108207. JZ acknowledges CONACyT and DGEP (UNAM) scholarships, and support by the CAS Research Fellowship for International Young Researchers. The work of PFV was supported by the ‘Macroproyecto de Tecnologí as para la Universidad de la Información y la Computación’ (Secretarí a de Desarrollo Institucional de la UNAM, Programa Transdisciplinario en Investigación y Desarrollo para Facultades y Escuelas, Unidad de Apoyo a la Investigación en Facultades y Escuelas). GMD is thankful for financial support from grants PIP-CONICET(Argentina) 6433 and PICT-ANPCYT(Argentina) 14018. We thank Stan Kurtz (CRyA, Morelia, UNAM) for reading this manuscript. We also would like to thank the computational team of ICN: Enrique Palacios and Antonio Ramí rez, for maintaining and supporting our Linux servers, and Martí n Cruz for the assistance provided.