Continuous Jets and Backflow Models for the Formation of W50/SS 433 in Magnetohydrodynamics Simulations

The basic structure in all runs has a dynamic behavior similar to that found in previous light-jet simulations (e.g., Ohmura et al. 2020). Therefore, we briefly summarize the basic structure of light jets and then describe the overall morphology of each run in detail. Figure 1 shows maps of the gas density (top) and z-direction velocity (bottom) for run L-R1 at t = 88 kyr. The basic structures of supersonic jets, such as a beam, backflow, cocoon, and thick shell (shocked ISM) are seen. Surrounding plasma mixes into the cocoon owing to Kelvin–Helmholtz (KH) instabilities at the contact discontinuity between the jet gas and shocked ISM. In particular, backflows from the tips of the jet interact with each other, and vortex motions are highly developed in the cocoon. The background density spatially increases at z < −20 pc in the western direction, and the advance of the jet in the western direction thus slows appreciably. The magnetic field distribution and its detail are reported in Section 3.3.

Zoom In Zoom Out Reset image size

Figure 2 shows the shape of the jet for all runs when the contact discontinuity between the jet and ISM reaches 100 pc. Because the speed of advance of jets depends on the injection condition, the times at which the jets reaches 100 pc are different; i.e., H-R3, H-R1, L-R3, and L-R1 are at t = 16, 29, 55, and 88 kyr respectively. We here define the outer shape of jets as the contact discontinuity rather than a bow shock. However, it is difficult to identify the contact discontinuity between the jet plasma and shocked ISM because the gases are highly mixed by the KH instability. We therefore define the outer shape as the maximum radius of each z at which the jet tracer function f is less than 10⁻⁵. An analogous definition has been widely used in previous works on jet propagation (e.g., Gaibler et al. 2009). Here we denote the shell radius and the lengths from the origin to the eastern/western edges of the jet as l_R, l_east, and l_west.

Zoom In Zoom Out Reset image size

We find that all runs have a circular shell centered on the origin when jets reach ∼100 pc. The ratio of l_east to l_west is 100:75 for all runs. Lighter jets such as L-R1 and L-R3 form a larger shell. A lower density ratio corresponds to a lower lateral expansion velocity (see Equation (2)). However, because the kinetic energy luminosity of light jets is less than that of heavy jets, the advance of the light jet becomes much slower than that of the heavy jets. Lighter jets thus have enough time to form a larger circular shell while they propagate a distance of 100 pc. Run L-R1 therefore reproduces the characteristic morphology of W50, including a circular shell and elongated wings. Note that the wings of L-R3 are much wider than those of W50.

In this section, we describe the time evolution of the jet morphology (Figure 3). Top and bottom panels show L-R1 and L-R3, respectively. Colors denote the shapes at t = 8, 18, 28, 38, 48, 58, 68, 78, 88 kyr for run L-R1 and t = 5, 15, 25, 35, 45, 55 kyr for run L-R3. In the early period of t ≲ 30 kyr, the shapes of L-R1 and L-R3 are spheroids rather than circular shells. The spheroidal shape is well known from previous simulations and is similar to that of FR II radio sources, such as 3C 452 (e.g., Harwood et al. 2015). We find that a circular shell and wings then form. The region of the wing is filled with back-flowing plasma and it is thus difficult for the region to expand rapidly like the shell. Once the shell and wings have formed, they expand as self-similar evolution.

Zoom In Zoom Out Reset image size

We next consider what physical parameters determine the speed of advance of the jet. We easily estimate the speed of advance of the jet in the positive z-direction by adopting Equation (1). Meanwhile, the density distribution in the western direction exponentially increases, and the gradient effect should thus be included in the expression for the speed of advance:

$$\begin{eqnarray}&&{v}_{\mathrm{head}}=\displaystyle \frac{{dz}}{{dt}}=\displaystyle \frac{\sqrt{\eta (z)}}{1+\sqrt{\eta (z)}}{v}_{\mathrm{beam}},\,\eta (z)={\eta }_{0}\exp \left(\displaystyle \frac{z-20}{{Z}_{{\rm{d}}}}\right).\end{eqnarray} \tag{ 11 }$$

