An Extremely Young Protostellar Core, MMS 1/ OMC-3: Episodic Mass Ejection History Traced by the Micro SiO Jet

Figure 1 shows the 1.3 mm continuum emission toward MMS 1 obtained with ALMA. We have detected centrally concentrated continuum emission both in the low- and high-angular resolution images. Note that the origin of the 1.3 mm continuum emission observed toward the protostellar sources in OMC-3 is considered to be significantly dominated by the thermal dust emission. The emission attributed to the free-free emission, tracing ionized jets, is likely negligible as estimated by Takahashi et al. (2006, 2009, 2013a, 2019).

In order to characterize the source structure, particularly the centrally concentrated emission, presumably tracing the dust disk, we performed two-component two dimensional (2D) Gaussian fitting on the high-angular resolution image. The residual level of the fitting result is less than 5.6% with respect to the observed peak flux (i.e., the residual level is less than S/N = 2.5). The source structure was fitted by an extended (${\sim}750{\pm}47\,{\rm au}{\times}510{\pm}33$ au) and a compact (${\sim}43{\pm}5.1\,{\rm au}{\times}26{\pm}4.7$ au) Gaussian component, both of which are listed in Table 2. The position angle²²2The position angle is measured from the north to the east with respect to the blue-shifted lobe. (p.a.) of the major axis of both extended and compact components ($150^{\circ}\pm 7.5^{\circ}$ and $156^{\circ}\pm 12^{\circ}$, respectively) are roughly aligned to the direction of the large-scale filament in this region (i.e., Chini et al. 1997; Johnstone & Bally 1999; p.a.${\sim}135^{\circ}$). The orientation of the compact component is almost perpendicular to the molecular jet detected with ALMA in the SiO emission (Section 3.3). Assuming that the detected compact 1.3 mm continuum emission traces the dust disk, an inclination angle of the disk ($i$) is estimated to be ${\sim}53^{\circ}$ using the axis ratio of the 1.3 mm continuum compact component. Note that the millimeter continuum fitting was also performed by Liu et al. (2023) as a part of the full polarization data analysis. The data were obtained in the 1.1 mm band. A slightly higher angular resolution than that presented here of $0$″$.14$ was achieved with a factor of 6.8 better sensitivity for the full polarization data presented. Despite the significant sensitivity differences, the total flux measured within the compact component is consistent within the factor of $\sim$1.2. The measured size of the disk-like structure in this paper is factor of $\sim$4.5 larger (in terms of the surface area) than those estimated in Liu et al. (2023). This is likely due to the factor of $\sim$1.7 larger beam size (in terms of the beam surface area) in the data set presented here compared with that presented in Liu et al. (2023) to ApJ).

The mass of the circumstellar material can be estimated from the 1.3 mm continuum emission using the following equation,

where $F_{\lambda}$ is the total 1.3 mm flux, $d$ is the distance to the source, ${\kappa}_{\lambda}$ is the dust opacity (absorption coefficient per unit gas mass). Here, a gas-to-dust mass ratio of 100 is assumed. $T_{d}$ is the dust temperature, and $B_{\lambda}(T_{d})$ is the Planck function at a temperature of $T_{d}$. Assuming, $d$ = 393 pc (distance to the OMC-3 region measured from the $Gaia$ DR2 data; Tobin et al. 2020; McBride & Kounkel 2019; Gaia Collaboration et al. 2018), ${\kappa}_{\rm{1.3mm}}$ = 0.011 cm² g^-1 from the dust coagulation model of the MRN (Mathis et al., 1977) with thin ice mantles at a number density of 10⁸ cm^-3 (Ossenkopf & Henning, 1994), $T_{\rm{dust}}$ = 30 K – 100 K (e.g., Nomura & Millar 2005; Tomida et al. 2017), and given the measured total fluxes listed in Table 2, the masses of the compact and extended components are estimated to be $M_{\rm{H_{2}(comp.)}}{\sim}$0.0051–0.020 $M_{\odot}$ and $M_{\rm{H_{2}(ext.)}}{\sim}$0.041–0.16 $M_{\odot}$, respectively. Note that the dust temperature assumption was discussed by Liu et al. (2023) in detail based on recent studies toward Orion protostars by Tobin et al. (2020) and Xu & Kunz (2021) considering the stellar radiation and disk accretion heating, respectively. The minimum and maximum temperatures adopted here are comparable to the temperature assuming the stellar radiation expected from the given radius corresponding to the extended and compact structures.

