Continuous Jets and Backflow Models for the Formation of W50/SS433 in Magnetohydrodynamics Simulations
Takumi Ohmura, Kojiro Ono, Haruka Sakemi, Yuta Tashima, Rikuto Omae, Mami Machida
DDraft version February 16, 2021
Typeset using L A TEX twocolumn style in AASTeX62
Continuous Jets and Backflow Models for the Formation of W50/SS433 in Magnetohydrodynamics Simulations
T. Ohmura, K. Ono, H. Sakemi, Y. Tashima, R. Omae, and M. Machida Department of Physics, Faculty of Sciences, Kyushu University, 744 Motooka, Nishi-ku, Fukuoka 819-0395, Japan Division of Science, National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan (Received 2020 November 23; Revised 2021 January 29; Accepted 2021 February 10)
Submitted to ApJABSTRACTThe formation mechanism of the W50/SS433 complex has long been a mystery. We propose a newscenario in which the SS433 jets themselves form the W50/SS433 system. We carry out magneto-hydrodynamics simulations of two-side jet propagation using the public code CANS+. As found inprevious jet studies, when the propagating jet is lighter than the surrounding medium, the shockedplasma flows back from the jet tip to the core. We find that the morphology of light jets is spheroidalat early times, and afterward, the shell and wings are developed by the broadening spherical cocoon.The morphology strongly depends on the density ratio of the injected jet to the surrounding medium.Meanwhile, the ratio of the lengths of the two-side jets depends only on the density profile of thesurrounding medium. We also find that most of the jet kinetic energy is dissipated at the obliqueshock formed by the interaction between the backflow and beam flow, rather than at the jet terminalshock. The position of the oblique shock is spatially consistent with the X-ray and TeV gamma-rayhotspots of W50.
Keywords:
ISM: individual objects (W50/SS433) – methods: numerical – magnetohydrodynamics(MHD) – galaxies jets INTRODUCTIONAstrophysical jets emitted from gravitational objectssuch as black holes, neutron stars, and protostars areuniversal phenomena observed in various layers of theMilky Way. The jets interact with the surrounding in-terstellar medium (ISM), and the kinetic energy of thejets is therefore transformed into ISM thermal and tur-bulent energies. Additionally, some of the kinetic en-ergy of jets is given to energetic, non-thermal particles.Microquasar jets are therefore a candidate for cosmic-ray acceleration sites (e.g., Heinz & Sunyaev 2002; Bor-das et al. 2009). In particular, Cooper et al. (2020)argued that the cosmic-ray component generated in mi-croquasar jets makes a non-negligible contribution toobserved cosmic-ray spectra. However, physical quanti-ties and characteristics of galactic jets are currently notwell understood because the angular size of many galac-
Corresponding author: T. [email protected] tic jets is small compared with the spatial resolution ofobservations.The radio nebula W50 is located near the galacticplane, and it is a rare galactic object for which high-quality radio observations have been made (Dubneret al. 1998; Gao et al. 2011; Farnes et al. 2017; Broder-ick et al. 2018). The radio morphology of W50 has twocharacteristics, namely a circular shell with a diameterof 58 arcmin and lateral extensions, which are calledwings, in eastwest directions. The X-ray binary SS433launches sub-relativistic precessing jets at the center ofthe shell structure of W50 (Blundell & Bowler 2004;Roberts et al. 2008; Blundell et al. 2018). The jet ve-locity of SS433 is v jet = 0 . c , where c is the speedof light, and the precession period is 164 days (Margon1984; Margon & Anderson 1989). Additionally, the ki-netic energy luminosity of jets is L kin , jet ∼ erg s − (Brinkmann et al. 2005). The wings are usually associ-ated with the jets because the major axis of W50 is wellaligned with jets from SS433. Polarization observationsindicate the existence of a helical magnetic field in theeastern wing (Farnes et al. 2017; Sakemi et al. 2018a,b). a r X i v : . [ a s t r o - ph . H E ] F e b Ohmura et al.
Non-thermal X-ray and gamma-ray observations indi-cate that the W50/SS433 system is an efficient site forcosmic-ray acceleration, and the emissions have beendetected along the major axis of the jets (Brinkmannet al. 1996a,b, 2007; Safi-Harb & ¨Ogelman 1997). Sev-eral prominent regions have been identified as X-rayhotspots and are labeled e1, e2, and e3 on the easternside and w1 and w2 on the western side. In particular,very-high-energy γ -rays ( >
25 TeV) have been detectedfrom around e1, e2, and w1 hotspots by the High Al-titude Water Cherenkov observatory (Abeysekara et al.2018). The spectral energy distribution at the hotspotsis consistent with the model for the leptonic scenario.This scenario is that relativistic electrons scatter cos-mic microwave background photons to TeV energies viathe inverse Compton process. In addition, no TeV γ -ray signal from the central binary system has been de-tected. The hotspots should therefore be in-situ particleacceleration sites associated with the SS433 jets. More-over, several studies provided observational evidence of γ -ray emissions at MeV to GeV energy levels from notonly around X-ray hotspots but also around the spheri-cal shell of W50 (e.g., Xing et al. 2019; Sun et al. 2019).Detailed analysis of over 10 years of GeV gamma-raydata recorded by the Fermi Gamma-ray Space Telescoperevealed that the precessional period of the SS433 jetsis linked with the time variability of γ -ray excess in thenorthern shell of W50 (Li et al. 2020). These observa-tions also strongly support that W50/SS433 has plau-sible cosmic-ray acceleration sites. However, it remainsunclear whether the outstanding acceleration is due tothe interaction between the jet and the supernova rem-nant (SNR) or just the SS433 jets themselves.W50 has a unique radio morphology, called themanatee nebula, and the formation mechanism of theW50/SS433 system has been debated for many years.One scenario is that W50/SS433 formed through in-teraction between the jets of SS433 and wind bubbles,which were swept up by the stellar wind of the compan-ion star of SS433 (Begelman et al. 1980; Konigl 1983).In addition, Asahina et al. (2014) suggested that theshell of W50 was formed by winds from an accretiondisk with a mass accretion rate exceeding the Edding-ton limit. These scenarios pose the problem that it takesa long time for a wind explosion to reach 40–50 pc. Themost popular scenario is that the jets and SNR shell in-teract. The circular shell of W50 is the shell of an SNR,resulting from a supernova explosion at SS433’s birth.The interaction between the jets and SNR shell form the elongated two-side wings . A structure protrudingfrom a supernova shell is often found in core-collapsedSNRs and looks like the elongated structure of W50(e.g., Gaensler et al. 1998). Grichener & Soker (2017)estimated the additional energy used to create the wingsof the core-collapse SNR to be 1–10 percent of the explo-sion energy of the supernova under simple geometricalassumptions. However, these wings are shorter thanthe radius of their own SNR. In contrast, the length ofthe eastern wing of W50 is much longer than the radiusof the shell. For the wing length to be comparable tothe shell radius, the additional kinetic energy neededto extend the W50 wing must be the same as the su-pernova explosion energy. Simply assuming the kineticenergy luminosity of the SS433 jets is 10 erg s − andthe explosion energy is 10 erg, the long-term activityof SS433 for 10–100 kyrs is essential.The interaction between the SS433 jets and the SNRshell has been studied by conducting hydrodynamicsimulations. Vel´azquez & Raga (2000) conducted thehydrodynamic simulation of jet propagation inside anSNR, as modeled using the analytical Sedov solution.Additionally, they showed that the cocoon gas, which isheated at the jet’s terminal shock, slightly increases theshell expansion velocity near the interaction region. Athree-dimensional simulation conducted by Zavala et al.(2008) revealed that the exponential density profile of anambient ISM identified by HI observations reproducedthe length ratio of the easternwestern wings. The rea-son for the difference in the wing lengths is that a jetpropagating into a region of increasing density drasti-cally decelerates. Goodall et al. (2011a) reported two-dimensional hydrodynamical simulation results for mod-els of several jets propagating in an SNR shell; i.e., pre-cession jets, cone jets, and fireball jets. They suggestedthat multi-episodic jet activity is needed for the forma-tion of W50. Recently, Millas et al. (2019) performed ahigh-resolution adaptive-mesh refinement hydrodynam-ical simulation of precessing jets propagating in an SNR.Their results indicate that precession does not affect thelarge-scale morphology of W50, while the jets are wellcollimated at ∼ .
068 pc from SS433 and become hol-low (Monceau-Baroux et al. 2015; Millas et al. 2019).Note that Monceau-Baroux et al. (2014) reported thatprecession would be dynamically important from theSS433 scale to the sub-parsec scale because the precess-ing jets transport more kinetic energy into the ISM thando straight jets owing to the increase in the area of inter-action between the jets and ISM. These hydrodynamical Some papers refer to ’wings’ as ’ears’. models for the interaction of the jets and SNR can ex-plain the morphology of W50 within a realistic time,which is 10–100 kyrs, in contrast with the wind model.However, there are problems. Although the jet termina-tion regions are expected to be sites of efficient particleacceleration in the SNR and jet scenarios, prominent X-ray emissions are not detected in the terminate regionsof the wings. Moreover, there is no physical explana-tion of where X-ray hotspots form in the intermediateregion of SS433 jets. Goodall et al. (2011a) pointed outdifferences between simulation results and observations.In particular, the connections between the jet and SNRshell are a smooth transition in contrast with simulationresults (see Figure 10 in Goodall et al. 2011a).In terms of the morphology of continuous jets pro-duced by an active galactic nucleus (AGN), some numer-ical simulations of light jets, which have a lower densitythan the surrounding medium, reproduced the structureof the shell and wings (Gaibler et al. 2009; Hodges-Kluck& Reynolds 2011; Rossi et al. 2017). When the jet den-sity is lower than the ISM density, we can use the simpleformula of the light jet to estimate the advance speedof the jet head and the lateral expansion velocity (e.g.,Norman et al. 1982; Krause 2003). The important pa-rameter for the advance speed is the mass ratio of theinjected jet density to ISM density, η = ρ jet /ρ ISM . Theadvance speed of light jets v head can be determined an-alytically from the momentum balance at the contactdiscontinuity between the jets and ISM. The advancespeed is approximately (Krause 2002) v head = √ Aη √ Aη v beam , (1)where A and v beam are respectively the ratio of the beamto head area and the beam velocity. Equation 1 indicatesthat the advance speed also depends on the ratio of thebeam to head. The ratio of the beam to head is closeto unity for a heavy jet ( η ∼
1) and decreases to 0.1 fora light jet ( η (cid:28)
1) (Krause 2002). Thus, light jets de-celerate rapidly, and the propagation enters a non-linearstage. Meanwhile, the lateral expansion of light jets iswell described by the blast-wave equation of motion af-ter the bowshock reachs at ∼ R jet / η . , where R jet isthe jet beam radius. (Krause 2003, 2005). Light jetshave strong backflows whose plasma continuously goestoward the core, and the jets become supersonic expan-sion in the lateral direction. When we simply assumeconstant energy injection E = L kin , jet t and the densityof the homogeneous ISM ρ ISM , the radius of lateral ex- pansion can be written as R lateral = (cid:18) L kin , jet t πρ ISM (cid:19) / = (cid:18) R ηv t (cid:19) / , (2)where L kin , jet = πR ρ jet v is the kinetic energy lu-minosity of the jets. Therefore, v beam , R jet , and η arequantities important to the morphology of light jets.In this paper, we propose new modeling of the mag-netohydrodynamics (MHD) of W50/SS433 generated byjets and backflows. The jets injected during the period20–90 kyrs are two-sided, light, supersonic, and magne-tized. We conduct four simulations for jets with differ-ent radii and density ratios, and we present a mechanismfor the formation of a shell and wings like those of W50.The remainder of the paper is structured as follows. Wedescribe the basic equations and numerical setup of thejets and ISM in Section 2. In section 3, we present ourfour models for MHD simulations. We discuss the re-sults and other formation scenarios of W50/SS433 inSection 4 and give a summary in Section 5. NUMERICAL SETUPWe solve ideal MHD equations in the axisymmetriccylindrical coordinate system (
R, φ, z ). The equationsare ∂ρ∂t + ∇ · ( ρ v ) = 0 , (3) ρ (cid:18) ∂ρ∂t + v · ∇ v (cid:19) = −∇ (cid:18) p + B π (cid:19) + 14 π ( B · ∇ ) B , (4) ∂∂t (cid:18) e + ρv B π (cid:19) + ∇ · (cid:20) ( e + p + ρv v − ( v × B ) × B π (cid:21) = 0 , (5) ∂ B ∂t = ∇ × ( v × B ) , (6) ∇ · B = 0 , (7) ∂f∂t + ∇ · ( f v ) = 0 , (8) where ρ, v , p, B , e, and f are respectively the density,velocity, pressure, magnetic field, internal energy den-sity, and jet tracer function. To distinguish the jetsand ISM, the jet tracer function is set equal to unityfor the injected jet flow and equal to zero for the ISM.We adopt the ideal gas law and give the internal en-ergy density as e = ( γ − p , where γ = 5 / Ohmura et al. magnetic field (Dedner et al. 2002). The number of cellsin our simulation box is ( N R , N z ) = (1200 , < R <
60 pc and − < z <
120 pc. We thus adopt a uniform grid inboth directions, ∆ R = ∆ z = 0 .
05 pc. We apply a sym-metric boundary condition at R = 0 and an outflowboundary condition, which sets a zero gradient acrossthe boundary, at other boundaries.In this work, we assume the distance to SS433 as 5.5kpc (Lockman et al. 2007). The ISM around SS433 thusseems to have a galactic exponential profile (Dehnen &Binney 1998). We adopt a modified exponential densityprofile for the initial density of the ISM: ρ amb = µm p n ( z ≥ − µm p n exp (cid:18) − RR d − z + 20 Z d (cid:19) ( z < − , (9)where R d = 5 . Z d = 30 pc is the scale height of the galactic disc, n = 0 . − , m p is the proton mass, and SS433 islocated at z = 0 and R = 0. Note that the HI intensitymap of the GALFA-HI survey shows that the constantdensity distribution of the ISM is a reasonable assump-tion in the eastern region of SS433 (Peek et al. 2011).The western and eastern directions respectively corre-spond to z < z > K at the location of SS433. Weassume that the uniform magnetic field is parallel tothe z-axis and that the plasma β ≡ p/ ( B / π ) is 10( B z ∼ . µ G).The two-side jets are injected by a cylindrical nozzlealong the z-axis at the origin. We adopt a velocity v jet =0 . × cm s − and initial Mach number M = 10for the jets in all runs. We make calculations for fourmodels, which have different density contrasts η = 10 − and 10 − and different jet radii R jet = 1 and 3 .
33 pc.The jet nozzle length is 2 pc ( | z | < L kin = πR ρ jet v = πR ηρ v (10) ∼ . × η (cid:18) R jet (cid:19) erg s − . We inject a toroidal field B φ = B jet sgn ( z ) sin ( πR/R jet ),where sgn is the sign function, in the jet nozzle. We set β = 1 at R = 0 . R jet and hence B jet = 140 µ G for allruns. The parameters are summarized in Table 1 and2. NUMERICAL RESULTS
Table 1.
Common physical parameters in our runsISM density at the origin ρ . × − g cm − ISM temperature at the origin T KISM magnetic field B z, ISM µ GJet velocity v jet . × cm s − Jet sonic Mach number M jet β β jet Table 2.
Parameters in runs for different jetsRun η (= ρ jet /ρ ) R jet [pc] L kin [erg s − ]H-R3 10 − . × H-R1 10 − . × L-R3 10 − . × L-R1 10 − . × Morphology
The basic structure in all runs has a dynamic behav-ior similar to that found in previous light-jet simulations(e.g., Ohmura et al. 2020). Therefore, we briefly sum-marize the basic structure of light jets and then describethe overall morphology of each run in detail. Figure 1shows counters of the gas density ( top ) and z-directionvelocity ( bottom ) for run L-R1 at t = 88 kyrs. The basicstructures of supersonic jets, such as a beam, backflow,cocoon, and thick shell (shocked ISM) are seen. Sur-rounding plasma mixes into the cocoon owing to Kelvin–Helmholtz (KH) instabilities at the contact discontinu-ity between the jet gas and shocked ISM. In particular,backflows from the tips of the jet interact with eachother, and vortex motions are highly developed in thecocoon. The background density spatially increases at z < −
20 pc in the western direction, and the advanceof the jet in the western direction thus slows apprecia-bly. The magnetic field distribution and its detail arereported in section 3.3.Figure 2 shows the shape of the jet for all runswhen the contact discontinuity between the jet and ISMreaches 100 pc. Because the advance speed of jets de-pends on the injection condition, the times at which thejets reaches 100 pc are different; i.e., H-R3, H-R1, L-R3,and L-R1 are at t = 16 , , , and 88 kyrs respectively.We here define the outer shape of jets as the contact dis-continuity rather than a bow shock. However, it is diffi-cult to identify the contact discontinuity between the jetplasma and shocked ISM because the gases are highlymixed by the KH instability. We therefore define theouter shape as the maximum radius of each z at which Figure 1.
Density ( top ) and z-direction velocity ( bottom ) maps for run L-R1 at t = 88 kyrs. the jet tracer function f is less than 10 − . An analogousdefinition has been widely used in previous works on jetpropagation (e.g., Gaibler et al. 2009). We here denotethe shell radius and the lengths from the origin to theeastern/western edges of the jet as l R , l east , and l west .We find that all runs have a circular shell centeredon the origin when jets reach ∼
100 pc. The ratio of l east to l west is 100 : 75 for all runs. Lighter jets suchas L-R1 and L-R3 form a larger shell. A lower densityratio corresponds to a lower lateral expansion velocity(see equation 2). However, because the kinetic energyluminosity of light jets is less than that of heavy jets,the advance of the light jet becomes much slower thanthat of the heavy jets. Lighter jets thus have enoughtime to form a larger circular shell while they propagatea distance of 100 pc. Run L-R1 therefore reproduces thecharacteristic morphology of W50, including a circularshell and elongated wings. Note that the wings of L-R3are much wider than those of W50. 3.2. Time evolution
In this section, we describe the time evolution of thejet morphology (Figure 3). Top and bottom panels showL-R1 and L-R3, respectively. Colors denote the shapesat t = 8 , , , , , , , ,
88 kyrs for run L-R1and t = 5 , , , , ,
55 kyrs for run L-R3. In theearly period of t (cid:46)
30 kyrs, the shapes of L-R1 and L-R3are spheroids rather than circular shells. The spheroidalshape is well known from previous simulations and issimilar to that of FR II radio sources, such as 3C 452(e.g., Harwood et al. 2015). We find that a circularshell and wings then form. The region of the wing isfilled with back-flowing plasma and it is thus difficultfor the region to expand rapidly like the shell. Once theshell and wings have formed, they expand as self-similarevolution.We next consider what physical parameters determinethe advance speed of the jet. We easily estimate the ad-
Ohmura et al.
Figure 2.
Shapes of jets for all runs. Runs H-R3, H-R1, L-R3, and L-R1 are at t = 16 , , , and 88 kyrs, respectively. Wedenote 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 . vance speed of the jet in the positive z direction adoptingequation 1. Meanwhile, the density distribution in thewestern direction exponentially increases, and the gra-dient effect should thus be included in the expression forthe advance speed: v head = dzdt = (cid:112) η ( z )1 + (cid:112) η ( z ) v beam , η ( z ) = η exp (cid:18) z − Z d (cid:19) . (11)Here, we simply assume that A = 1 and ρ , v beam areconstant values.Figure 4 ( left ) shows the positions of the tips ofthe eastern/western jets as a function of the time foreach run ( solid lines ) and analytic lines ( dashed lines η = 10 − , dotted lines η = 10 − ) of equation 1 ( z > z <
0, towards to the galactic plane).All of the lines at z > A is decreasingin equation 1 because vortices grow downstream of theterminal shock. The jets, which have a large radius, havea small deceleration rate for both runs H and L because A remains at unity for a longer time than in the case ofjets having a small radius. When a jet propagates to-ward the galactic plane (in the western direction), whichmeans the ambient density increases as the distance in-creases ( z < right ) showsthe shell radius l R as a function of time for each run( solid lines ). Dashed lines are the time evolution of theanalytical model (equation 2) whose kinetic energy lu-minosity is L kin , jet = 10 , , and 10 erg s − . Allruns have l R ∝ t / in accordance with the analytical Figure 3.
Time evolution of the shape of jets of runs L-R1 ( top ) and L-R3 ( bottom ). Colors denote the shape at t = 8 , , , , , , , ,
88 kyrs for run L-R1 and t =5 , , , , ,
55 kyrs for run L-R3. solution. Note that equation 2 describes the length ofthe bow shock along the R axis instead of the contactdiscontinuity. l R is defined by the position of the con-tact discontinuity, and the expansion rate in all runs isthus lower than that for equation 2. Note that the ap-proximation for equation 2 is not valid when the densityratio of the ISM to jets is close to unity. Furthermore,jets with larger radii have more kinetic energy luminos-ity, and the speed of radial expansion thus depends onthe jet’s radius.To discuss the formation mechanism of shells andcocoon, we focus on the length ratio of the two-sidejets and the length ratio of the jets and shell radius.In this work, the shell radius of W50 and the lengthsfrom SS433 to the tips of eastern and western wingsare respectively R W50 ∼
48 pc, L east ∼ . L west ∼ . 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 forall runs. The dashed and dotted lines are analytic solu-tions obtained when η = 10 − and 10 − . The numericalresults for all runs and analytic solutions show a linearevolution over time. In this work, we perform calcula-tions until the tip of one side jet reaches 100 pc, but thelength ratio, l east /l west , increases with simulation time.It is thus expected that the length ratio becomes 1.4when l east = 120 pc, which is consistent with radio ob-servations. Interestingly, the length ratios are roughlythe same for all runs when the positions of jets are fixed.These results indicate that l east /l west is determined onlyby the density profile of the ISM. We next investigate the time evolution of the ratioof length from the center to the tip of the jet propaga-tion in the eastern direction to the shell radius l east /l R (right panel of Figure 5). For all runs, l east /l R saturatesat certain values. The lower-density runs of L-R1 andL-R3 have a constant value of l east /l R for a long time.These results mean that outer shapes seem to undergoself-similar expansion in the later phase (see Figure 3).The ratio of l east and l R strongly depends on η . Al-though the advance speed of the jet head is proportionalto η , the speed of radial expansion does not depend on η . Although the advance speed is strongly dependenton η , the speed of radial expansion is weakly dependenton η . Therefore, l east /l R is greater when η is closer tounity. The ratio of the eastern wing to the shell radiusof W50/SS433 is L east /R W50 ∼ .
5, which is consistentwith the L-R1 run. In this work, we assume axisym-metric coordinates. We note that the advance speed inthree-dimensional simulations can be higher than thatexpected from axisymmetric simulations (Perucho et al.2019). 3.3.
Flow and Magnetic Fields
Figure 6 ( top ) presents the distributions of physicalquantities on the jet axis that show the characteristicsof shock waves. The red and blue curves present thepressure ¯ p and velocity gradients ( − ∂ ¯ v z /∂z ) averagedaccording to ¯ f = 1 πR (cid:90) R jet πRf dR. (12)The velocity gradient in the direction of propagationroughly represents the magnitude of the shock wave.The bottom panel of Figure 6 shows stream lines of theaverage velocity. The velocity gradient (blue curve) hasmany sharp spikes while the pressure distribution haswide-tailed peaks. The two are generally correlated andthe peaks are considered shock waves, such as reconfine-ment shocks, oblique shocks, and terminal shocks. Asan example, we see reconfinement shocks at z ∼ − z ∼ −
70, 90 pc) and oblique shocks ( z ∼ −
45, 60 pc)in both side jets. A terminal shock is a location of en-ergy exchange from jet kinetic energy to thermal energyand it produces backflows. When the propagation speeddecreases, the backflow becomes turbulent and has cir-cular and arc-like motion (see the bottom panel of Figure6). The backflow therefore interacts with the beam itself
Ohmura et al.
Figure 4.
Positions of the tips of the eastern/western jets ( left ) and the shell radius ( right ) as a function of time for each run.Dashed and dotted lines in the left panel show analytic results of equation 1 and 11 for runs L and H, respectively. Dashed linesin the right panel show analytic results obtained using equation 2.
Figure 5. l east /l west ( left ) and l east /l R ( right ) as functions of time for each run. Symbol colors denote the position of the tip ofthe eastern jet ( left ) and the shell radius ( right ). Dashed and dotted lines in the left figure show analytic values obtained usingequations 1 and 11 for runs H and L, respectively. and drives the strong oblique shock (Mizuta et al. 2004).We estimate the size of the backflow vortex in section4.1. In our calculation, the oblique shock in the easternjet releases more kinetic energy than the terminal shock.Meanwhile, the western jet shows the opposite correla-tion. The background ISM density distribution increasesin the western direction, and the conversion efficiency ofthe jet kinetic energy to thermal energy at the terminalshock in the western direction is therefore higher thanthat in the eastern direction. The red curve shows thatthe pressure maximum of the jet head at z = −
75 pc is3 times that of the jet in the eastern direction ( z = 100pc). Figure 7 ( left ) shows the distribution of toroidal mag-netic fields and magnetic energy for run L-R1 at t = 88kyrs ( left ) and the distribution of current density, j = ∇ × B ( right ). The magnetic energy is increased byshock compression; in particular, its energy is stronglyenhanced downstream of the terminal shock of the west-ern jet. We find three areas in the cocoon in Figure 7( left-top ): an area of a positive field (blue), an area of anegative field (red), and an area of a weak magnetic field(white). Both backflows, which are antiparallel toroidalfields, interact in the cocoon, and the gas is mixed byvortices. The magnetic fields are therefore dissipatedby magnetic reconnection in this mixing region. Fig-ure 7 ( right ) shows that the interaction of both back- Figure 6. ( top ) Cut along the z-axis of beam-averaged val-ues of − ∂ ¯ v z /∂z (blue) and ¯ p/p inj (red) for run L-R1 at t = 88kyrs. ( bottom ) Velocity stream lines and the outer shape areillustrated in purple and black for run L-R1 at t = 88 kyrs. flows generates current sheets, where particles can beefficiently accelerated by magnetic reconnection. Ourresults would be affected by an artefact of the axisym-metric condition. If we conduct three-dimensional sim-ulations, turbulence would be generated in the cocoon,which would readily convert toroidal and poloidal fields. DISCUSSION4.1.
Dynamics of Backflow
We here make a simple estimation of the interval ofthe oblique shocks. Note that the oblique shocks areformed by the interaction between the beam and back-flow. The interval thus corresponds to the size of thebackflow vortex l bf . For the eastern jet of L-R1, weassume that the radius r bf is about 10 pc (see Figure3). To verify this assumption, we make the simple ar-gument of mass flux conservation. The simulation re-sult gives a z -direction velocity and density of the back-flow of about 0 . v jet and 10 − ρ jet , respectively. Con-servation of mass flux in a hollow cylindrical symme-try gives πR ρ jet v jet = π ( r − r ) ρ bf v bf . We obtain r bf ∼ R jet = 8 pc, and the assumption is reasonable.We assume that the R -direction velocity is about thespeed of sound c bf . To obtain c bf , we first calculate thespeed of sound at the hotspot c h using the Rankine–Hugoniot condition, c = 5 M + 14 M − M c , (13)where M = 10 is the Mach number of the jets while c jet = M − v jet is the speed of sound of the jets. Note that we ignore the advance speed because the advancespeed is much lower than v jet and c h . We can estimate c bf simply using the adiabatic condition: c = (cid:32) R r (cid:33) γ − c , (14)where we assume that the radius of the hotspot is thesame as the radius of the jet. The R -direction velocityis therefore evaluated as c bf = v jet (cid:118)(cid:117)(cid:117)(cid:116)(cid:32) R r (cid:33) γ − M + 14 M − M ∼ . v jet . (15)This value is close to the z -direction velocity of backflow v bf , which is about 0 . v jet (see Figure 1). The oscillat-ing time of the R -direction backflow is τ = r bf /c bf , andthe interval of oblique shocks is thus evaluated as l bf = 2 τ v bf ∼
25 pc . (16)This value is roughly consistent with our numerical re-sult. The bottom panel of Figure 6 shows that the back-flow departs from z = 90 pc and arrives at z = 65pc, where a strong oblique shock is excited. However,the actual dynamics of the beam and backflow are morecomplex and many weak shocks are thus excited alongthe beam. 4.2. Synchrotron Emission
We create a pseudo synchrotron image of W50/SS433using physical values from MHD simulations. To calcu-late the synchrotron emissivity, we convert the two-dimensional axisymmetric cylindrical coordinates tothree-dimensional Cartesian coordinates and assumea viewing angle of 78 . ◦ . We calculate and integrate theemissivity at each cell along the line-of-sight direction: I = (cid:90) LoS J sync ds, (17)where I and J sync are respectively the synchrotron in-tensity and emissivity. Our simulations do not model theevolution of non-thermal electrons, and hence some as-sumptions are needed to calculate the synchrotron emis-sivity. One is that the spectral index is uniform every-where. The other is that the number density of non-thermal electrons is proportional to the gas density andthermal energy density. Under these assumptions, thesynchrotron emissivity is given by (Jun & Norman 1996) J sync = C ( ν ) ρ − α p α B α +1 ⊥ , (18)where C , B ⊥ , ν , and α are a normalized constant, themagnetic field component in the plane of the sky, the0 Ohmura et al.
Figure 7. ( left ) Toroidal component of magnetic field ( top ) and magnetic energy ( Bottom ) maps of run L-R1 at t = 88 kyrs.( right) Current density maps of run L-R1 at t = 88 kyrs. The plot area corresponds to the green rectangle in the left panel.The two backflows, which are antiparallel toroidal fields, interact in the cocoon, which generates current sheets. Figure 8. ( left ) Synchrotron intensity map overlaid with the outer shape (white) of the jet in run L-R1 on a linear scale. ( right )Map of the distribution of jet particles in run L-R1 on a linear scale. Both maps are plotted for a viewing angle of 78 . ◦ frequency of radiation, and the spectrum index, respec-tively. We adopt that α = 0 .
52 because some radioobservations provide the flatter spectrum at the easternwing of W50 (Dubner et al. 1998; Gao et al. 2011). Thecalculation of polarized emissions is beyond the scope ofthis paper due to axisymmetry.Figure 8 ( left ) shows the synchrotron intensity mapat t = 88 kyrs for run L-R1. The synchrotron powerstrongly depends on magnetic energy, and a weak emis-sion is thus observed along with the shell (see Figure7). The bright region in the western wing is around theoblique shocks and the termination shock. Bright small spots on the jet axis are formed by oblique shock whosespacing is defined by equation 16.We can see the wide bright area both in the west andeast wings which consist of thin filaments. The accu-mulated magnetic field along the contact discontinuityis the origin of these bright regions. We consider thatthese bright spots correspond to SS433 X-ray hot spots.For example, the magnetic filament around (R, z) = (10,80) is the origin of the bright filament closest to the east-ern end. On the other hand, because the magnetic fieldaround the terminal shock is diffused, there is no X-rayradiation around that region.1The synchrotron intensity map obtained from numer-ical simulation is similar to X-ray observations of W50(Brinkmann et al. 1996a,b, 2007; Safi-Harb & ¨Ogelman1997). In particular, the positions of X-ray hotspots ofW50 (e1, e2, w1, and w2) are consistent with those in thenumerical results. The oblique shocks are thus thoughtto be rather efficient particle accelerators in contrastwith the terminal shocks in the long-term evolution ofjets. Sudoh et al. (2020) constructed an analytic modelfor non-thermal leptonic emission from knots in SS433jets, and their results suggest that the X-ray observationcan be reproduced by a synchrotron emission when theacceleration process is efficient.We next discuss the distribution of electrons emit-ting radio waves in W50/SS433. One main differencebetween microquasar jets and AGN jets is the out-burst age; i.e., the age of relativistic electrons in thejet cocoon. The synchrotron cooling time t c of non-thermal electrons in the microquasar jets is t c ∼ . × B − E − s, where E e is the energy of electrons. Themagnetic field strength of the eastern and western wingsis about 1 − − µ G (Safi-Harb & ¨Ogelman 1997) andthe cooling time scale is thus estimated as t c ∼
100 kyrswhen the frequency is 1 GHz and the magnetic field is100 µ G. The plasma created around the terminal oroblique shock is thus not cooled to a few GeV and emitsradiation at a centimeter wavelength. We calculate theshocked particle content in W50/SS433, F , to integratethe tracer function, f : F = (cid:90) LoS f ds. (19)Figure 8 ( right ) is the same as Figure 8 ( left ) but wedisplay a map of the distribution of jet particles . Theplasma accumulates in the cocoon after passing the ter-minal shock and oblique shocks. The lifetime of elec-trons for synchrotron radiation is much longer than thejet active time (simulation time), and relativistic elec-trons in the shell can thus emit continually. The in-tensity is lower in both wings than in the shell becauserelativistic electrons are advected by backflow so thatthere are fewer electrons in both wings. This character-istic is also found for the radio emission of the easternwing. Another radio characteristic of the eastern wingis an extended radio filament in the terminal region. Inaddition, Brinkmann et al. (2007) found an X-ray ring inthis region, which corresponds spatially to the terminalshock of the SS433 jet. However, our numerical resultindicates that the shock size is only a few parsecs. Thefilament is thus formed by the backflow of the termi-nal shock. Additionally, Sakemi et al. (2018a) reportedthe existence of helical fields that coil the eastern wings.
Figure 9.
Radial distance for SNRs ( dashed ) and jets ( solid )as a function of time in the case of ISM gas densities ρ ISM =0 . m p (blue) , . m p (green) , and 1 m p (red) g cm − . Weassume that the kinetic energy luminosity of jets and theinitial explosion energy of SNRs are L k , = 10 erg s − and E SNR , = 10 erg, respectively. This result is consistent with our backflow scenario. Wenext discuss the shell emissions of W50. We find thatmagnetic reconnection generates current sheets in theshell (see Figure 7) and re-accelerated electrons in cur-rent sheets provide a prominent radio emission. This ra-dio emission corresponds to the filament structure in thesouthern region of W50; otherwise, the MeV γ -ray emis-sion observed in the northern region of W50 is explainedby this acceleration scenario. In this simulation, we as-sume an axisymmetric condition and the magnetic fieldin the shell of W50 is thus exhausted by magnetic re-connection. However, if we carry out three-dimensionalsimulations, the magnetic field can be stored in the shellbecause the flow can escape easily in the 4 π direction.The synchrotron emission from the shell would thereforebe enhanced. 4.3. SNR Evolution
Previous studies have suggested that the SNR pro-duced by the compact object SS433 formed the shellstructure of W50, and we thus discuss the differencein radial expansion between our backflow scenario andSNR scenarios. We next compare the features of ra-dio emissions of W50/SS433 with those of other SNRs.We mention again the shell diameter is D = 96 pc.The typical evolution of an SNR is described in theSedov phase (i.e., the adiabatic phase). An expandingSNR is then prescribed by the snowplow phase when ra-diative losses become important. We assume that theshell of W50 is in the Sedov phase and radiative cool-ing is thus negligible. Note that the radial expansion2 Ohmura et al. of the Sedov phase is faster than that of the snowplowphase. The shell radius of the SNR develops in the Se-dov phase as R ∝ ( E SNR t /ρ ISM ) / , where E SNR is theinitial explosion energy of the SNR. The dashed linesin Figure 9 show the radius of the SNR as a functionof time for three different values of the ISM density.We simply assume that the ISM has a constant den-sity ρ ISM = 1 m p , . m p , and 0 . m p g cm − and E SNR , = 10 erg. When the ISM density is greaterthan 0 . m p , it takes more than 100 kyrs to reach a shellradius of 48 pc. The ISM is thus about 0 . m p , which isa very low density with which to form the W50/SS433shell. However, we note that W50/SS433 is locatednear the galactic plane, and HI emissions are observedaround W50 (Peek et al. 2011, Sakemi et al. submittedto PASJ). Such a low-density environment is thus un-likely. In another case, we believe that the W50 nebulaformed in a high-power explosion ( E SNR (cid:29) erg),such as a hypernova or multiple supernovas. In con-trast, we plot the time evolution of the radial expan-sion for the light jets for jet kinetic energy luminosity L kin , = 10 erg s − in Figure 9. The approximatesolution of the radial expansion of jets is R jet ∝ t / .We mention again the Sedov solution is proportionalto t / . The radius obtained using the jet model be-comes larger than that obtained using the SNR modelwhen the total energy exceeds the SN explosion energy, E SNR = 10 erg. The time is about t ∼
40 kyrs forthese model settings. In the case of jets instead of anSNR, the shell of W50 reaches 48 pc within 100 kyrswhen the ISM is 0 . m p g cm − . Therefore, jets havinglong-term activity, for which a lifetime of >
10 kyrs isexpected, are a good candidate for the formation mech-anism of W50/SS433.We next discuss radio observations of typical SNRsand W50/SS433. Galactic SNRs have a radio obser-vational relationship between the radio surface bright-ness and the diameter of the SNR, the so-called ra-dio 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 andgreen diamonds) and W50/SS433 (red star). The figureshows that W50/SS433 is prominent and has a large sizecompared with galactic SNRs. We note the observationsof four prominent galactic SNRs (CTB 37A, Kes97, CTB37B, and G65.1+0.6) interacting with molecular clouds,which enhances radio emissions. The wings of W50 havea radio surface brightness comparable to that of the shellof W50 (Dubner et al. 1998). Previous hydrodynamicalsimulations clearly showed that the jetted gas and thegas-related supernova explosions are separated; in par-ticular, the jetted gas is not distributed around the shell Figure 10.
Surface brightness vs. diameter at 1 GHz forshell SNRs (blue circles and green diamonds) in Pavlovicet al. (2014) and W50 (red star). The solid line is the orthog-onal regression in Pavlovic et al. (2014). Green diamonds arefour prominent galactic SNRs (CTB 37A, Kes97, CTB 37B,and G65.1+0.6), which interact with molecular clouds. region (e.g., Goodall et al. 2011a). Although some ob-servations suggest that there are molecular clouds asso-ciated with jets in the western wing of W50 (Yamamotoet al. 2008; Liu et al. 2020), there is no critical evidenceof interacting molecular clouds, which would acceleratenon-thermal particles. Additionally, the jetcloud inter-action is not strong and it thus does not enhance radioemissions throughout the western wing. Therefore, theorigin of the prominent emission from the shell of W50is not the same as that of other prominent SNRs.Finally, we mention that the shell of W50 maintainsa nearly perfectly circle-like structure. Stafford et al.(2019) investigated the relationship between the radiomorphology and the sizes (ages) of SNRs. They foundthat larger (older) SNRs are more elliptical and/or elon-gated because the inhomogeneous ISM, turbulent struc-ture, and molecular clouds affect their morphology. Aswe mentioned above, W50 is in a large size (old) cate-gory if the shell of W50 comprises SNRs.4.4.
Comparison with other W50/SS433 formationscenarios
In this section, we compare our results with those ofGoodall et al. (2011a). We first mention that the radialexpansion of the SNR is too fast in the simulations con-ducted by Goodall et al. (2011a) compared with the Se-dov solutions and hydrodynamical simulations of SNRs.They assumed the kinetic energy of SNRs is ∼ ergand the ISM density is ∼ . − , and the shell ra-dius of SNRs reaches 50 pc within 20 kyrs. They mighthave underestimated the normalization of the kinetic en-3ergy of SNRs or overestimated the density unit of theISM density by an order of magnitude. We also mentionthe advance speed of jets in Goodall et al. (2011a). Asan example, model Jet1 in Goodall et al. (2011a) hasa young age of 2.2 kyrs when the jet crosses the tip ofthe eastern wing. The average advance speed of modelJet1 almost maintains the launch velocity of 0 . c fromSS433; i.e., 120 pc / . ∼ . c . Thus, the ad-vance of model Jet1 does not decelerate. Note that theresults of kinematical studies of the radio observationsof W50 indicate that an upper limit to the proper mo-tion of the radio filaments of the eastern wing is 0 . c (Goodall et al. 2011b). Jets that do not decelerate mustbe comparable or greater than the ISM in density. Someprevious numerical simulations have modeled protostarjets for heavy jets (e.g., Ustamujic et al. 2016). How-ever, numerical and observational studies of AGN jetsand microquasar jets suggest that extensive cocoons de-velop when the surrounding medium is lighter than themedium of the jets.Panferov (2017) constructed a model combined withthe theory of the dynamics of SNRs, wind bubbles, andjets for the formation of W50. Their jet model was basedon turbulent jets in Landau & Lifshitz (1987). In thecase of a wind bubble, if the binary system of SS433 isa super-Eddington accretion, the age of the wind bub-ble needs to exceed 300 kyrs to produce the large shellof W50. They also modeled the dynamics of SNRs inthe Sedov phase and snowplow phase, describing thedynamics using the fitting function of Kim & Ostriker(2015). They argued that the age of W50 as an SNR isat least 97 kyrs, which agrees with our results in Figure9. The main result of Panferov (2017) is that the shellof W50 formed in a supernova explosion about 100 kyrsago, and the SS433 jets produce both wings and pushedup the pre-existing SNR within 27 kyrs. However, hy-drodynamical studies showed that interaction betweenthe jets and SNR does not push up the shell of W50in the radial direction (Goodall et al. 2011a) Therefore,some other mechanism is needed for the model of Pan-ferov (2017) to exchange energy efficiently between thejets and SNR shell. SUMMARYWe conducted a series of axisymmetric MHD simula-tions of two-side jets to study the parameter dependencefor the formation of shell and wings of the W50/SS433system. We argue that light and supersonic jets canproduce the shell and wing morphology in a backflowscenario. The main results of this work are as follows. • At an early time ( t (cid:46)
30 kyrs), the morphology oflight jets ( η < − ) had a spheroidal shape. Af- terward, the shell and wings formed and they ex-panded under self-similar evolution. In contrast,slightly heavy jets ( η ∼ − ) cannot form a largeshell because the advance is too fast comparedwith the radial expansion. • Jets having a large radius do not decelerate likejets having a smaller radius. Moreover, larger jetsprovide much kinetic energy luminosity, and theradial expansion of the cocoon thus becomes fast.We can therefore produce wings and a shell of ar-bitrary size by appropriately selecting parameters,such as the jet radius and density contrast. • The exponential ISM profile produces the asym-metric wings observed for W50/SS433 in previousworks (Zavala et al. 2008; Goodall et al. 2011a;Millas et al. 2019). We find that the length ratiobetween the eastern and western wings, l east /l west are roughly the same for all runs when the posi-tions of jets are fixed. Thus, l east /l west determinedonly by the density profile of ISM. Furthermore,our runs using the ISM profile of our model pro-vide an asymmetric shape that is consistent withobservations. • The main challenge to the jet scenario is the re-quirement for the long-period activity of the jets.Run L-R1 requires a long activity lasting morethan 100 kyrs to result in the observed size of W50.The advance speed and the speed of radial expan-sion respectively depend on the density contrastand density of the ISM. Therefore, if the ISM isslightly lower than 0 . − , our scenario over-comes the issues noted above. • A large amount of jet kinetic energy is convertedinto thermal energy for jets propagating in thegalactic plane because these jets interact with thehigh-density ISM. The interaction between thebackflows and beams drive strong oblique shocksat the intermediate point of beams. The positionsof oblique shocks correspond to X-ray hotspotsof W50. Meanwhile, highly turbulent flows drivemagnetic reconnection, which generates currentsheets in the shell region where the two backflowshave interacted. • We discussed the jet scenario and SNR scenario.For an SNR interpretation, an extremely-low-density atmosphere, n ISM (cid:28) . − , or a high-energy explosion, E SNR (cid:29) erg, is requiredfor the SNR’s froward shock to reach about 50 pc.Moreover, the shell of W50 is larger and much4 Ohmura et al. more prominent than that of other galactic SNRs.Although the jet scenario strongly depends on theactivity time and the total amount of energy, theshell of W50 could have been formed by back-flow within a period of 100 kyrs. Notice that ournumerical results can explain the morphology ofW50/SS433, but they do not entirely rule out theSNR scenario.Future work should treat some of the key physicsthat we ignored. First, three-dimensional effects affectthe dynamics of jets. The advance speeds in three-dimensional simulations are higher than those in ax-isymmetric simulations, and turbulence is generated bynon-axisymmetric motions. Second, we ignored windsfrom the companion star and/or accretion disk. The in-teraction between winds and jets might positively affectthe formation of the shell of W50. If SS433 is an ultra-luminous X-ray source, then the state of the accretiondisk is a supercritical accreting flow. Radiation hy-drodynamics simulations of supercritical accretion flowsuggest that strong winds ( ∼ erg s − ) can be drivenby radiation pressure (Kawashima et al. 2009). From an observational viewpoint, some ultraluminous X-raysources have bubbles on the scale of a few tens of par-secs, similar to the case for W50 (e.g., Pakull et al. 2010;Berghea et al. 2020). Finally, we need to implement theprecession of jets. However, we mention again that theprecessions do not affect the large-scale morphology ofW50 because the jets are well collimated within 0.1 pcfrom SS433 and they become hollow (Millas et al. 2019).We thank the anonymous referee for comments thathave improved the clarity of this paper. We are gratefulto Dr. H., Yamamoto for useful discussion. This workwas supported by JSPS KAKENHI Grant NumbersJP20J12591 (T.O.), JP19K03916, JP20H01941 (M.M.)and JP20J13339 (H.S.). Our numerical computationswere carried out on the Cray XC50 at the Centerfor Computational Astrophysics of the National Astro-nomical Observatory of Japan. The computation wascarried out using the computer resource by ResearchInstitute for Information Technology, Kyushu Univer-sity. We thank Glenn Pennycook, MSc, from EdanzGroup (https://en-author-services.edanzgroup.com/ac)for editing a draft of this manuscript.REFERENCES Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2018,Nature, 564, E38, doi: 10.1038/s41586-018-0688-8Asahina, Y., Ogawa, T., Kawashima, T., et al. 2014, ApJ,789, 79, doi: 10.1088/0004-637X/789/1/79Begelman, M. C., Sarazin, C. L., Hatchett, S. P., McKee,C. F., & Arons, J. 1980, ApJ, 238, 722,doi: 10.1086/158029Berghea, C. T., Johnson, M. C., Secrest, N. J., et al. 2020,ApJ, 896, 117, doi: 10.3847/1538-4357/ab9108Blundell, K. M., & Bowler, M. G. 2004, ApJL, 616, L159,doi: 10.1086/426542Blundell, K. M., Laing, R., Lee, S., & Richards, A. 2018,ApJL, 867, L25, doi: 10.3847/2041-8213/aae890Bordas, P., Bosch-Ramon, V., Paredes, J. M., & Perucho,M. 2009, A&A, 497, 325,doi: 10.1051/0004-6361/200810781Brinkmann, W., Aschenbach, B., & Kawai, N. 1996a, A&A,312, 306—. 1996b, A&A, 312, 306Brinkmann, W., Kotani, T., & Kawai, N. 2005, A&A, 431,575, doi: 10.1051/0004-6361:20041768Brinkmann, W., Pratt, G. W., Rohr, S., Kawai, N., &Burwitz, V. 2007, A&A, 463, 611,doi: 10.1051/0004-6361:20065570 Broderick, J. W., Fender, R. P., Miller-Jones, J. C. A., et al.2018, MNRAS, 475, 5360, doi: 10.1093/mnras/sty081Cooper, A. J., Gaggero, D., Markoff, S., & Zhang, S. 2020,MNRAS, 493, 3212, doi: 10.1093/mnras/staa373Dedner, A., Kemm, F., Kr¨oner, D., et al. 2002, Journal ofComputational Physics, 175, 645,doi: 10.1006/jcph.2001.6961Dehnen, W., & Binney, J. 1998, MNRAS, 294, 429,doi: 10.1046/j.1365-8711.1998.01282.xDubner, G. M., Holdaway, M., Goss, W. M., & Mirabel,I. F. 1998, AJ, 116, 1842, doi: 10.1086/300537Farnes, J. S., Gaensler, B. M., Purcell, C., et al. 2017,MNRAS, 467, 4777, doi: 10.1093/mnras/stx338Gaensler, B. M., Green, A. J., & Manchester, R. N. 1998,MNRAS, 299, 812, doi: 10.1046/j.1365-8711.1998.01814.xGaibler, V., Krause, M., & Camenzind, M. 2009, MNRAS,400, 1785, doi: 10.1111/j.1365-2966.2009.15625.xGao, X. Y., Han, J. L., Reich, W., et al. 2011, A&A, 529,A159, doi: 10.1051/0004-6361/201016311Goodall, P. T., Alouani-Bibi, F., & Blundell, K. M. 2011a,MNRAS, 414, 2838,doi: 10.1111/j.1365-2966.2011.18388.xGoodall, P. T., Blundell, K. M., & Bell Burnell, S. J.2011b, MNRAS, 414, 2828,doi: 10.1111/j.1365-2966.2011.18809.x5