Here, we simply assume that A = 1 and ρ₀, v_beam are constant values.

Figure 4 (left) shows the positions of the tips of the eastern/western jets as a function of the time for each run (solid lines) and analytic lines (dashed lines η = 10⁻², dotted lines η = 10⁻³) of Equation (1) (z > 0) and Equation (11) (z < 0, toward to the Galactic plane). All of the lines at z > 0 (eastern) are not in good agreement with analytic values because the speed of advance decreases rapidly. The area of the jet head increases during the propagation of the jets; i.e., A is decreasing in Equation (1) because vortices grow downstream of the terminal shock. The jets, which have a large radius, have a small deceleration rate for both runs H and L because A remains at unity for a longer time than in the case of jets with a small radius. When a jet propagates toward the Galactic plane (in the western direction), which means the ambient density increases as the distance increases (z < 0), Equation (11) well describes the results of runs H-R3 and L-R3.

Zoom In Zoom Out Reset image size

We next describe the results of the time evolution of radial expansion for each run. Figure 4 (right) shows the shell radius l_R as a function of time for each run (solid lines). Dashed lines are the time evolution of the analytical model (Equation (2)) whose kinetic energy luminosity is L_kin,jet = 10³⁹, 10⁴⁰, and 10⁴¹ erg s⁻¹. All runs have l_R ∝ t^2/5 in accordance with the analytical solution. Note that Equation (2) describes the length of the bow shock along the R axis instead of the contact discontinuity. l_R is defined by the position of the contact discontinuity, and the expansion rate in all runs is thus lower than that for Equation (2). Note that the approximation for Equation (2) is not valid when the density ratio of the ISM to jets is close to unity. Furthermore, jets with larger radii have greater kinetic energy luminosity, and the speed of radial expansion thus depends on the jet's radius.

To discuss the formation mechanism of shells and cocoon, we focus on the length ratio of the two side jets and the length ratio of the jets and shell radius. In this work, the shell radius of W50 and the lengths from SS 433 to the tips of eastern and western wings are respectively R_W50 ∼ 48 pc, L_east ∼ 121.5 pc, and L_west ∼ 86.5 pc according to radio observations (Dubner et al. 1998; Goodall et al. 2011a). Thus, L_east/R_W50 and L_east/L_west are 2.5 and 1.4, respectively. Figure 5 (left) shows the time evolution of the length ratio l_east/l_west for all runs. The dashed and dotted lines are analytic solutions obtained when η = 10⁻² and 10⁻³. The numerical results for all runs and analytic solutions show a linear evolution over time. In this work, we perform calculations until the tip of one side jet reaches 100 pc, but the length ratio, l_east/l_west, increases with simulation time. It is thus expected that the length ratio becomes 1.4 when l_east = 120 pc, which is consistent with radio observations. Interestingly, the length ratios are roughly the same for all runs when the positions of jets are fixed. These results indicate that l_east/l_west is determined only by the density profile of the ISM.

Zoom In Zoom Out Reset image size

We next investigate the time evolution of the ratio of length from the center to the tip of the jet propagation in the eastern direction to the shell radius, l_east/l_R (right panel of Figure 5). For all runs, l_east/l_R saturates at certain values. The lower-density runs of L-R1 and L-R3 have a constant value of l_east/l_R for a long time. These results mean that outer shapes seem to undergo self-similar expansion in the later phase (see Figure 3). The ratio of l_east and l_R depends strongly on η. Although the speed of advance of the jet head is proportional to η, the speed of radial expansion does not depend on η. Although the speed of advance is strongly dependent on η, the speed of radial expansion is weakly dependent on η. Therefore, l_east/l_R is greater when η is closer to unity. The ratio of the length of the eastern wing to the shell radius of W50/SS 433 is L_east/R_W50 ∼ 2.5, which is consistent with the L-R1 run. In this work, we assume axisymmetric coordinates. We note that the speed of advance in three-dimensional simulations can be higher than that expected from axisymmetric simulations (Perucho et al. 2019).