Figure 2 presents comparisons of the missing flux measured in the CO and SiO images obtained from the low resolution (${\theta}{\sim}1{\arcsec}.5$) and high resolution (${\theta}{\sim}0{\arcsec}.2$) images. The comparisons show that the fluxes obtained from the two data sets are almost the same. Thus, there is no significant missing flux in the high angular resolution data set. The missing flux comparisons imply that a majority of the flux originating from the outflow and the jet, traced by CO (2–1) and SiO (5–4), came from compact structures that are packed within the MRS of the high-angular resolution data (${\sim}2{\arcsec}.0$). Since no significant missing flux was reported particularly, in the mid- and high-velocity ranges for both CO and SiO emissions, we only use the high angular resolution images in this paper.

We note discrete flux enhancements both in CO (2–1) and SiO (5–4) at the high velocity range of $|v_{\rm{LSR}}-v_{\rm{sys}}|{\approx}50$ km s^-1. These line profiles show clear evidence of the extremely high velocity (EHV) flow which was first discovered by Bachiller & Cernicharo (1990) toward NGC 1333 (HH 7-11), and then was reported for several low- and intermediate-mass Class 0 sources (Bachiller et al., 1991a, b, 2000; Zapata et al., 2005; Hirano et al., 2010; Gómez-Ruiz et al., 2013, 2019; Tafalla et al., 2017; Lee et al., 2017; Matsushita et al., 2019; Lee, 2020). Together with recent studies of the compact EHV flow toward MMS 5 by Matsushita et al. (2019), and MMS 6 Takahashi et al. (2012), Takahashi et al. (2019), and Takahashi et al. (2024 in prep.), the EHV flow associated with MMS 1 is one of the most compact EHV flows ever reported (see Section 3.3). Note that, hereafter, the terminology of “jet (or SiO jet)” will be used to refer to the EHV flow.

We present the ALMA CO (2–1) and SiO (5-4) emission arising from the outflow and jet located in MMS 1 in Figure 3. A compact molecular outflow is traced by the CO emission, and a collimated high-velocity jet within the outflow is traced by the SiO emission. Both outflow and jet associated with MMS 1 were, for the first time, detected thanks to the ALMA’s sensitivity and high-angular resolution capabilities. This clearly indicates that MMS 1 contains a protostar. The blue- and red-shifted gas are located in the southwest (SW) and northeast (NE) directions with respect to the millimeter source peak, respectively, with the position angle of -137^∘. As seen in Figure 3b, the axis of the CO outflow and SiO jet is aligned more or less perpendicular to the disk-like structure traced in the 1.3 mm continuum emission. As clearly seen in Figure 3a, the jet traced in the SiO emission shows a wiggled structure.

Figures 4a, b, and c present moment maps obtained from the CO (2–1) emission. The CO blue- and red-shifted lobes are extended up to ${\sim}3{\arcsec}.2$ ($\sim$1260 au) and ${\sim}5{\arcsec}.1$ ($\sim$2000 au), respectively. The CO emission is detected with the LSR velocity range between -35 and 60 km s^-1. Figure 4c shows that the majority of the gas detected in CO (2–1) emission has a velocity dispersion of $|v_{\rm{LSR}}-v_{\rm{sys}}|{\lesssim}7$ km s^-1. An X-shaped structure associated with central protostar, tracing the outflow cavity is clearly seen. The width of the outflow, where the outflow has the widest width in both lobes, is measured to be ${\sim}1{\arcsec}.3$ ($\sim$510 au).

As presented in Figures 4d, e, and f, the collimated jet is detected in SiO (5–4) emission with the projected length of ${\sim}2{\arcsec}.3$ (${\sim}900$ au; blue-shifted emission) and ${\sim}1{\arcsec}.3$ (${\sim}510$ au; red-shifted emission), respectively. The SiO emission was detected with the LSR velocity between -60 and 70 km s^-1. Velocity tendency is consistent with that observed in CO, while the SiO emission mainly tracing the high-velocity collimated jet rather than the wide-opening angle outflow seen in CO. Note that the SiO emission is located close to the central star, in particular the blue-shifted component shows a width comparable to the CO outflow (Figure 3b and 4e). This wide width SiO emission distribution is only seen in the relatively low-velocity (the LSR velocity between -15 km s^-1 and 5 km s^-1), tracing the outflow cavity. Within the outflow cavity, we see a high-velocity collimated jet.

