Multiple Jets in the Bursting Protostar HOPS 373SW

We present the outflows detected in HOPS 373SW, a protostar undergoing a modest 30% brightness increase at 850 μm. Atacama Large Millimeter/submillimeter Array observations of shock tracers, including SiO 8–7, CH3OH 7k–6k, and 12CO 3–2 emission, reveal several outflow features around HOPS 373SW. The knots in the extremely high-velocity SiO emission reveal the wiggle of the jet, for which a simple model derives a 37° inclination angle of the jet to the plane of the sky, a jet velocity of 90 km s−1, and a period of 50 yr. The slow SiO and CH3OH emission traces U-shaped bow shocks surrounding the two CO outflows. One outflow is associated with the high-velocity jets, while the other is observed to be close to the plane of the sky. The misaligned outflows imply that previous episodic accretion events have either reoriented HOPS 373SW or that it is an unresolved protostellar binary system with misaligned outflows.


INTRODUCTION
Variable and episodic accretion is recognized as a standard feature of the star formation process (Fischer et al. 2022).This expectation is supported by observations across a broad range of diagnostic signatures, including FU Orionis and EX Lup type brightening events (e.g., Herseokholee@kasi.re.kr lee.jeongeun@snu.ac.kr big 1966; Hartmann & Kenyon 1996;Audard et al. 2014), water and CO snow lines that extend beyond expectations for present-day luminosities (Lee 2007;Jørgensen et al. 2015;Hsieh et al. 2019), and multiple clumps within individual jets (Lee et al. 2010) and outflows (Plunkett et al. 2015).Furthermore, near-infrared, midinfrared, and sub-mm monitoring observations provide direct evidence of variable accretion as a common protostar phenomenon (Rice et al. 2015;Park et al. 2021).
Protostellar jets and outflows are indirect records of the variable accretion history, being natural byproducts of the protostellar accretion process (e.g., Frank et al. 2014;Lee et al. 2020;Fischer et al. 2022).Furthermore, the mass ejection rates are seen to be proportional to the mass accretion rates (e.g., Lee et al. 2020).Within the outflows, high-velocity jets show quasi-episodic knots from which we can derive the timescales of the variability (Plunkett et al. 2015;Lee et al. 2020;Jhan et al. 2022).
HOPS 373 is a protostellar binary system located in the Orion B cloud (e.g.Johnstone et al. 2001) at a distance of 428 pc (Kounkel et al. 2018) showing a modest 30% brightness variability at 850 µm in the James Clerk Maxwell Telescope (JCMT) Transient Monitoring Survey (Yoon et al. 2022).HOPS 373 consists of two continuum sources with a separation of 3. ′′ 6 (1500 au at the distance of 428 pc): referred to as the northeast (NE) and southwest (SW) sources (Tobin et al. 2015(Tobin et al. , 2020;;Yoon et al. 2022).Based on the spectral energy distribution, the estimated bolometric luminosity and the bolometric temperature of the HOPS 373 system are 5.3 L ⊙ and 37 K, respectively (Kang et al. 2015;Furlan et al. 2016).
Complex outflow features are observed in HOPS 373.The 12 CO 2-1 emission shows both a large-scale southeast (SE)-northwest (NW) outflow and a small-scale east-west outflow (Tobin et al. 2016;Yoon et al. 2022).The dynamical age and the inclination of the largescale CO outflow are 10 4 years and 40 • , respectively (Furlan et al. 2016;Nagy et al. 2020).The small-scale CO outflow is associated only with HOPS 373SW.In contrast, no small-scale molecular emission has been detected in HOPS 373NE, except for the C 17 O 3-2 line (Yoon et al. 2022;Lee et al. 2023a).Yoon et al. (2022) argued that the large-scale SE-NW outflow either is associated with the NE source or that wiggles within the outflow from the SW source changed its outflow direction from a blueshifted large-scale SE outflow to a red-shifted small-scale east outflow.
High spatial resolution ALMA observations show prominent jets and bow shocks ejected from HOPS 373SW.In this work, we investigate the properties of these observed jets and outflows and derive the possible link with the observed variable accretion.The observations and archive data sets are described in Sections 2. In Sections 3-6, the properties of SiO and CH 3 OH jets/outflows are investigated.We discuss the stellar mass, episodic accretion, and the origin of the two outflows in Section 7 before summarizing our results in Section 8.

OBSERVATIONS
We investigated the jets and outflows in HOPS 373SW using Cycle 7 ALMA observations (2019.1.00386.T, PI: Jeong-Eun Lee) and Cycle 4 ALMA archive data (2015.1.00041.S, PI: John J. Tobin).Detailed observations and data reduction can be found in Lee et al. (2023a) for the Cycle 7 data and in previous works (Tobin et al. 2020;Yoon et al. 2022) for the Cycle 4 data.
The Cycle 7 observations were carried out on 2021 August 2 and 6 with thirty-eight 12-m antennas, and a longest baseline of 5.6 km.The total integration time on the source was 61 minutes.The jet tracers SiO 8-7 (347.330631GHz) and CH 3 OH 7 k -6 k (338.34-338.73GHz) (Tafalla et al. 2010;Tychoniec et al. 2021) were covered with this observation.The SiO and CH 3 OH lines were cleaned using the natural weighting with the uvtaper of 0. ′′ 2 (∼455 kλ) to increase the sensitivity, resulting in the beam size of ∼0.′′ 23 with a rms of ∼3 mJy beam −1 over the 1 km s −1 spectral resolution.In addition, two CH 3 OH lines at 338.408698 GHz and 338.344588GHz were averaged to increase the signal-to-noise ratio (S/N).A high resolution (∼0.′′ 08) SiO image with a rms of ∼2 mJy beam −1 over the 1 km s −1 was also cleaned using natural weighting without uvtaper in order to study the jets near the continuum peak.
12 CO 3-2 (345.7959899GHz) was covered by the Cycle 4 archive data, which were observed on 2016 September 3 and 4 and on 2017 July 19.The 12 CO line was cleaned using the natural weighting with the uvtaper of 500 kλ, resulting in a beam size of 0. ′′ 25 with a noise level of 20 mJy beam −1 over 1 km s −1 spectral resolution.
The Atacama Compact Array (ACA) archive data toward HOPS 373 were also used to check the extended features.These observations were carried out in March 2019 as a part of the program 2018.1.01565.S (PI: Tom Megeath).Lines of C 18 O 2-1 and N 2 D + 3-2 were covered with a beam size of 7. ′′ 4 × 4. ′′ 9 and a spectral resolution of 0.34 and 0.16 km s −1 , respectively.

MORPHOLOGY OF JETS AND OUTFLOWS
HOPS 373SW shows prominent jets and outflows.Figure 1 presents the SiO channel maps over the peak intensity map of CH 3 OH.In this work, we adopt the systemic velocity (V sys ) of 10.13 km s −1 derived from the 13 CH 3 OH lines (Lee et al. 2023a).All velocities are presented using the relative velocity with respect to this systemic velocity (∆V = V − V sys ).Typically for outflows, the CH 3 OH emission favors a low-velocity gas, while the SiO emission is more abundant in the fastest gas (Tafalla et al. 2010).For HOPS 373SW, the CH 3 OH emission has only slow components (|∆V | <4 km s −1 ) as shown in Figure 2. In contrast, the SiO emission has multiple velocity components: slow (|∆V | < 15 km s −1 ), fast (15 km s −1 < |∆V | < 40 km s −1 ), and extremely high velocity (EHV, |∆V | > 40 km s −1 ).Note that the observation window covers only up to −70 km s −1 from the systemic velocity to the blue, while the red-shifted SiO emission is detected up to ∼90 km s −1 .Thus, we cannot rule out that the high-est velocity component in the blue-shifted SiO emission has been missed.Previously, in low mass protostars, the SiO emission has been observed in clumps (Tychoniec et al. 2021) rather than in the continuous outflow denoted by the CO, and these clumps have different peak velocities from the underlying outflow along with narrow line widths (<15 km s −1 ).Due to this situation, the distribution of outflow shocks is more clearly revealed by the peak intensity (moment 8 in the CASA task 'immoments', the maximum value of the spectrum in each pixel) maps as shown in Figures 2 and 3.
The EHV and fast SiO emission components trace the collimated NE-SW jet revealed by the EHV CO emission (see Figure 8 in Yoon et al. 2022, and the bottom panel of Figure 4).Figure 5 shows that the NE-SW SiO emission is also associated with Spizter/IRAC 4.5 µm intensity tracing the H 2 emission in the NIR (green contours in Figures 1 and 5, Tobin et al. 2007;Velusamy et al. 2011).The H 2 emission is detected only in the blue lobe (Yoon et al. 2022).The red-shifted (NE) infrared emission might be extincted by the optically thick protostellar envelope and, therefore much fainter than the blue-shifted (SW) emission.The EHV and fast SiO emission consists of several knots in the red (HR1, HR2,...) and blue (HB1, HB2, and HB5) lobes, as shown in Figure 1.The EHV SiO knots in the red lobe are well developed up to a terminal knot, HR5.On the other hand, the EHV SiO knot HB1 in the blue lobe is found only near the protostar, and the fast SiO knot HB2 and HB5 in the middle panel of Figure 1 appear as counterparts to the red-shifted EHV SiO knots HR2 and HR5, respectively.The redshifted SiO jet shows a wiggle from north to south as reported in the EHV CO map (Yoon et al. 2022), for which we provide a plausible jet precession solution in Section 7.3.The SiO and CO emission within the blue lobe shows a point  symmetry to those in the red lobe with respect to the continuum peak.
Figure 3 shows the peak intensity and the absolute peak velocity maps for HR1 and HB1 near the continuum peak using the high spatial resolution images (∼0 ′′ .08).They are derived after smoothing with a boxcar of 5 km s −1 .The knots HR1 and HB1 consist of several subknots.In the red lobe, the sub-knot HR1a exhibits a velocity gradient in the southeast direction.HR1b and HR1c appear to be connected, following a curved trajectory.There is an increase in velocity from HR1b to HR1c.HR1d and HR1e are weak compared to HR1a and HR1b.Counterparts in the blue lobe are identified by their distance from the protostar and peak velocity.
The slow SiO and CH 3 OH emission trace the bow shocks and two distinct outflows as shown in the top panel of Figure 1 and 2, respectively.The NE-SW outflow (hereafter Out-flow1) is confined by U-shaped bow shocks, LR1 and LB1 (see the black curves in Figures 1, 2,  and 4), which could be related to the above jet.The SE-NW outflow (hereafter Outflow2) is also surrounded by U-shaped bow shocks, LR2 and LB2.Similar U-shapes of CH 3 OH and SiO emission are seen in L1157 and they trace the wings of the bow shock (Gómez-Ruiz et al. 2013;Codella et al. 2020).
The distribution of SiO and CH 3 OH emission is different between the two HOPS 373SW outflows.In Outflow1, the wings of LR1 are mainly traced by the SiO emission, while those of LB1 are traced by the CH 3 OH emission.In contrast, only SiO emission (e.g., HR5, HR5L, HB5, and HB5L) is detected at the apexes of LR1 and LB1 without corresponding CH 3 OH emission.In Outflow2, the CH 3 OH emission is detected at both the apexes and in the wings.Furthermore, slow SiO emission is clearly detected in LB2 but only marginally in LR2, while the tip of LB2 shows both red-and blue-shifted slow components of SiO and CH 3 OH emissions im- plying that the tip might be close to the plane of the sky.
Figure 4 shows the CO channel maps overlaid on the SiO and CH 3 OH peak intensity (moment 8) maps.The EHV CO emission (bottom panel) traces the jets shown by the EHV SiO emission.The fast (middle panel) and slow (top panel) blue-shifted CO outflow components are confined by the wings of LB1, implying an association with Outflow1.The red-shifted Outflow1 can be traced by the fast red-shifted CO emission and the slow SiO emission (the wings of LR1).
The CO outflow associated with Outflow2 might be resolved out, except for the observed slow SE (red) and NW (blue) CO outflow components (A and B in the top panel of Figure 4).The velocities may therefore be close to the systematic velocity and Outflow2 oriented close to the plane of the sky, as revealed by the CH 3 OH emission at the tip of LB2.Note that we cannot rule out that all of the slowest CO emission in the map, except for the emission around the apexes of LR1 and LB1, might be related to Outflow2 because the slow CO emission in the north (south) wall within ∼3 ′′ of the red (blue) lobe seems to connect to the wings of LR2 (LB2).In this case, the outflow velocity might have been decelerated through interaction with the ambient medium and have slowed to close to the systemic velocity near the apex of LR2 and LB2.

POSITION VELOCITY DIAGRAMS OF JETS AND OUTFLOWS
Figure 6 shows the position velocity diagram along the R.A. (jet) direction within |∆Dec.|<1.′′ 5.The knots in the red-shifted jet have an observed jet velocity of ∼55 km s −1 and a terminal velocity of ∼45 km s −1 at HR5.On the other hand, the knots in the blue-shifted jet have decreasing velocities with ∼-55 km s −1 at HB1, ∼-35 km s −1 at HB2, and a terminal velocity of ∼-20 km s −1 at HB5. Velocity asymmetries in the red and blue lobes are also observed in the SiO jets of L1448C (Yoshida et al. 2021) and L1157 (Podio et al. 2016), as well as in the optical/infrared jets from Class II sources (Hirth et al. 1994).This asymmetric velocity could originate from the launching mechanism itself, for example, differences of launching velocity on each side of the disk by the magnetic field (Dyda et al. 2015) or the accreting mass (Liu & Shang 2012), or as a result of interaction with an asymmetric environment (Podio et al. 2016).In addition, this could originate due to different inclination angles in the blue and red lobes due to the wiggling of jet.However, the distance of terminal knot HB5 is closer to the protostar compared to HR5; thus, the velocity of the blue-shifted jet should be slower than that of the red-shifted one.The HR1 knots could have been ballistically ejected.The yellow lines in the upper insert of Figure 6 indicate velocity gradi-ents along the distance from the origin with 270 km s −1 arcsec −1 and 120 km s −1 arcsec −1 .Therefore the sub-knots HR1a and HR1b would have been ejected 7.5 tan(i) and 17 tan(i) years ago, respectively, where i is the inclination of jet from the plane of sky.The red lines indicate the velocity gradients for ejection events that occurred previously with a similar time interval (9.4 tan(i) years).It is plausible that subknots HR1d and HR1e might have been ballistically ejected with the same cadence.On the other hand, sub-knots in HB1 do not show the velocity gradient along the distance despite their proximity to the velocity gradient associated with each counterpart.The absence of the velocity gradient in the HB1 could be attributed to the limitations in velocity coverage and/or the weaker emission in the blue lobe compared to the red lobe.The insert panels are zoomed in the high-velocity jet near the continuum source in the high angular resolution image.The emission is averaged within |∆Dec.|<0.′′ 5, and the start and step of the contour levels are 3 σ (1σ= 0.4 mJy beam −1 ).Both the SiO spectra are smoothed with a boxcar of 3 km s −1 to increase the S/N.The emission within the green boxes at the continuum peak are associated with complex organic molecules (see Lee et al. 2023a).The yellow and red lines in the insert panels are the slopes for the ballistic ejection (see the Section 4).
Figure 7 shows the position velocity diagram of CO 3-2 (image and white contours), SiO 8-7 (red contours), and CH 3 OH (green contours) along the Declination (perpendicular to the jet axis).Near the protostar at ±1. ′′ 1 (bottom panels), both CO and SiO emission connect the high-velocity components (HR1e and HB1e) to the ambient envelope material, implying that the ambient material is being accelerated by the jet (Rabenanahary et al. 2022) as also shown in HH 211 (Jhan & Lee 2021) and Cep E (de A. Schutzer et al. 2022).The emission in the velocity range of 20-50 km s −1 is seen only in the southern region in the red-shifted lobe while no emission is seen in the northern region.The fast CO emission (the middle panel of Figure 4) also shows that the red-shifted fast CO emission is extended along southern part while stunted along the northern part near the protostar.The blue-lobe shows an opposite trend.Further- more, the red-shifted EHV SiO knots reveal a wiggling pattern from north to south as shown in the red curve in the bottom panel of Figure 1.Therefore, it is plausible that the environment toward the northern part of the red-shifted flow was previously entrained by the earlier north-east jet and outflow, producing the observed asymmetric emission.
Around the HB2 knot (-2.′′ 6, middle right panel of Figure 7), the SiO emission is similar to that in HB1e, but with smaller jet velocity.The north (∆Dec.=+0.′′ 3) and south (∆Dec.=-0.′′ 4) CH 3 OH emission (green contours) may trace the wings of the bow shocks LB1 and LB2, respectively, with velocities similar to the systemic velocity.Around HR2 (+2.′′ 6), the knot HR2a might accelerate the north wall traced by the fast SiO and CO emission (see the fast CO outflows in the middle panel of Figure 4).The HR2b knot may be interpreted as a southern accelerated component.It is also possible, however, for HR2b to be associated with Outflow2 because it is on the trajectory of the bow shock LR2, as shown in Figure 1.
The southern SiO component indicates the apex of LB1(-4.′′ 1, top right panel of Figure 7), and is surrounded by two CH 3 OH components at the systemic velocity.Northern SiO emission appears to be associated with LB2.Around the apex of LR1 (+5.′′ 2, top left panel), we also see the terminal knot HR5 and entrained components.Weak SiO knots toward the south are associated with LR2.

JET MODELS
A simple jet model provides an estimate of the jet physical parameters (e.g., Raga et al. 1990;Lee & Sahai 2004;Jhan et al. 2022).The regular spaced knots can be produced by a periodic velocity variation in which fast-jet material catches up with slow-jet material.This is in agreement with the position velocity diagram which shows higher velocities toward the central source as shown in HR2a, HR4, HB2, then HB5L in Figure 6.
A time dependent periodic jet velocity V (t) with a fixed period (P ), mean jet velocity (V j ), and velocity variation amplitude(∆V j ) can be written as: With this model, the distance (D) to the first knot is and the interval (∆D) between the knots is The observed distance to the first knot (D obs ) and the observed interval of the knots (∆D obs ) are respectively, where i is the inclination angle of the jet to the plane of the sky.
The observed jet velocity (V obs ) and the observed sideways-ejection velocity (∆V obs ) are respectively, when no energy dissipation is assumed and the sideways-ejection velocity is equal to ∆V j .Therefore, the parameters of the jets (V j , ∆V j , P , and i) can be derived from the observed parameters (V obs , ∆V obs , D obs , and ∆D obs ).For example, the inclination is derived as The jet model is applied to the red-shifted EHV SiO knots HR2a and HR4.Both knots exhibit the velocity pattern predicted by the jet model, with the higher velocity in the upstream side and lower velocity in the downstream side, as previously mentioned.The knot HR3 is ignored because it does not show such a characteristic velocity pattern.In addition, the intensity of HR3 is lower than that of HR2a and HR4.Note that the knot HR1d seems to have the velocity pattern predicted by the jet model similar to HR2a and HR4, and might therefore be considered as a potential candidate for the first knot.However, the distance to the HR1d (∼0.′′ 8) is shorter than ∆D obs (1.′′ 8) while the Lee et al.D obs is generally longer than ∆D obs (Jhan et al. 2022).In addition, the knot HR1 could be formed by the periodic ballistic ejection as mentioned in Section 4. Therefore, the knot HR1d is not used as the first knot in our jet model.The observed parameters provide the jet parameters as shown in Table 1.Thus, knots HR2a and HR4 appear to have been ejected ∼80 and ∼130 years ago, respectively.If HR1a (0.20 ′′ ), HR1b (0.48 ′′ ), and HR5 (5.5 ′′ ) have similar inclination and jet velocity, they were ejected ∼6, ∼13, and ∼150 years ago, respectively.The former two time scales are consistent with the dynamical time scales of HR1a (5.7 years) and HR1b (13 years) derived from the ballistic ejection in Section 4. Note that when the HR3 is also considered in the jet model, the inclination (20 • ) and the period (13 years) decrease while the jet velocity (155 km s −1 ) increases.The dynamical time scales of knots are reduced by half compared to the jet model without HR3.Then HR1a was ejected in the minimum phase of flux variation (see Figure 4 of Yoon et al. 2022).Thus, we prefer the jet model without HR3.
The knots HB2 and HB5 could be counterparts of the red-shifted SiO jet.Although their velocity is decreasing toward the terminal knot HB5, the projected distances from the continuum peak are similar to those of the counterpart.When the inclination and jet velocity in the red lobe are adopted, HB2 and HB5 were ejected 72 and 120 years ago, respectively.

SiO and CO
The CO 3-2 and SiO 8-7 lines constrain the physical condition within the HR1a knot, which lies 0. ′′ 2 from the continuum peak and appears to have been ejected ∼6 years ago.Using the CASA task, 'imfit', the high-resolution SiO 8-7 image (see Figure 3) shows a spatial width of the HR1a knot of 39 ± 10 mas (16.7 ± 4.3 au at the distance of 428 pc), a similar width as that seen in other jets (see Figure 7 in Lee et al. 2020).
We derive the SiO and CO column densities using a Large Velocity Gradient (LVG) model (van der Tak et al. 2007;Lee et al. 2013).The radiative and collision coefficients for CO and SiO are provided by the LAMDA database (Schöier et al. 2005;Yang et al. 2010;Balança et al. 2018).
A simple approach tells us that the HR1a knot appears to be dense clump.The SiO 8-7 spectrum at the location of HR1a shows a peak intensity of 60 K and a line width (FWHM) of 20 km s −1 (see the red spectrum in the top panel of Figure 8).Convolved to the beam size of the CO 3-2 image (0. ′′ 27), the intensity of the SiO emission (blue spectra) is a factor of 5 lower than that of the original (red) SiO spectrum due to beam dilution.The EHV gas traced by CO 3-2 and SiO 8-7 lines show similar spatial distribution in HH211 (Jhan & Lee 2021), so we assumed that the CO 3-2 and SiO 8-7 lines Figure 8. Spectra of CO 3-2 and SiO 8-7 at knots HR1a (top), HR1b (middle), and HB5L (bottom).The low-resolution (0. ′′ 26) SiO spectra (blue) are the results of convolving the high-resolution (0. ′′ 08) SiO spectra (red) with the same beam as the CO 3-2 spectra (black) for comparison.The blue and black spectra are multiplied by 5, 4, 2.5 from top to bottom.
trace the same gas.The EHV CO 3-2 emission from the knot, corrected by the beam dilution factor using the ratio of the red and blue SiO 8-7 lines, has an intensity of ∼50 K as shown by the black spectrum in the top panel of Figure 8. Assuming dense (∼10 7 cm −3 ) and warm gas (100-500 K) conditions, the CO column density is 1.1-2.9×10 18cm −2 .Thus, if the CO abundance is 10 −4 and the length of the line of sight is the same as the knot width (16.7 au), then the derived gas density is 0.4-1.2× 10 8 cm −3 , which is consistent with our initial assumption.Further, the mass of HR1a is 0.4-1.1 ×10 −5 M ⊙ and the implied mass ejection rate is 1.2-3.2×10 −6 M ⊙ yr −1 , when the length of HR1a (123 ± 13 mas), the inclination (37 • ), and the jet velocity (91 km s −1 ) are adopted.
The abundance ratio of SiO with respect to CO is different between the EHV and the slow components.At HR1a, the EHV SiO column density is 2.7-3.2 × 10 15 cm −2 , and the abundance ratio of SiO to CO is 14-5 × 10 −3 (100 K -500 K).At HR1b (the middle panel of Figure 8), where the line width is ∼10 km s −1 , the column density of SiO is 1.3-1.6 × 10 15 cm −2 , and the abundance ratio of SiO to CO is 2-10 × 10 −3 (100 K -500 K).The EHV knots (HR2a to HR5) show a similar line width (∼10 km s −1 ) and a similar intensity ratio between SiO 8-7 and CO 3-2 within a factor of two.
For the slow component (the bottom panel of Figure 8), the abundance ratio decreases by an order of magnitude.The gas density and temperature of the slow component might be between those of the EHV component and those derived by the excitation analysis of CH 3 OH (see Section 6.2).We applied the same approach used in the EHV SiO component above since the emitting area of the slow SiO component and CO emission at HB5L are similar.The gas temperature of the slow component is likely higher than 50 K because the peak intensity of the beam-corrected CO is ∼ 60 K (see the bottom panel of Figure 8).Thus, we adopt the gas density of 10 7 cm −3 and the same gas temperature of 100-500 K as in HR1a.With this setup, the SiO column density is 2.8-4.8 × 10 14 cm −2 and the abundance ratio of SiO to CO is 1.3-4.6 × 10 −4 .We note that the red-shifted slow SiO component HR5L is not considered in this analysis because the CO data are too noisy.The EHV component has a higher abundance ratio of SiO to CO than the slow components, which is also seen in Ser-emb 8 (N) (Tychoniec et al. 2019), and the ratios are similar to what is observed in other protostars (Hirano et al. 2010;Podio et al. 2021;Dutta et al. 2022a,b).
The SiO could be formed in the EHV gas launched from a dust-poor or dust-free inner disk, where the Si in the dust grains has al-ready been fully liberated into the gas phase, as appears to be the case for other protostellar sources (Lee et al. 2020;Podio et al. 2021).In that case, the abundance ratio of SiO to CO of ∼5 × 10 −3 could be reproduced if the mass ejection rate is ∼10 −6 M ⊙ yr −1 for the dust-free case or ∼10 −7 M ⊙ yr −1 for the dust depletion (see Figure 12 in Tabone et al. 2020).These ejection rates are consistent with rates estimated above for HOPS 373SW (∼10 −6 M ⊙ yr −1 , see the mass ejection rate derived with HR1a).
The abundance of SiO in the slow-moving gas is a factor of 10 lower than that in the EHV gas.For the gas density and temperature derived from the excitation condition, it takes >3000 years for the liberated SiO abundance to decrease by a factor of 10 (De Simone et al. 2022) due to freeze-out onto grain, which is much longer than the dynamical time scale (∼100 years).Thus, the slow-moving SiO component is unlikely to have been originally produced within the jet, where the SiO abundance was high.Instead, the observed slowmoving SiO likely originated from grain sputtering within the bow shock.The grain sputtering in a C-shock could produce the slow-moving SiO emission (Gusdorf et al. 2008).This low SiO abundance could be produced within a few tens years according to C-shock models (Gusdorf et al. 2008).

CH 3 OH
Figure 9 shows the CH 3 OH 7 k -6 k spectra at the CH 3 OH emission peaks (CR1-3 and CB1-3) presented as white crosses in Figure 2. The multiple CH 3 OH lines are fitted by the Large Velocity Gradient (LVG) model (van der Tak et al. 2007;Lee et al. 2013) using the radiative and collision coefficients for CH 3 OH provided by the LAMDA database (Schöier et al. 2005;Rabli & Flower 2010).We assume that the abundances of A-and E-type CH 3 OH are the same.The results for the six CH 3 OH peaks are listed in Table 2. Except for CB1 and CR1, most CH 3 OH emitting regions have a CH 3 OH column density of ∼10 15 cm −2 , a gas density of 10 7 cm −3 , and a gas temperature of ∼50 K.
a The unit of the emitting area is relative scale with respective to the beam size of CH 3 OH.
The protostellar mass is an important parameter needed to understand the separation and period of any potential binary.The protostellar mass is derived from the Keplerian rotation profile of the gas emission.A rotation feature is shown in the 13 CH 3 OH (Lee et al. 2023a).However, the derived stellar mass (<0.001M ⊙ ) is much lower than the theoretical mass of the hydrostatic core (0.01-0.05 M ⊙ Larson 1969;Penston 1969;Masunaga et al. 1998;Saigo & Tomisaka 2006) when the Keplerian motion and the inclination of 37 • are adopted.Thus, the observed velocity features near HOPS 373SW might trace the infalling rotating envelope, which has a velocity profile with a steeper power law index (1/r) than that of a Keplerian disk (1/ √ r), and HOPS 373SW might have a very small, as yet unobserved, disk (< a few au).
The protostellar mass can be roughly estimated using an indirect method.Class 0 protostars have a similar ratio of mass ejection rate to mass accretion rate, Ṁj / Ṁacc , (∼0.19, Nisini et al. 2015;Lee et al. 2020;Podio et al. 2021).The recent mass ejection rate derived in Section 6.1 is 1.2-3.2×10 −6 M ⊙ yr −1 , which results in a mass accretion rate of 0.6-1.7 ×10 −5 M ⊙ yr −1 .Taking this value for the mass accretion rate, the stellar mass can be estimated from where the bolometric luminosity of 5.3 L ⊙ is adopted (Kang et al. 2015;Yoon et al. 2022).The stellar mass is then found to be between 0.025-0.07M ⊙ .
The derived stellar mass provides a consistent condition launching the SiO jets.The jet-launching radius (R jet ) is estimated as (Lee et al. 2020), where G is the gravitational constant.The derived R jet is ∼0.03 au, which is consistent with that in other jets (Lee et al. 2020).The jetlaunching radius is much smaller than the dust sublimation radius (R s ) of 0.16 au, implying that the jets are launched from the dust-free region as mentioned in Section 6.1.The dust sublimation radius is calculated as (Monnier & Millan-Gabet 2002) (11) where Q R is the ratio of dust absorption efficiency at the color temperature of the incident and re-emitting radiation fields and T sub is the dust sublimation temperature.Here, we adopt Q R = 4 and T sub = 1500 K (Dullemond et al. 2001).Note that the values of M star , R jet , and R s could be smaller than the values estimated here if the bolometric luminosity of HOPS 373SW itself is lower than 5.3 L ⊙ , which is the total luminosity of HOPS 373 including the SW and NE sources.

Dynamical time scales for the two outflows
Outflow1 producing the U-shaped bow shocks LR1 and LB1 is younger than Outflow2 producing LR2 and LB2 (see Figure 4).The EHV jets and the H 2 emission in the near-infrared are associated with Outflow1 (see Figure 5).The jet model shows the knot HR5 at the tip of LR1 in the end of Outflow1 was ejected 150 years ago (see Section 5).
The dynamical time scale of Outflow2 can be derived indirectly.The peak intensity and velocity of the SiO knot at the LB2 tip are ∼40 K and -3.5 km s −1 , respectively, and those of the marginally detected weak SiO knot at the LR2 tip are ∼19 K and 4.5 km s −1 , respectively.The peak temperature of 19-40 K is produced by a ∼30 km s −1 shock if the pre-shock density is 10 6 cm −3 (Gusdorf et al. 2008).In addition, in LR2 and LB2 the velocity of the peak SiO emission is slightly faster than that of the peak CH 3 OH emission by 0.6 km s −1 and 2.4 km s −1 , respectively.The CH 3 OH is quickly destroyed through reactions with H atoms after sputtering from the grain surface in the high-velocity shock (>25 km s −1 , Suutarinen et al. 2014;Nesterenok 2018).Thus, the SiO and CH 3 OH knots would trace the tip and wings of the bow shocks, respectively, as shown in the bow shock LB1, and they might be overlapping each other due to the slow shock velocity.Using the above estimated shock velocities the dynamical time scale of the Outflow2 can be estimated.For a shock velocity of 30 km s −1 , the inclination of the outflows is derived as 7-9 • from the velocities at the peak intensities of LR2 (4.5 km s −1 ) and LB2 (-3.5 km s −1 ), and the LR2 and LB2 would have been ejected ∼500 and ∼370 years ago, respectively.

Episodic accretions
The EHV SiO jet shows two periodic variation timescales.The jet model reveals a period of 50 years with an observed separation of 1. ′′ 9 along the HR2a and HR4 knots.The higher resolution image shows that three sub-components lie within HR1 and HR2a (as shown in Figure 1) where the separation is 0. ′′ 3, a factor of 6 smaller than the larger scale variation and implying a period of ∼8 years when the jet velocity and the inclination are assumed the same as the jet model.
These periodic variations may be related to a periodical perturbation of the accretion in the disk by a close binary companion (Lee et al. 2020).Adopting a protostellar mass of 0.035 M ⊙ , the variation with ∼8 years and ∼50 years could be induced by events at 1.3 au and 4.4 au, respectively, which cannot be resolved by the present observations where the major axis of the continuum emission (FWHM) is 71.1 ± 2.6 mas (30.4 ± 1.1 au).L1157 MMS, with a separation of 16 au, is the tightest known binary among the directly observed Class 0 protostars (Tobin et al. 2022).Such small scales required detection via VLA continuum observations.Interestingly, the wiggle in the Class 0 protostellar jet HH 211 implies a binary separation of ∼4.6 au (Lee et al. 2010).Theoretical modeling finds that the circumstellar disk is truncated by interaction with a binary such that the disk radius is 0.25-0.5 the binary semi-major axis distance Lee et al. (Terquem et al. 1999).Thus, for HOPS 373SW the 1.3 au length scale might be associated with the disk radius while the 4.4 au scale represents the separation of the binary.
This putative close binary could reproduce the wiggle pattern of the SiO EHV emission shown in Figure 1.When a protostar has a mass of M p , a disk radius of r d , and a companion is separated with a distance of a and has a mass of M c , then the ratio of precession period τ p to orbital period τ o is given by : ) where β is the half-openning angle (Terquem et al. 1999;Lee et al. 2020).Thus, τ p is ∼15 times longer than τ o assuming M p /M c =1.The wiggling pattern is described by the parametric equation using β, V j • τ p (the spatial period), the phase angle (ϕ), and the inclination (see Eqs (1)-(3) in de A. Schutzer et al. 2022).The wiggle motion is roughly reproduced by the precession model assuming M p /M c =1 and using the parameters V j , P = τ o , and i from the jet model (Table 1), β= 8 • , and ϕ= 50 • as shown in the red and blue curved lines in the bottom panel of Figure 1.
According to the relic record provided by the observed knots, many burst events should have occurred in the past.Figure 10 shows an anti-correlation between the C 18 O 2-1 and N 2 D + emission.
The C 18 O emission has the Full Width Half Maximum (FWHM) of 4. ′′ 63 ± 0. ′′ 61 which is measured by the CASA task, 'imfit'.The expected FWHM from the current luminosity (5.3 L ⊙ ) is 1. ′′ 8 when the CO sublimation temperature is 21 K (Frimann et al. 2017).Thus, the current CO emission is a factor of 2.6 more extended than that expected from the current luminosity.In addition, the peak of the N 2 D + is 6. ′′ 3-7.′′ 4 off from the center of HOPS 373SW, implying the luminosity during a past burst reached ∼40-60 L ⊙ (Hsieh et al. 2019;Murillo et al. 2022).
We note that HOPS 373SW is very likely the dominant trigger for the extended CO snowline because the distribution of N 2 D + emission implies the heating source is closer to the HOPS 373SW rather than HOPS 373NE, and the C 18 O 2-1 emission peak is also close to the HOPS 373SW.Furthermore, the orientation of the two N 2 D + peaks is approximately perpendicular to the large-scale CO outflow in the SE-NW direction, implying that the large-scale CO outflow is also associated with HOPS 373SW rather than HOPS 373NE.
The freeze-out time scales for CO and water provide an upper limit for when this large burst occurred.The number density of molecular hydrogen around the observed CO snowline is ∼5 × 10 5 cm −3 when the density profile for HOPS 373 presented in Furlan et al. ( 2016) is adopted.At this density, the freeze-out time scale for CO is ∼10 4 years (e.g., Lee et al. 2004).On the other hand, the water snowline, which is derived by the boundary of emission distribution of complex organic molecules (∼0.′′ 1), is similar to that derived from the current luminosity (Lee et al. 2023a).The water snowline for a large burst (40-60 L ⊙ ) is ∼0.′′ 23 (100 au at the distance of 428 pc), where the density is higher compared to that at the CO snowline and the freeze-out time scale for the water is ∼100 years (van 't Hoff et al. 2022;Lee et al. 2023a).Thus, we expect that this recent burst occurred between the freeze-out time scales for the CO and the water (between ∼10 2 and ∼10 4 years).The dynamical time scales of Outflow1 (150 years), Outflow2 (∼500 years), and the large scale CO outflow in the SE-NW direction (∼10 4 years, Nagy et al. 2020) are all within the estimated time range of the accretion burst that induced the observed CO snowline.Thus, the strongest burst associated with one of the three outflows likely shaped the observed CO snowline.Fur- ther observations and studies are needed to determine which outflow is associated with this strongest burst event.

Origin of multiple Outflows
The SiO, CH 3 OH, and CO lines toward HOPS 373SW reveal two outflows: Outflow1 (NE-SW) and Outflow2 (SE-NW).In addition, the large-scale CO outflow in the SE-NW direction was also ejected from HOPS 373SW.A possible explanation for these features is that the three outflows were successively ejected from a single protostar.Initially, the large-scale CO outflow had been ejected ∼10 4 years ago (Nagy et al. 2020).Then, Outflow2 was ejected ∼370-500 years ago and lastly Outflow1 was formed ∼150 years ago after the system went through a reorientation, as appears to be the case for IRAS 15398-3359 (Okoda et al. 2021).Given that HOPS 373 is an early stage Class 0 protostar (Kang et al. 2015;Furlan et al. 2016;Stutz et al. 2013), it expected that the angular momentum axis of the accreting gas might change over time due to in-homogeneities in the star-forming core (Sakai et al. 2018;Zhang et al. 2019;Gaudel et al. 2020;Okoda et al. 2021).In addition, the outflow reorientation also occurs when the rotation axis of the core is misaligned with the global magnetic field direction (e.g., Machida et al. 2020).The extended CO snowline suggests that a strong accretion burst occurred 10 2 -10 4 years ago.Therefore, the three observed misaligned outflows may have been induced by accretion bursts accompanied by reorientation.
A second possible explanation for the two outflows is that they are evidence of a tight protostellar binary.In this case as well, the change in the outflow axis from the large scale CO outflow in the SE-NW direction to Outflow1 or Out-flow2 occurred by an accretion burst accompanied by reorientation.The EHV SiO jet shows the wiggling motion, which could be due the precession induced by a close-binary system as mentioned above.Quasi-episodic knots with a period of 50 years also support this hypothesis (Lee et al. 2020).As mentioned above, it is possible that the knot HR2b is associated with Outflow2 and that the CH 3 OH emission peaks CR2, CR1, and CB1 are on the trajectory of the bow shocks LR2 and LB2.If those CH 3 OH peaks are associated with Outflow2, then CR1 and CB1 would have be ejected ∼130 years ago, similar to the time scale of the terminal knots (∼150 years) in Outflow1.In addition, the two outflows have different spatial properties.They are misaligned by ∼30 • and the jet velocities are ∼90 km s −1 and ∼30 km s −1 , respectively.Therefore, each outflow might have been ejected by a different member of a close-in protostellar binary.
Such a misaligned close binary could be formed through core fragmentation.Misaligned outflows with an angle of 70 • have been reported in the Class 0 protostellar binary system VLA 1623 A, where the binary separation is 36 au (Hara et al. 2021).The binary separation in the HOPS 373SW is estimated to be 4.4 au (see above), which cannot be resolved with the present high-resolution continuum image (FWHM ∼34 au).Such a misaligned close binary might be formed initially at larger separations via turbulent fragmentation and have later migrated inward (e.g., Lee et al. 2017Lee et al. , 2019;;Offner et al. 2022).We note that magnetohydrodynamic simulations (Lee et al. 2019;Saiki & Machida 2020) show that close binaries (< 30 au) can be formed via core fragmentation.A recent hydrodynamical simulation for a multiple protostellar system also showed that the interaction between a close binary could also result in misaligned disks (Lee et al. 2023b).

SUMMARY
We have investigated the jets and outflows in HOPS 373SW using high-resolution ALMA observations of the SiO 8-7, CH 3 OH, and CO 3-2 lines.The results are: (1) The slow moving SiO and CH 3 OH show two distinct bow shocks which surround two separate CO outflows: Outflow1 (NE-SW) and Outflow2 (SE-NW).
(2) The extreme high velocity (EHV) and fast moving SiO emission show several knots, which could be formed in quasi-episodic jets with an inclination of 37 • , a jet velocity of 91 km s −1 , and a period of 50 years.These knots also show a wiggle motion previously noted in the CO and H 2 emissions by Yoon et al. (2022).
(3) The SiO jets show two periodic variation timescales with periods of 50 years and 8 years.The stellar mass is estimated to be 0.035 M ⊙ in which case the longer the periodic variation might be induced by a binary with a separation of 4.4 au.
(4) The observed extended CO snow line implies that a previous outburst with a luminosity of 40-60 L ⊙ occurred between ∼10 2 and ∼10 4 years ago.It is possible that the observed CO snowline is associated with either the bow shock LR1 and LB1 (∼150 years ago) or the large scale CO outflow in the SE-NW direction (∼10 4 years ago).
(5) The two observed misaligned outflows appear to have either be ejected from a single protostar, due to a reorientation of the angular momentum axis, or from protostellar binary, where each outflow is related to a different binary component.The necessary tight binary would have been formed thorough core fragmentation with a wider separation and after which it migrated inward.
Figure 1.Maps of SiO 8-7 emission at slow (top), fast (middle), and extremely high velocity (EHV, bottom).The red and blue contours indicate the red-and blue-shifted SiO emissions integrated over the velocity ranges respective to the systemic velocity presented in each panel.The start and step of contour levels are 10σ (1σ= 13, 17, and 19 mJy beam −1 km s −1 for the slow, fast, and EHV components, respectively).The background image represents peak intensity (moment 8) map of CH 3 OH, and CH 3 OH emissions below 5σ are clipped out for clarity.The green contours in the bottom panel indicate the H 2 emission (Yoon et al. 2022).The origin and black cross indicate the continuum peak positions of HOPS 373SW and 373NE, respectively.The locations of the knots are marked in the red (HR1, HR2a, ..., HR5, HR5L) and blue (HB1, HB2, HB5, HB5L) lobes.The black half-ellipses indicate the bow shocks, LB1, LR1, LB2, and LR2.The beam sizes of CH 3 OH and SiO are presented by black and red ellipses in the lower right corner.The red and blue curves in the bottom panel indicate an example wiggling pattern due to the precession of a putative close binary (see Section 7.3).

Figure 2 .
Figure 2. Map of CH 3 OH emission.The red and blue contours indicate the red-and blue-shifted CH 3 OH emissions integrated over the velocity ranges of 0.5 < |∆V | <3.5 km s −1 , and the start and step of contour levels are 10σ (1σ= 3.0 mJy beam −1 km s −1 ).The background image indicates the peak intensity (moment 8) map of SiO 8-7 emission, and the SiO emissions below 5σ are clipped out for clarity.The positions of the strong CH 3 OH emission are marked by white crosses (CR3, CR2, CR1, CB1, CB2, and CB3).The CH 3 OH emission at the origin is originated by a hot corino as well as the red-shifted outflow (Lee et al. 2023a).The beam sizes of SiO and CH 3 OH are presented by black and red ellipses in the lower right corner.

Figure 3 .
Figure3.Peak intensity (contour) and absolute velocity at the peak intensity (color) maps for the SiO 8-7 using the high angular resolution image.The maps are created as moments 8 and 9 maps using the CASA task 'immoments' after smoothing with a boxcar of 5 km s −1 .The SiO emission toward the west (∆R.A. < 0 ′′ ) is blue-shifted, but the absolute velocity is plotted for comparison with the red-shifted lobe.The start and step of contour levels are 6σ (1σ= 1.1 mJy beam −1 ).The knots HR1 and HB1 consist of several sub-knots.The emission near the origin is originated from complex organic molecules rather than SiO.The red contours indicate the ALMA Band 7 continuum image.The contours are 10, 20, 40, 80, 160, and 320σ (1σ = 0.11 mJy beam −1 ).The beam sizes of the continuum and SiO images are plotted in the red and black ellipses in the lower right corner.

Figure 4 .Figure 5 .
Figure 4. Maps of 12 CO 3-2 emission at slow (top), fast (middle), and extremely high velocity (EHV, bottom).The red and blue contours indicate the red-and blue-shifted CO emissions integrated over the velocity ranges respective to the systemic velocity presented in each panel.The contours start at 10σ with steps of 5 σ (1σ = 0.04, 0.07, 0.08 Jy beam −1 km s −1 from top to bottom panels, respectively).The background image and green contours indicate the peak intensity (moment 8) maps of the SiO and CH 3 OH emission, respectively, with 10, 20, and 30 σ levels (1σ= 1.6 mJy beam −1 ).Only the SiO emissions over 5 σ are plotted for the clarity.The black crosses indicate the continuum peaks of NE and SW sources, though the contours of the central CH 3 OH peak mostly hide the SW cross.The black curves indicate the U-shaped bow shocks, LR1 and LB1, induced by Outflow1, and LR2 and LB2 by Outflow2.

Figure 6 .
Figure 6.Position velocity diagram for SiO emission along the R.A. (jet) direction.Here the emission within |∆Dec.|<1.′′ 5 is averaged.The start and step of the contour levels are 3σ (1σ= 0.7 mJy beam −1 ).The insert panels are zoomed in the high-velocity jet near the continuum source in the high angular resolution image.The emission is averaged within |∆Dec.|<0.′′ 5, and the start and step of the contour levels are 3 σ (1σ= 0.4 mJy beam −1 ).Both the SiO spectra are smoothed with a boxcar of 3 km s −1 to increase the S/N.The emission within the green boxes at the continuum peak are associated with complex organic molecules (seeLee et al. 2023a).The yellow and red lines in the insert panels are the slopes for the ballistic ejection (see the Section 4).

Figure 7 .
Figure 7. Position velocity diagram of CO (white), SiO (red), and CH 3 OH (green) emission along the Dec. direction.The relative distance to the protostar (∆R.A.) is indicated in each panels.The image and white contours indicate the CO 2-1 emission.The red and green contours indicate the SiO and CH 3 OH emissions, respectively.The start and step of the contour levels are 3 σ (1σ = 11, 1.7, 1.6 mJy beam −1 for CO, SiO, and CH 3 OH, respectively).The vertical dashed lines indicate the source velocity.The images are smoothed with a boxcar of 3 km s −1 .

Figure 9 .
Figure 9. CH 3 OH spectra for named positions (white crosses in Figure 2) relative to the local continuum peak.The red and blue vertical lines indicate the frequencies of the A-and E-type CH 3 OH, respectively.The upper energy levels are presented near the vertical lines.

Table 2 .
LVG models for CH 3 OH emission in HOPS 373SWPosition Emitting area a N(CH 3 OH) n H 2 T gas