Figure 6 (top) presents the distributions of physical quantities on the jet axis that show the characteristics of shock waves. The red and blue curves present the pressure $\bar{p}$ and velocity gradients ($-\partial {\bar{v}}_{z}/\partial z$) averaged according to

$$\begin{eqnarray}&&\bar{f}=\displaystyle \frac{1}{\pi {R}_{{\rm{jet}}}^{2}}{\int }_{0}^{{R}_{{\rm{jet}}}}2\pi {Rf}\,{dR}.\end{eqnarray} \tag{ 12 }$$

The velocity gradient in the direction of propagation roughly represents the magnitude of the shock wave. The bottom panel of Figure 6 shows streamlines of the average velocity. The velocity gradient (blue curve) has many sharp spikes while the pressure distribution has wide-tailed peaks. The two are generally correlated and the peaks are considered as shock waves, such as reconfinement shocks, oblique shocks, and terminal shocks. As an example, we see reconfinement shocks at z ∼ −10, 10, and 20 pc. After the jet propagates far away, it is difficult to identify the origin of each peak because many reflection waves propagate in various directions. However, when including the density and velocity distributions (Figure 1), we can identify the terminal shocks (z ∼ −70, 90 pc) and oblique shocks (z ∼ −45, 60 pc) in both side jets. A terminal shock is a location of energy exchange from jet kinetic energy to thermal energy and it produces backflows. When the propagation speed decreases, the backflow becomes turbulent and has circular and arc-like motion (see the bottom panel of Figure 6). The backflow therefore interacts with the beam itself and drives the strong oblique shock (Mizuta et al. 2004). We estimate the size of the backflow vortex in Section 4.1. In our calculation, the oblique shock in the eastern jet releases more kinetic energy than the terminal shock. Meanwhile, the western jet shows the opposite correlation. The background ISM density increases in the western direction, and the conversion efficiency of the jet kinetic energy to thermal energy at the terminal shock in the western direction is therefore higher than that in the eastern direction. The red curve shows that the pressure maximum of the jet head at z = −75 pc is three times that of the jet in the eastern direction (z = 100 pc).

Zoom In Zoom Out Reset image size

Figure 7 (left) shows the distribution of toroidal magnetic fields and magnetic energy for run L-R1 at t = 88 kyr (left) and the distribution of current density, j = ∇ × B (right). The magnetic energy is increased by shock compression; in particular, its energy is strongly enhanced downstream of the terminal shock of the western jet. We find three areas in the cocoon in Figure 7 (top left): an area of a positive field (blue), an area of a negative field (red), and an area of a weak magnetic field (white). The two backflows, which are antiparallel toroidal fields, interact in the cocoon, and the gas is mixed by vortices. The magnetic fields are therefore dissipated by magnetic reconnection in this mixing region. Figure 7 (right) shows that the interaction of the two backflows generates current sheets, where particles can be efficiently accelerated by magnetic reconnection. Our results would be affected by an artifact of the axisymmetric condition. If we conduct three-dimensional simulations, turbulence would be generated in the cocoon, which would readily convert toroidal and poloidal fields.

Zoom In Zoom Out Reset image size

We here make a simple estimation of the interval of the oblique shocks. Note that the oblique shocks are formed by the interaction between the beam and backflow. The interval thus corresponds to the size of the backflow vortex l_bf. For the eastern jet of L-R1, we assume that the radius r_bf is about 10 pc (see Figure 3). To verify this assumption, we make the simple argument of mass flux conservation. The simulation result gives a z-direction velocity and density of the backflow of about 0.15v_jet and 10⁻¹ ρ_jet, respectively. Conservation of mass flux in a hollow cylindrical symmetry gives $\pi {R}_{\mathrm{jet}}^{2}{\rho }_{\mathrm{jet}}{v}_{\mathrm{jet}}\,=\pi ({r}_{\mathrm{bf}}^{2}-{r}_{\mathrm{jet}}^{2}){\rho }_{\mathrm{bf}}{v}_{\mathrm{bf}}$. We obtain r_bf ∼ 8R_jet = 8 pc, and the assumption is reasonable. We assume that the R-direction velocity is about the speed of sound c_bf. To obtain c_bf, we first calculate the speed of sound at the hotspot c_h using the Rankine–Hugoniot condition,