The detected high-velocity SiO emission traces two components: (i) emission associated with the central region, showing a collimated structure and (ii) bullet-shaped emission (i.e., bow-shocks as explained using Figure 5 in the following paragraph) located at the tips of the CO outflow. The first component traces a collimated jet and is slightly brighter on the SW side (blue-shifted emission) than on the NE side (red-shifted emission). The spatially resolved image (Figure 4d) reveals that the jet appears to wiggle, and that the SW side of the jet contains at least three bright knots (and an additional faint knot as presented in Figure 6). The second component, showing the bullet-shaped structure, is located at ${\sim}6{\arcsec}.2$ ($\sim$580 au; blue-shifted emission) and ${\sim}2{\arcsec}.0$ ($\sim$790 au; red-shifted emission) away from the 1.3 mm continuum peak position and brighter in the NE side (red-shifted emission) than in the SW side (blue-shifted emission). The location of the SiO bullet-shaped structure roughly coincides with the tips of the outflow lobes observed in the CO emission. The associated gas velocity is lower than that observed in the jet. The blue-shifted emission is detected in the LSR velocity range between -10 and -40 km s^-1, and the red-shifted emission is detected in the LSR velocity range between 10 and 30 km s^-1 (Figure 4 and 5). Velocity dispersions of ${\Delta}v{\sim}$10 km s^-1 and ${\Delta}v{\sim}$4 km s^-1 are observed toward the blue-shifted and the red-shifted bullet-shaped emission, respectively (Figure 4f).

Figure 5 presents the position-velocity diagram (PV-diagram) obtained from the CO (2–1) and SiO (5–4) data cubes cut along the outflow axis (p.a. = -137^∘). Clearly seen in the PV-diagram, there are two distinct components: (i) a collimated jet and (ii) the bullets located at the tips of the CO outflow. Regarding the first component (i), SiO emission is concentrated in the central $2{\arcsec}$ region, showing a linear velocity increase with respect to distance from the center. The right panel of Figure 5 presents the zoomed image of the jet, which clearly shows that there are at least three components (as denoted by colored lines) showing similar respective velocity increases. This is particularly clear in the blue-shifted jet (i.e., SW side). The terminal point of each component coincides with an emission peak, which also corresponds to individual knots indicated in Figure 4d. Regarding the second component (ii), the SiO emission is located at the tips of the CO blue- and red-shifted outflows. The spatial and velocity patterns in the PV-diagram are consistent with those explained with the jet bow-shock model presented in Figure 2 of Arce et al. (2007). Consistent with what we see in Figure 4b and c, the CO emission in Figure 5 mainly traces the low- to intermediate-velocity gas components ($|v_{\rm LSR}-v_{\rm sys}|{\lesssim}$20 km s^-1). These components are spatially extended and mostly associated with the outflow cavity and lobes. The CO emission is also detected toward the bow-shocks, and is particularly bright in the red-shifted emission component (i.e., NE side).

In this subsection, we estimate two types of timescales: (i) outflow and jet dynamical timescales, which represent how young the protostar is, and (ii) dynamical timescales to each knot, which is possibly related to the time interval of the episodic mass ejection within the jet.

We showed evidence for jet wiggling in Section 3.3 (Figure 3a). In general, jet wiggling can be explained either by orbital motion of a binary system (Masciadri & Raga, 2002) or by jet precession (or precession of the jet-driving object; Frank et al. 2014). Both models produce a “garden-hose” effect, resulting in jet wiggling (Raga et al., 1993). The observed MMS 1 jet morphology allows us to distinguish the origin of the jet wiggling (Frank et al., 2014; de A. Schutzer et al., 2022, and references therein). In the case of jet wiggling due to binary orbital motion, the jet is expected to bend in the direction opposite to the location of the secondary (binary) component. This produces a V-shaped jet in the vicinity of the protostar. The previously ejected gas propagates within the cone-shaped structure, and the side-to-side wiggling over time produces a W-shaped structure in the projected jet locus. Alternatively, in the case of jet precession, the jet is expected to be ejected symmetrically. Side-to-side moving gas within the cone therefore produces an S-shaped locus over time.

The SiO jet ejected from MMS 1 shows an S-shaped wiggle. Furthermore, at the current angular resolution of ${\sim}0^{\prime\prime}.15$ ($\sim$60 au in the linear size scale), we did not find any evidence that MMS 1 is a binary system. These facts support that the jet wiggling associated with MMS 1 can be caused by the jet precession rather than the orbital motion of a binary system.

We fitted the wiggling jet observed in MMS 1 with a three dimensional spiral morphology using the following parametric equations (3), (4), and (5), as described byde A. Schutzer et al. (2022) and originally based on Raga et al. (1993) and Masciadri & Raga (2002);

where we consider the situation that the jet precesses inside a cone of main axis $z$. Here, $\beta$ is the jet half-opening angle, $\lambda$ is the spatial period, and $\phi$ is the phase angle. The precession axis tilts at an inclination angle of $i$ with respect to the plane of the sky. The projected coordinates (${\alpha}^{\prime}$, ${\delta}^{\prime}$) in this plane are given by

We adopted the inclination angle of the jet as 45^∘ which is estimated from the disk axis ratio, as discussed in Section 3.4.1. Then, the other parameters of the jet position angle $\lambda$, $\beta$, and the phase angle $\phi$ were fitted simultaneously by eye. During the fitting process, we found that the northern and southern jets have different position angles of 35^∘ and 51^∘, respectively, requiring an axis misalignment of 16^∘. Hence, the fitting was done independently for the northern and southern jets. The fitting results and parameters are presented in Figure 7 and Table 3, respectively. Although the axis is misaligned between the northern and southern jets, the fitted parameters are similar, indicating that the jet shows a more or less symmetric ejection.

Although the cause of the jet precession is not fully understood, a few possible scenarios have been proposed. The first one is the tidal interaction between the associated protostellar disk and the noncoplanar binary companion, creating an S-shaped wiggle jet in the vicinity of protostars (Terquem et al., 1999). However, we did not find any apparent substructures within the central 100 au scale with the current angular resolution, as shown in Figure 1b. Thus, we can exclude the tidal interaction scenario.

Alternatively, recent simulation studies have demonstrated that the misalignment of the core rotation axis with respect to the global magnetic field can produce an S-shaped wiggling jet (Hirano & Machida, 2019; Hirano et al., 2020; Machida et al., 2020). In their core-collapse simulations, the precession of the protostar and disk occurs when the prestellar core has a rotation axis that is not parallel with the global magnetic field. In our observations, a slight misalignment between the large-scale magnetic field (p.a. $\sim$45^∘) and both the disk rotation axis (i.e., perpendicular to the disk major axis, p.a. $\sim$65^∘) and the observed SiO jet axis (p.a. 35^∘ for the northern jet and p.a. 51^∘ for the southern jet) has been observed toward MMS 1, which could explain the jet precession observed in MMS 1.

To further validate the acceleration of the jet in MMS 1 seen in Figure 5, we compare a recent numerical simulation of protostar evolution with our results. The data in Figure 8 were taken from the core collapse simulation performed by Machida & Basu (2019).

We first explain the numerical settings and detailed simulation results by Machida & Basu (2019), for which a spherical core with a Bonnor-Ebert density profile was adopted as the initial state. The radius $R_{\rm cl}$ and mass $M_{\rm cl}$ of this initial cloud core are $R_{\rm cl}=1.2\times 10^{4}$ au and $M_{\rm cl}=1.7\,M_{\odot}$. A uniform magnetic field with $B_{0}=4.5\times 10^{-5}$ G and a rigid rotation of $\Omega_{0}=1.1\times 10^{-13}$ rad s^-1 were added to the initial state. The gravitational collapse of the core was calculated using a nested grid method (Machida et al., 2004; Machida & Hosokawa, 2013), in which the resistive magnetohydrodynamics equations (see, eqs. [1]–[7] in Machida & Matsumoto 2012) were solved. The finest spatial resolution of the simulation is $5.6\times 10^{-3}$ au. As a result, both the star-forming core ($\sim 10^{4}$ au) and the protostar ($\sim 0.01$ au) were spatially resolved in the simulation. Starting from the prestellar core, the evolution of the system is calculated for 2000 yr after protostar formation.