$$\begin{eqnarray}&&{c}_{{\rm{h}}}^{2}=\displaystyle \frac{5{{ \mathcal M }}^{4}+14{{ \mathcal M }}^{2}-3}{16{{ \mathcal M }}^{2}}{c}_{\mathrm{jet}}^{2},\end{eqnarray} \tag{ 13 }$$

where ${ \mathcal M }=10$ is the Mach number of the jets while ${c}_{\mathrm{jet}}={{ \mathcal M }}^{-1}{v}_{\mathrm{jet}}$ is the speed of sound of the jets. Note that we ignore the speed of advance because it is much lower than v_jet and c_h. We can estimate c_bf simply using the adiabatic condition:

$$\begin{eqnarray}&&{c}_{\mathrm{bf}}^{2}={\left(\displaystyle \frac{{R}_{\mathrm{jet}}^{2}}{{r}_{\mathrm{bf}}^{2}}\right)}^{\gamma -1}{c}_{{\rm{h}}}^{2},\end{eqnarray} \tag{ 14 }$$

where we assume that the radius of the hotspot is the same as the radius of the jet. The R-direction velocity is therefore evaluated as

$$\begin{eqnarray}&&{c}_{\mathrm{bf}}={v}_{\mathrm{jet}}\sqrt{{\left(\displaystyle \frac{{R}_{\mathrm{jet}}^{2}}{{r}_{\mathrm{bf}}^{2}}\right)}^{\gamma -1}\displaystyle \frac{5{{ \mathcal M }}^{4}+14{{ \mathcal M }}^{2}-3}{16{{ \mathcal M }}^{4}}}\sim 0.12{v}_{\mathrm{jet}}.\end{eqnarray} \tag{ 15 }$$

This value is close to the z-direction velocity of backflow v_bf, which is about 0.15v_jet (see Figure 1). The oscillating time of the R-direction backflow is τ = r_bf/c_bf, and the interval of oblique shocks is thus evaluated as

$$\begin{eqnarray}&&{l}_{\mathrm{bf}}=2\tau {v}_{\mathrm{bf}}\sim 25\,\mathrm{pc}.\end{eqnarray} \tag{ 16 }$$

This value is roughly consistent with our numerical result. The bottom panel of Figure 6 shows that the backflow departs from z = 90 pc and arrives at z = 65 pc, where a strong oblique shock is excited. However, the actual dynamics of the beam and backflow are more complex and many weak shocks are thus excited along the beam.

We create a pseudo-synchrotron image of W50/SS 433 using physical values from MHD simulations. To calculate the synchrotron emissivity, we convert the two-dimensional axisymmetric cylindrical coordinates to three-dimensional Cartesian coordinates and assume a viewing angle of 788. We calculate and integrate the emissivity at each cell along the line of sight direction:

$$\begin{eqnarray}&&{ \mathcal I }={\int }_{\mathrm{LoS}}{{ \mathcal J }}_{\mathrm{sync}}{ds},\end{eqnarray} \tag{ 17 }$$

where ${ \mathcal I }$ and ${{ \mathcal J }}_{\mathrm{sync}}$ are respectively the synchrotron intensity and emissivity. Our simulations do not model the evolution of nonthermal electrons, and hence some assumptions are needed to calculate the synchrotron emissivity. One is that the spectral index is uniform everywhere. The other is that the number density of nonthermal electrons is proportional to the gas density and thermal energy density. Under these assumptions, the synchrotron emissivity is given by (Jun & Norman 1996)

$$\begin{eqnarray}&&{{ \mathcal J }}_{\mathrm{sync}}={C}_{1}(\nu ){\rho }^{1-2\alpha }{p}^{2\alpha }{B}_{\perp }^{\alpha +1},\end{eqnarray} \tag{ 18 }$$

where C₁, B_⊥, ν, and α are a normalized constant, the magnetic field component in the plane of the sky, the frequency of radiation, and the spectral index, respectively. We adopt α = 0.52 because some radio observations provide the flatter spectrum in the eastern wing of W50 (Dubner et al. 1998; Gao et al. 2011). The calculation of polarized emissions is beyond the scope of this paper due to axisymmetry.

Figure 8 (left) shows the synchrotron intensity map at t = 88 kyr for run L-R1. The synchrotron power depends strongly on magnetic energy, and a weak emission is thus observed along with the shell (see Figure 7). The bright region in the western wing is around the oblique shocks and the termination shock. Small bright spots on the jet axis are formed by oblique shocks, whose spacing is defined by Equation (16).

Zoom In Zoom Out Reset image size

We can see the wide bright area both in the west and east wings, which consist of thin filaments. The accumulated magnetic field along the contact discontinuity is the origin of these bright regions. We consider that these bright spots correspond to X-ray hotspots of SS 433. For example, the magnetic filament around (R, z) = (10, 80) is the origin of the bright filament closest to the eastern end. On the other hand, because the magnetic field around the terminal shock is diffused, there is no X-ray radiation around that region.

The synchrotron intensity map obtained from numerical simulation is similar to X-ray observations of W50 (Brinkmann et al. 1996, 2007; Safi-Harb & Ögelman 1997). In particular, the positions of X-ray hotspots of W50 (e1, e2, w1, and w2) are consistent with those in the numerical results. The oblique shocks are thus thought to be rather efficient particle accelerators, in contrast with the terminal shocks in the long-term evolution of jets. Sudoh et al. (2020) and Kimura et al. (2020) constructed an analytic model for nonthermal leptonic and hadronic emission from knots in SS 433 jets, and their results suggest that the X-ray observation can be reproduced by synchrotron emission when the acceleration process is efficient.

We next discuss the distribution of electrons emitting radio waves in W50/SS 433. One main difference between microquasar jets and AGN jets is the outburst age, i.e., the age of relativistic electrons in the jet cocoon. The synchrotron cooling time t_c of nonthermal electrons in the microquasar jets is ${t}_{{\rm{c}}}\sim 1.9\times {10}^{13}{B}^{-2}{E}_{{\rm{e}}}^{-1}\,{\rm{s}}$, where E_e is the energy of electrons. The magnetic field strength of the eastern and western wings is about 1–100 μG (Safi-Harb & Ögelman 1997) and the cooling timescale is thus estimated as t_c ∼ 100 kyr when the frequency is 1 GHz and the magnetic field is 100 μG. The plasma created around the terminal or oblique shock is thus not cooled to a few GeV and emits radiation at a centimeter wavelength. We calculate the shocked particle content in W50/SS 433, ${ \mathcal F }$, to integrate the tracer function, f:

$$\begin{eqnarray}&&{ \mathcal F }={\int }_{{\rm{LoS}}}f\,{ds}.\end{eqnarray} \tag{ 19 }$$