The simulation produces episodic jets driven from the inner disk region (see Fig. 2 of Machida, 2014). These jets with speeds $\sim 100$ km s^-1 appear repeatedly every $\sim 1-10$ yr until the end of the simulation 2000 yr after protostar formation. In the following analysis, we consider the simulation result between 200.21-201.69 yr after protostar formation, during which the high-speed jet appears three times, as shown in Figure 8. Although the full simulation reveals a wide range of mass ejection rates due to the jets ($\sim 10^{-6}-10^{-9}$ M_⊙ yr^-1, see Fig. 7 of Machida & Basu 2019), the mass ejection rate during the epoch shown in Figure 8 is about $3\times 10^{-6}$ M_⊙ yr^-1.

In Figure 8, the density and velocity distributions on the $y=0$ plane (left) and the PV diagram along the jet or $z-$axis (right) at three different epochs, 200.21, 200.58, and 201.69 yr after protostar formation are plotted. As described above, the simulation spatially resolved the protostar and circumstellar disk and reproduced the outflow and jet. Thus, we can qualitatively compare the kinematics of the observed jet with the simulation. The left panels of Figure 8 show that the high-velocity flow, with several local density peaks, is surrounded by the low-velocity flow (for details, see Machida & Basu 2019). The right panel of Figure 8 plots the PV-diagram at each epoch, in which we only used the simulation data around the z-axis as -5 au ${\leq}x{\leq}$ 5 au and -5 au${\leq}y{\leq}$ 5 au to create the PV-diagram (the white dashed line in panels (a), (c), and (e)) in order to focus on the high-velocity component (i.e., the jet).

In the second and fourth quadrants of the PV diagram in Figure 8b, we can confirm several spine-like components. Each spine corresponds to a clump ejected from the region near the protostar. We can see a weak density peak at $z\simeq{\pm}20$ au within the white dashed line of Figure 8a which corresponds to the clump produced during a previous mass ejection event occurring $197.18$ years after protostar formation. In the PV diagram (Figure 8b), we can also confirm the jet intermittent components at $z=\pm 20$ au, corresponding to the clump made during the previous ejection event, which is labeled as [1]. In addition, a strong mass ejection is currently occurring at the root of the jet. A sharp spine around $z\simeq 0$ in the PV diagram (Figure 8a) labeled as [2] corresponds to the current mass ejection event close to the protostar. The color in the PV diagram representing the launching radius of each velocity component indicates that the jet launching radii are widely distributed across the range of 0.01-1 au. We confirmed that various velocity components in the range $<150$ km s^-1 appear at every mass ejection event. The mass ejection occurs at different disk radii and the gas ejected at different radii has different velocities. Thus, a sharp spine appears around $z=0$ in the PV diagram immediately after the mass ejection occurs. Note that mass ejection recurrently occurs in the main accretion phase (Machida & Basu, 2019).

The spine labeled as [2] widens and is inclined from the vertical axis in Figure 8(d). In addition, we can confirm that the jet velocity labeled as [2] linearly increases as the distance from the protostar increases in Figure 8(f), which is similar to the Hubble-like velocity structure or Hubble flow. As described above, the ejected gas has different velocities and the gas with different velocities propagates along the $z-$axis. Thus, the highest-velocity gas is located farthest from the protostar as time proceeds. On the other hand, the low-velocity gas remains close to the protostar because their initial velocity (i.e., launching velocity) is low. As a result, the inclined spine-like structure appears in the PV diagram. In addition, the ejected gas interacts with the envelope material and thus decelerates. Therefore, a sharp spine widens with time and the propagation or jet velocity decelerates to have a velocity plateau, as seen in the right panels of Figure 8.

In Figure 8(f), we can see a new spine labeled as [3] in the PV diagram which corresponds to a new mass ejection event occurring at $201.69$ yr after the protostar formation (Figure 8e). The spline inclines as time in the PV diagram. Reflecting the past mass ejection event, we can see, at least, three spike-like structures [1]-[3] in Figure 8(f). The spine gradually inclines and widens with time. In other words, the inclination angle and width of the spine at an old mass ejection event are more significant than that at a new mass ejection event.

Some morphological similarities between the observation (Figure 6) and simulation (Figure 8 left panels) are apparent. In addition, the PV diagrams made by the observation (Figure 5) and simulation (Figure 8) are qualitatively the same. However, the velocity distribution and the interval of the mass ejection events are not quantitatively the same between the observation and the simulation. Since the spatial resolution of the simulation is restricted, the timescale over which the PV diagram varies is extremely short (months). We thus can not quantitatively compare the observation with the simulation.

Following the simulation, we interpret the observed PV diagram velocity distribution. Multiple mass ejection events produce the multiple spine-like structures seen in Figure 5. The different angles seen in the PV diagram can be explained by the time evolution of each jet component. The high-velocity gas component spreads fastest as the ejected gas moves away from the protostar. Thus, in the PV diagram the spine appears to incline with time. In addition, the interaction between the ejected gas and envelope material decelerates the jet and produces a velocity plateau. In Figure 5, for component [1] the farthest part of the spine reaches a constant velocity (i.e., velocity plateau), potential evidence for the decelerated gas. Such a velocity plateau has been reported in other protostellar systems (e.g., Wang et al. 2014; Matsushita et al. 2019; Morii et al. 2021b).

The PV diagram obtained from the simulation shows similar characteristics to those obtained from the SiO observation. However, the simulation PV diagram seems to differ from that found by previous outflow studies using CO, which mainly traces an older and static outflow, such as reported by Santiago-García et al. (2009) and Plunkett et al. (2015). Hence, we consider that the chemical processes in this extremely young jet are very active and may not be in a steady state. In this subsection, we discuss possible origins of the SiO emission observed toward MMS 1. We also discuss how they possibly relate to the mass ejection phenomena produced by the dynamical model presented in Subsection 4.1.

Three different paths are considered to release silicates into the gas phase and subsequently form the SiO molecules observed in jets (Cabrit et al., 2012; Podio et al., 2021). The first path is the sputtering of silicate grains in C-shocks, which requires a velocity differential in the range of 10 and 40 km s^-1 to produce SiO in the gas phase (e.g., Schilke et al. 1997; Caselli et al. 1997; Panoglou et al. 2012). The second path is the release of silicon into the gas phase by vaporization in grain-grain collisions, which is expected to occur in slow J-shocks ($\lesssim$50 km s^-1) with the amount of silicon released into the gas phase on the order of a few percent (Guillet et al., 2009). The third path is the release of silicon into the gas phase in the region where the gas temperature reaches the dust sublimation temperature (Glassgold et al., 1991; Tabone et al., 2020).

The first case (C-shocks) has been discussed mainly for forming SiO in previously observed sources (Schilke et al., 1997; Gusdorf et al., 2008a, b) for which SiO is known as a shock tracer. The second case (J-shocks) can be expected to occur in specific locations, such as at the apex of bow shocks (e.g., Hartigan et al. 2004). The PV-diagrams in the right panels of Figure 8 indicate that the jet velocity is consistent with that required to release the silicate grains into the gas phase either through the C- or J-shocks. Moreover, the numerical simulation performed by Machida et al. (2020) reveals that the region close to the jet driving region has density and temperature discontinuities, indicating the presence of J-shocks. These considerations imply that the observed SiO emission might be produced by both the first and second pathways above, C-shocks or J-shocks.

Regarding the third pathway, the launching point of the high-speed jet ($v_{\rm{jet}}{\sim}$100 km s^-1) is expected to be several $\times$ 0.1 au (Figure 8 right panels). The temperature expected at these radii is a few $\times$ 1000 K, assuming a source bolometric luminosity of 1 L_⊙. This temperature is comparable to the dust sublimation temperature of 1500 K (e.g., Vaidya et al. 2009), implying that some of the observed high-velocity SiO components might be formed by this additional route.