Figure 8 (right) is the same as Figure 8 (left) but we display a map of the distribution of jet particles. The plasma accumulates in the cocoon after passing the terminal shock and oblique shocks. The lifetime of electrons for synchrotron radiation is much longer than the jet active time (simulation time), and relativistic electrons in the shell can thus emit continually. The intensity is lower in both wings than in the shell because relativistic electrons are advected by backflow so that there are fewer electrons in both wings. This characteristic is also found for the radio emission of the eastern wing. Another radio characteristic of the eastern wing is an extended radio filament in the terminal region. In addition, Brinkmann et al. (2007) found an X-ray ring in this region, which corresponds spatially to the terminal shock of the SS 433 jet. However, our numerical result indicates that the shock size is only a few parsecs. The filament is thus formed by the backflow of the terminal shock. Additionally, Sakemi et al. (2018a) reported the existence of helical fields that coil the eastern wings. This result is consistent with our backflow scenario. We next discuss the shell emissions of W50. We find that magnetic reconnection generates current sheets in the shell (see Figure 7), and reaccelerated electrons in current sheets provide prominent radio emission. This radio emission corresponds to the filament structure in the southern region of W50; otherwise, the MeV γ-ray emission observed in the northern region of W50 is explained by this acceleration scenario. In this simulation, we assume an axisymmetric condition and the magnetic field in the shell of W50 is thus exhausted by magnetic reconnection. However, if we carry out three-dimensional simulations, the magnetic field can be stored in the shell because the flow can escape easily in any direction. The synchrotron emission from the shell would therefore be enhanced.

Previous studies have suggested that the SNR produced by the compact object SS 433 formed the shell structure of W50, and we thus discuss the difference in radial expansion between our backflow scenario and SNR scenarios. We next compare the features of radio emissions of W50/SS 433 with those of other SNRs. We mention again that the shell diameter is D = 96 pc. The typical evolution of an SNR is described in the Sedov phase (i.e., the adiabatic phase). An expanding SNR is then described by the snowplow phase when radiative losses become important. We assume that the shell of W50 is in the Sedov phase and radiative cooling is thus negligible. Note that the radial expansion of the Sedov phase is faster than that of the snowplow phase. The shell radius of the SNR develops in the Sedov phase as $R\propto {({E}_{\mathrm{SNR}}{t}^{2}/{\rho }_{\mathrm{ISM}})}^{1/5}$, where E_SNR is the initial explosion energy of the SNR. The dashed lines in Figure 9 show the radius of the SNR as a function of time for three different values of the ISM density. We simply assume that the ISM has a constant density ρ_ISM = 1m_p, 0.1m_p, and 0.01m_p g cm⁻³ and E_SNR,51 = 10⁵¹ erg. When the ISM density is greater than 0.1m_p, it takes more than 100 kyr to reach a shell radius of 48 pc. The ISM is thus about 0.01m_p, which is a very low density with which to form the W50/SS 433 shell. However, we note that W50/SS 433 is located near the Galactic plane, and H i emissions are observed around W50 (Peek et al. 2011; Sakemi et al. 2021). Such a low-density environment is thus unlikely. In another case, we believe that the W50 nebula formed in a high-power explosion (E_SNR ≫ 10⁵¹ erg), such as a hypernova or multiple supernovas. In contrast, we plot the time evolution of the radial expansion for the light jets for jet kinetic energy luminosity L_k,39 = 10³⁹ erg s⁻¹ in Figure 9. The approximate solution of the radial expansion of jets is R_jet ∝ t^3/5. We mention again that the Sedov solution is proportional to t^2/5. The radius obtained using the jet model becomes larger than that obtained using the SNR model when the total energy exceeds the supernova explosion energy, E_SNR = 10⁵¹ erg. The time is about t ∼ 40 kyr for these model settings. In the case of jets instead of an SNR, the shell of W50 reaches 48 pc within 100 kyr when the ISM is 0.1m_p g cm⁻³. Therefore, jets having long-term activity, for which a lifetime of >10 kyr is expected, are a good candidate for the formation mechanism of W50/SS 433.

Zoom In Zoom Out Reset image size