Our gas-dynamic model (left panels of Figure 8) does not consider chemical networks. However, the physical conditions expected from this model are consistent with the conditions under which the SiO molecule could be formed, as predicted by the previous theoretical studies described above. Chemohydrodynamic analysis, where the chemical distribution is affected by the gas-dust dynamical evolution (e.g., Flower & Pineau des Forêts 2003; Godard et al. 2019; Castellanos-Ramírez et al. 2018), will be required to make complete comparisons in future studies.

As discussed in Section 4.1, both observations and simulations clearly show evidence of the non-steady mass ejection events within protostellar jets. Recent ALMA CO and SiO observations toward protostellar sources reveal chains of knots within protostellar jets, suggesting episodic mass ejection events (Plunkett et al. 2015; Matsushita et al. 2019; Jhan et al. 2022; Dutta et al. 2022, 2023, and this work). The dynamical timescales between these ejection events range from several years to thousands of years. At the same time, recent flux variability surveys performed at sub-millimeter and mid-infrared wavelengths suggest that the fraction of variable sources is high during this protostellar phase (Johnstone et al. 2018; Lee et al. 2021; Park et al. 2021). Two distinct forms of variability are observed: years-long secular variability and short-timescale stochastic variability. The range of timescales is attributed to a variety of processes including underlying changes in the disk accretion rate, geometric changes in the circumstellar disk, hydromagnetic interactions between the stellar surface and inner disk, magnetic reconnection of the stellar magnetosphere, and the number or size of spots on the stellar surface (Park et al., 2021; Fischer et al., 2023). In conclusion, these monitoring studies demonstrate the importance of understanding the non-steady nature the star formation activities.

Figure 9 presents 850 $\mu$m continuum flux monitoring results for MMS 1 obtained by the JCMT Transient Survey (Herczeg et al., 2017). The observations cover late 2015 through mid-2023, almost eight years of roughly monthly cadence. For bright sources, such as MMS 1, the relative flux calibration uncertainty between epochs is better than 2% (Mairs et al., 2017, Mairs et al. 2024 accepted by ApJ). No drastic flux change is detected from the 850 $\mu$m light curve, suggesting that there has not been a clear burst event within the last 8 years. There is, however, a distinct $\sim$ 0.5% per year gradual decline in the brightness - therefore a candidate linear variable, according to Mairs et al. (2024 accepted). The decline is similar, though shallower, to the observed decline in HOPS 383 (Lee et al., 2021), a protostar that had an observed mid-IR burst between 2004 and 2006 (Safron et al., 2015), about ten years prior to the start of the JCMT survey.

For MMS 1, the timescale of the episodic mass ejection events associated with the SiO knots presented in Figure 6 ranges from 5 to 50 years. Thus, the monitored time period might not be long enough to detect the associated episodic mass accretion phenomena. We note that 850 $\mu$m burst events have been observed by the JCMT Transient Survey for a few protostars allowing for detailed consideration of the events. HOPS 373 underwent a months-long burst that was analysed by Yoon et al. (2022) and EC53 (also known as V371 Ser) has continuing quasi-periodic variations observed at both near-IR (Hodapp et al., 2012) and sub-mm (Yoo et al., 2017) leading to a multi-wavelength analysis by Lee et al. (2020) and an interferometric analysis by Francis et al. (2022). Thus far, however, neither of those two sources have had features in their outflow connected to the observed accretion bursts. Circumstantial evidence for accretion variability correlating with episodic mass ejection has been presented for a few Orion outflow sources through comparison of the ejection timescales and the observed accretion variability timescale found from mid-IR or sub-mm monitoring (Jhan et al., 2022; Dutta et al., 2023). Continued monitoring of MMS 1 at 850 $\mu$m, combined with proper motion measurements of the SiO knots, will enable us to better determine the connection between episodic mass accretion events and episodic mass ejection.

Considering the current ALMA angular resolution of ${\sim}0{\arcsec}.2$ and the observed maximum jet velocity of ${\sim}$70 km s^-1, our observations are only sensitive to the timescale longer than $\sim$5 years. In order to detect less significant episodic events associated with a shorter time period than five years, higher angular resolution observations (${\lesssim}0{\arcsec}.2$) and mild shock tracers such as CH₃OH and H₂CO or S-bearing molecules might be useful (e.g., Lee et al. 2018; Tychoniec et al. 2019, 2021; Codella et al. 2020, 2021).