We next discuss radio observations of typical SNRs and W50/SS 433. Galactic SNRs have a radio observational relationship between the radio surface brightness and the diameter of the SNR, the so-called radio surface-brightness-to-diameter relationship (Σ_R–D) (e.g., Pavlovic et al. 2014). Figure 10 displays the Σ_R–D relationship at 1 GHz for Galactic SNRs (blue circles and green diamonds) and W50/SS 433 (red star). The figure shows that W50/SS 433 is prominent and has a large size compared with Galactic SNRs. We note the observations of four prominent Galactic SNRs (CTB 37A, Kes97, CTB 37B, and G65.1+0.6) interacting with molecular clouds, which enhances radio emissions. The wings of W50 have a radio surface brightness comparable to that of the shell of W50 (Dubner et al. 1998). Previous hydrodynamical simulations clearly showed that the jetted gas and the gas-related supernova explosions are separated; in particular, the jetted gas is not distributed around the shell region (e.g., Goodall et al. 2011a). Although some observations suggest that there are molecular clouds associated with jets in the western wing of W50 (Yamamoto et al. 2008; Liu et al. 2020), there is no critical evidence of interacting molecular clouds, which would accelerate nonthermal particles. Additionally, the jet–cloud interaction is not strong and it thus does not enhance radio emissions throughout the western wing. Therefore, the origin of the prominent emission from the shell of W50 is not the same as that of other prominent SNRs.

Zoom In Zoom Out Reset image size

Finally, we mention that the shell of W50 maintains a nearly perfectly circular structure. Stafford et al. (2019) investigated the relationship between the radio morphology and the sizes (ages) of SNRs. They found that larger (older) SNRs are more elliptical and/or elongated because the inhomogeneous ISM, turbulent structure, and molecular clouds affect their morphology. As we mentioned above, W50 is in a large size (old) category if its shell comprises SNRs.

In this section, we compare our results with those of Goodall et al. (2011a). We first mention that the radial expansion of the SNR is too fast in the simulations conducted by Goodall et al. (2011a) compared with the Sedov solutions and hydrodynamical simulations of SNRs. They assumed that the kinetic energy of SNRs is ∼10⁵¹ erg and the ISM density is ∼0.1 cm⁻³, and the shell radius of SNRs reaches 50 pc within 20 kyr. They might have underestimated the normalization of the kinetic energy of SNRs by an order of magnitude. We also mention the speed of advance of jets in Goodall et al. (2011a). As an example, model Jet1 in Goodall et al. (2011a) has a young age of 2.2 kyr when the jet crosses the tip of the eastern wing. The average speed of advance of model Jet1 almost maintains the launch velocity of 0.26c from SS 433, i.e., 120 pc/2.2 kyr ∼ 0.18c. Thus, the advance of model Jet1 does not decelerate. Note that the results of kinematical studies of the radio observations of W50 indicate that an upper limit to the proper motion of the radio filaments of the eastern wing is 0.0405c (Goodall et al. 2011b). Jets that do not decelerate must be comparable to or greater than the ISM in density. Some previous numerical simulations have modeled protostar jets for heavy jets (e.g., Ustamujic et al. 2016). However, numerical and observational studies of AGN jets and microquasar jets suggest that extensive cocoons develop when the surrounding medium is lighter than the medium of the jets.

Panferov (2017) constructed a model combined with the theory of the dynamics of SNRs, wind bubbles, and jets for the formation of W50. Their jet model was based on turbulent jets in Landau & Lifshitz (1987). In the case of a wind bubble, if the binary system of SS 433 is a super-Eddington accretion, the age of the wind bubble needs to exceed 300 kyr to produce the large shell of W50. They also modeled the dynamics of SNRs in the Sedov phase and snowplow phase, describing the dynamics using the fitting function of Kim & Ostriker (2015). They argued that the age of W50 as an SNR is at least 97 kyr, which agrees with our results in Figure 9. The main result of Panferov (2017) is that the shell of W50 formed in a supernova explosion about 100 kyr ago, and the SS 433 jets produced both wings and pushed up the pre-existing SNR within 27 kyr. However, hydrodynamical studies showed that interaction between the jets and SNR does not push up the shell of W50 in the radial direction (Goodall et al. 2011a). Therefore, some other mechanism is needed for the model of Panferov (2017) to exchange energy efficiently between the jets and SNR shell.