The CO Outflow Components Ejected by a Recent Accretion Event in B335

Protostellar outflows often present a knotty appearance, providing evidence of sporadic accretion in stellar mass growth. To understand the direct relation between mass accretion and ejection, we analyze the contemporaneous accretion activity and associated ejection components in B335. B335 has brightened in the mid-IR by 2.5 mag since 2010, indicating increased luminosity, presumably due to an increased mass accretion rate onto the protostar. Atacama Large Millimeter/submillimeter Array (ALMA) observations of 12CO emission in the outflow reveal high-velocity emission, estimated to have been ejected 4.6–2 yr before the ALMA observations and consistent with the jump in mid-IR brightness. The consistency in timing suggests that the detected high-velocity ejection components are directly linked to the most recent accretion activity. We calculated the kinetic energy, momentum, and force for the ejection component associated with the most recent accretion activity and found that, at least, about 1.0% of the accreted mass has been ejected. More accurate information on the jet inclination and the temperature of the ejected gas components will better constrain the ejected mass induced by the recently enhanced accretion event.


INTRODUCTION
Episodic accretion plays an important role in stellar mass assembly, including as one of the solutions to the luminosity problem (e.g., Offner & McKee 2011;Dunham & Vorobyov 2012;Fischer et al. 2017).Theoretical models and observational evidence both suggest that episodic accretion is more extreme and frequent in the Class0/I stages (e.g., Vorobyov & Basu 2015;Park et al. 2021;Zakri et al. 2022).Although a few protostellar outbursts have been observed via brightness variations that last for years or even decades (Kóspál et al. 2007;Safron et al. 2015;Park et al. 2021), most observational evidence for episodic accretion comes from knotty features in the outflow/jet (Hirano et al. 2006;Lee et al. 2009Lee et al. , 2010;;Plunkett et al. 2015), or chemical signatures Corresponding author: Jeong-Eun Lee lee.jeongeun@snu.ac.kr in envelopes (e.g., Lee 2007;Kim et al. 2012;Jørgensen et al. 2015).
Outflows have been ubiquitously found around young stellar objects (YSOs) (e.g., Hatchell & Dunham 2009) since the first detection 40 years ago (Snell et al. 1980).Two types of outflows are commonly observed: (1) wideopening slow (few to 10 km s −1 ) molecular outflows, typically observed with molecular lines such as CO, and (2) well-collimated fast (100 km s −1 ) outflows/jets, typically observed in the near-infrared with forbidden iron emission [Fe II] and H 2 or seen in shock-tracing molecular species such as CO and SiO (Bally 2016;Bjerkeli et al. 2019;Dutta et al. 2022a;Lee 2020).
In addition, well-collimated outflows/jets consist of gas launched from the YSO and often present a knotty appearance (e.g.Hirano et al. 2006;Lee et al. 2009Lee et al. , 2010;;Plunkett et al. 2015), indicative of variable, potentially discrete, ejection events presumably linked to an unsteady and possibly eruptive mass accretion process.Dutta et al. (2022a,b) reported the knot features of SiO and CO emission observed in the protostellar stages from early to evolved and calculated the jet mass-loss rate.Jhan et al. (2022) also showed the presence of knotty features of SiO emission in six Class 0 and I sources.
The accretion-related YSO variability is another observational evidence of episodic accretion.While historical searches for outbursts focused on optical wavelengths, YSOs at the Class 0/I stages are deeply embedded in the surrounding envelope, thus YSOs at these stages cannot be observed in the optical/nearinfrared.However, sustained mid-infrared monitoring by WISE/NEOWISE (Contreras Peña et al. 2020;Park et al. 2021) and sub-mm monitoring at JCMT (Lee et al. 2021) revealed the variability of YSOs at these stages.Park et al. (2021) performed an ensemble study to classify YSO variability using 6.5 years (14 epochs) of NE-OWISE light curves for 717 Class 0/I protostars, identified by previous Spitzer and Herschel surveys (Megeath et al. 2012;Dunham et al. 2015;Esplin & Luhman 2019).They found that ∼ 55% of Class 0/I protostars can be classified as variable stars, with the majority undergoing small amplitude -order unity -variations and a smaller fraction having much larger brightness changes.Furthermore, the mid-IR variability of protostars has been shown to be typically a reliable tracer of accretion luminosity changes in these systems (Contreras Peña et al. 2020).
Then, the protostellar brightness variability should be causally connected with the knotty outflow/jet features; the mass flow onto the forming star is intimately linked to accretion through the disk and angular momentum dispersal from the disk via an outflow in an unstable process.However, only tentative connections have been identified between an increase in protostellar outflow activity and a known luminosity burst (e.g., Ellerbroek et al. 2014).Definitive measurements of an outflow event coupled to a specific accretion burst would demonstrate an incontrovertible physical link.In this paper, we present a unique case study of B335, where the contemporaneous accretion and ejection events were revealed.
B335 is an isolated dark globule located at a distance of 164.5 pc (Watson 2020).Previous studies (Yen et al. 2010;Evans et al. 2015;Yen et al. 2015a,b) revealed infalling motions toward the center of B335, while no clear rotation motions are observed at large scales.The lack of clear rotation motion at 1000 AU implies that the size of the protoplanetary Keplerian disk is <16 AU if it exists (using the updated distance of 164.5 pc; Yen et al. 2010Yen et al. , 2015a,b),b).Imai et al. (2016Imai et al. ( , 2019) ) have reported a number of complex organic molecules (COMs) from B335, suggesting that B335 harbors a hot corino associated with the central Class 0 protostar IRAS 19345+0727.
The large-extended outflow of 12 CO J = 1 − 0 for B335 was revealed by Hirano et al. (1988), nearly 2 ′ (north to south) × 8 ′ (east to west).Yen et al. (2010) reported the high-velocity 12 CO J = 2−1 component in B335 that have size of 6 ′′ (north to south) × 10 ′′ (east to west).Recently, Bjerkeli et al. (2019) used 12 CO J = 2 − 1 of combined high spatial resolution ALMA data to uncover a hint of a recent accretion outburst as a high-velocity ejection component in the B335.However, direct evidence of the accretion event responsible for the ejection activity was absent from that analysis.Therefore, here we analyze the WISE/NEOWISE light curve of B335 to study its time-dependent variability and connect it with the fossil record of the time-dependent evolution of the outflow/jet observed with ALMA for the same high-velocity ejection component.
In this work, we aim to investigate the relationship between contemporaneous accretion and ejection activities in B335.Matching the properties of the ejecta with the recent accretion history greatly aids the determination of the mass accreted onto the protostar versus the mass dispersed away from the system.We present the observations and data reduction in Section 2. Our results are described in Section 3. We discuss the physical relationship between the derived contemporaneous mass accretion rate and mass ejection rate in Section 4. The source of uncertainty in ejection mass estimation is discussed in Section 5. Finally, conclusions are given in Section 6.

ALMA
The ALMA data archive contains several observational datasets of B3351 .For the analysis presented in this work we obtained the highest-angular resolution data available in the archive.A high-resolution observation was carried out with ALMA in band 6 between October 21 and 29, 2017 (2017.1.00288.S PI: Bjerkeli, Per).The data contained five execution blocks and cover five spectral windows (SPWs).The SPWs included the 12 CO J = 2 − 1 transition, 13 CO J = 2 − 1 transition, C 18 O J = 2 − 1 transition, SiO J = 5 − 4 transition, and continuum emission, respectively.In particular, the spectral resolution of the data, including 12 CO v = 0 J = 2 − 1 (230.538GHz), was set to 122.070 kHz (0.16 km s −1 , ∆v) while the total bandwidth was 117.187MHz.The continuum SPW was set to the spectral res- a The RMS noise level of continuum image was estimated from the emission-free region.
b The RMS noise level of the 12 CO maps was calculated per channel (channel width = 0.16 km s −1 ) and averaged over the line-free channels.
olution of 15.625 MHz and a total bandwidth of 2.000 GHz.The data were initially calibrated using CASA v5.1.1 (McMullin et al. 2007).
To obtain a better signal-to-noise ratio (SNR), we performed self-calibration for the continuum emission using CASA v6.2.1 (McMullin et al. 2007).The selfcalibrated continuum was imaged with natural weighting using the same CASA version of the CLEAN algorithm.The data cube with 12 CO had the continuum emission self-calibration solution applied and was then imaged from the continuum subtracted visibility data using tclean task at CASA v6.2.1.For the CO image, we tested three different weightings: uniform, natural, and Briggs with a robust parameter of 0.5.The highvelocity ejection emission was not detected in the uniform weighting image, and the SNR of the emission in the Briggs weighting image was too low.As a result, we adopted a natural weighting and a uvtaper of 0.06 ′′ to improve the SNR.

WISE/NEOWISE
The Wide-field Infrared Survey Explorer (WISE ) is a mid-IR survey of the entire sky using four bands, W1 (3.4 µm), W2 (4.6 µm), W3 (12 µm), and W4 (22 µm), with an angular resolution of 6.1 ′′ , 6.4 ′′ , 6.5 ′′ , and 12.0 ′′ , respectively, from January to September 2010 (Wright et al. 2010).The WISE mission continued to operate for an additional four months using W1 and W2 bands after the depletion of hydrogen for the cryostat of WISE.This was known as the NEOWISE Post-Cryogenic Mission (Mainzer et al. 2011).In September 2013, WISE was reactivated as the NEOWISE mission and has been performing observations in the W1 and W2 bands with a cadence of six months since 2013 (Mainzer et al. 2014).
In each visit to an area of the sky, WISE performed 10-20 observations over a period of 1-3 days.
B335 has been observed by the WISE /NEOWISE survey since 2010: in April and October 2010 and every six months since 2014.To construct the source light curve, single-epoch catalogs found at the NASA/IPAC Infrared Science Archive were queried using a 5 ′′ radius around 53000 54000 55000 56000 57000 58000 59000 60000 MJD (days)  We note that the mark for the velocity resolution is too small to be visible.
Figure 1 presents the 4.6 µm mid-IR light curve of B335.The data includes the WISE /NEOWISE W2 data described in Section 2.2, as well as Spitzer IRAC I2 band photometry.The latter value is derived by performing aperture photometry on mosaic images obtained from the NASA/IPAC Infrared Science Archive (IRSA).The aperture size was set at 12 ′′ following the guidelines from the IRAC instrument handbook.Finally, the IRAC I2 photometric point was corrected to match the WISE W2 passband following the equations of Antoniucci et al. (2014).
The figure shows that B335 brightened by 2.5 magnitudes between WISE and NEOWISE observations, reaching a peak brightness of W2= 9.2 mag.Due to the gap between WISE (55304 MJD) and the start of NEOWISE (56949 MJD) observations, it is difficult to interpret the shape of the outburst.The light curve appears to show that YSO brightened smoothly between 2010 and 2017-2018 when it reached its peak brightness, but a quiescent phase followed by a sudden outburst around 2013-2014 could also describe the shape of the light curve.The additional Spitzer point also adds some uncertainty in interpreting the light curve.The brighter photometric point could be interpreted as B335 going through periodic outbursts.Such behavior, which could result from dynamical perturbations from stellar or planetary companions, has been observed in an increasing number of YSOs (Hodapp & Chini 2015;Dahm & Hillenbrand 2020;Guo et al. 2022).
In spite of the difficulties in interpreting the shape of the outburst, the high-amplitude variability observed between WISE and NEOWISE observations is likely associated with a recent accretion outburst of B335.The ALMA observation was taken in October 2017 when the outburst in the YSO had already reached its maximum brightness.

The high-velocity 12 CO clumps ejected by the recent outburst
B335 is estimated to have a systemic velocity of 8.3-8.5 km s −1 (v sys = 8.3 km s −1 adopted here, Evans et al. 2005;Jørgensen et al. 2007;Yen et al. 2011;Mottram et al. 2014).The systemic velocity was taken into account for all velocities used in this paper.The outflow inclination is ∼10 • from the plane of the sky (Hirano et al. 1988).Stutz et al. (2008) reported the inclination of the disk of 87 • from the line of sight (nearly edge-on), corresponding to the outflow inclination of 3 • from the plane of the sky.We adopt 10 • as the standard to derive outflow properties but explore the effect of the smaller inclination in Section 5.3.−27.3 km s −1 .Although the SiO and 13 CO emission were not detected at the high-velocity, the same CO high-velocity component is also found in the analysis by Bjerkeli et al. (2019).Figure 3 shows the peak intensity map and the velocity map at the peak intensity, where the considered velocity range is −37.3 km s −1 to −17.3 km s −1 , which fully covers the high-velocity ejection component.It is also revealed that this highvelocity ejection component consists of multiple components in Figure 3.
Outflows commonly show a bipolar structure.We thus examined whether red-shifted high-velocity ejection components associated with the blue-shifted highvelocity ejection components could be found and confirmed that no red-shifted high-velocity ejection components were detected in the ALMA observations.
To calculate the time when the high-velocity components were ejected, their velocities and distances between these components and the protostar are needed, taking into account the inclination of the target.The earliest ejected high-velocity component has a velocity of 145.6 km s −1 , after adjusting for the 10 • inclination of the outflow.It is located at a distance of ∼144.4 au from the central source.The most recent ejected component has a velocity of 165.7 km s −1 and a distance of ∼71.0 au (see the red and blue arrows, respectively, in the right panel of Figure 3).We estimate the kinematic timescale of both components by dividing the distance by the velocity.Therefore, we infer that the high-velocity ejec-tion components were ejected 4.6 to 2.0 years before the ALMA observation date.The gray-shaded region in Figure 1 represents the timescale over which the highvelocity ejection clumps were ejected.This is consistent with the period of increase when B335 brightened significantly.
The outflow inclination of 10 • was provided by Hirano et al. (1988) without a measurement error, making it challenging to estimate the uncertainty in the ejection period of high-velocity clumps.For an estimation, we assume a 30% uncertainty in the outflow inclination, i.e., 10±3 • and adopt the velocities and distances for the earliest and latest high-velocity ejection components.The ejection period ranges from 6.0 to 2.6 yr before the ALMA observation date for the outflow inclination of 13 • and from 3.3 to 1.4 yr for the outflow inclination of 7 • .As a result, the uncertainty of the ejection period could be up to 3 years given the 30% uncertainty of the outflow inclination.

Properties of ejected gas components
We compute the properties (mass, M ejection ; kinematic energy, E ejection ; momentum, P ejection ; and force, F ejection ) of the high-velocity ejection clumps following Dunham et al. (2014).Since the ejection properties depend on the mass and the H 2 column density is required to estimate the mass of ejected gas, we first calculate N(H 2 ) v,pixel using the 12 CO data for each pixel and each velocity channel within a selected range, where we assume that the ejection components are optically thin and follow local thermal equilibrium (LTE) conditions.The column density of the H 2 then can be estimated using: where f (J, T ex , X mol ) is a function of the quantum number of the lower rotation state, J, the molecule abundance relative to H 2 , X mole , and the excitation temperature of the ejected material, T ex .T mb ∆v is the temperature integrated for each pixel and velocity.
In Equation 2, we adopt a standard CO abundance relative to H 2 , X CO = 1.7 × 10 −4 (Lacy et al. 2017).The rest frequency of 12 CO J = 2 − 1 corresponds to ν = 230.538GHz, the 12 CO dipole moment has a value of µ = 0.11 Debye, and the upper energy level, E J+1 , equals 2.291 × 10 −15 erg.The partition function, Q(T ex ), and the degeneracy of the lower state, g J , are expressed as ∑ ∞ J=0 g J e −E J /kT ex and 2J + 1, respectively.The excitation temperature, T ex , is an unknown parameter when calculating N(H 2 ).Previous studies have found evidence for warm gas in outflows (e.g.Hatchell et al. 1999b,a;Nisini et al. 2000;Giannini et al. 2001).For example, Hatchell et al. (2007) adopted a value of T ex = 50 K by considering the evidence of warm gas in outflow, while others (Green et al. 2013;Lee et al. 2014;Je et al. 2015) found a higher T ex , close to 350 K, for the warm CO gas using higher J levels.
To model outflows, Dunham et al. (2010) chose values of T ex ranging between 10 and 100 K.In particular, they used T ex = 17.6 K to calculate lower limits of the outflow properties.We follow a similar approach, finding the minimum value of f (J, T ex , X mol ) and the corresponding excitation temperature T ex for the 12 CO J = 2 − 1 transition line.Figure 4 shows the value of the f (J, T ex , X mol ) versus T ex for the 12 CO J = 2 − 1 transition line with the minimum value of the f (J, T ex , X mol ) when T ex is 17.5 K.The right vertical axis of Figure 4 shows the enhancement factor of ejection mass to the minimum mass as a function of T ex .For instance, the mass of ejected gas with T ex = 100 K is 2.5 times larger than the derived minimum mass with T ex = 17.5 K.
Since the high-velocity 12 CO components were induced by the recent accretion activity, the T ex could be higher than 17.5 K. Additionally, as mentioned above, the presence of warm gas in the outflow indicates that the value of T ex might also greatly exceed the adopted value.Nevertheless, We adopt T ex = 17.5 K to obtain the minimum value for the ejection mass.Normalized f(J, T ex , X mole ) Figure 4. Left y-axis: value of f (J, T ex , X mol ) as a function of T ex for the 12 CO J = 2 − 1 transition line.The red dashed line indicates T ex = 17.5 K where f (J, T ex , X mol ) has a minimum value, f min (J, T ex = 17.5 K, X mol ).Right y-axis: Nomarlized f (J, T ex , X mol ) against the f min (J, T ex = 17.5 K, X mol ) for the 12 CO J = 2 − 1 transition line.The gray dashed line repersents the value of f min (J, T ex = 17.5 K, X mol ).
The ejection mass for each pixel and each velocity channel is calculated as where µ H 2 is the mean molecular weight per hydrogen molecule, with µ H 2 = 2.809 (Evans et al. 2022), m H is the mass of a hydrogen atom, and A pixel is the area of a pixel.The total H 2 column density, N(H 2 ), and total mass of ejected gas, M ejection , are calculated by summing N(H 2 ) v,pixel and M v,pixel over the pixels and velocity channels within a selected region and velocity range.
To calculate the total column density and total ejection mass over the high-velocity ejection clumps, we select a 1 ′′ ×1 ′′ region corresponding to 1.0 ′′ to 0 ′′ for ∆α and -0.5 ′′ to 0.5 ′′ for ∆δ (the orange dashed boxes in Figure 3).We also adopt the velocity range from -37.3 km s −1 to -17.3 km s −1 .We use only the emission above 5σ within the selected region and velocity range to obtain N(H 2 ) v,pixel and M v,pixel .Finally, the total H 2 column density and the total mass of the high-velocity ejection components are 2.1 × 10 23 cm −2 and 7.5 × 10 −8 M ⊙ , respectively.
The kinetic energy, momentum, and force of ejected gas for each pixel and velocity channel are calculated following (Dunham et al. 2010): where τ v,pixel is the dynamical time for each pixel and each velocity channel.The value of τ v,pixel can be estimated by dividing the distance between the given position and the continuum peak position by the ejection velocity.
The total value of each ejection property (E ejection = Σ E v,pixel , P ejection = Σ P v,pixel , F ejection = Σ F v,pixel ) is derived by summing over the pixels and velocity channels.Following this approach, we estimated E ejection , P ejection , and F ejection to be 1.9 ×10 40 erg, 1.2 ×10 −5 M ⊙ km s −1 , and 3.8 ×10 −6 M ⊙ km s −1 yr −1 , respectively.
For the high-velocity 12 CO outflow with the propagation velocity of ∼160 km s −1 over the size of 1600 au by 1000 au, Yen et al. (2010) reported F outflow = 4.1 × 10 −5 M ⊙ km s −1 yr −1 when scaled by the distance of 164.5 pc.This is approximately one order of magnitude higher than our result due to the larger coverage of the outflow.On the other hand, Hirano et al. (1988) obtained F outflow = 1.8 × 10 −4 M ⊙ km s −1 yr −1 , when scaled by the distance of 164.5 pc, for the lowvelocity 12 CO outflow of B335.This low-velocity outflow is even more extended than the high-velocity 12 CO outflow, making the outflow force almost two orders of magnitude higher than our result.In addition to the different scales, F outflow reported by previous studies must have been accumulated over a long period of time, while the F ejection calculated in this study is a value corresponding to one specific accretion event.

Mass ratio between ejection and accretion
In the early embedded stages, the central luminosity is mainly produced by accretion.Thus, the mass accretion rate, Ṁacc , can be obtained from the following equation: where L acc is the accretion luminosity, R ⋆ is the radius of the protostar, M ⋆ is the mass of the protostar, G is the gravitational constant, and f acc is the fraction of the accretion energy that is radiated (Hartmann et al. 2011).
We first estimate L acc to obtain the mass accretion rate during the timescale over which the highvelocity ejection clumps were ejected.We convert the W2 magnitude into W2 flux density (see right y-axis in Figure 1) as a proxy for the central luminosity.
The relation between the central luminosity and W2 flux density of B335 has been explored based on the 3D continuum modeling by Evans et al. (2023); L/L ⊙ = [1.837+0.845(Sν (W 2)×0.67)+0.001(Sν (W 2)×0.67) 2 ].The maximum luminosity, i.e., the outburst luminosity, is predicted by the model to be a factor of 5-7 higher than the quiescent luminosity of 3 L ⊙ .Alternatively, Contreras Peña et al. ( 2020) found a relation between the flux density of the mid-infrared and central luminosity for embedded protostars, F IR ∝ L 1.5 .According to this relation, the outburst luminosity at MJD 58391 is enhanced by a factor of ∼5 over the quiescent luminosity at MJD 55304.This result is consistent with the luminosity calculated using the relation found in Evans et al. (2023).Therefore, we adopt the relation between W2 flux density and central luminosity of B335 by Evans et al. (2023) to calculate the variation of central luminosity over the recent burst.
To calculate the mass accretion rate, we need some additional properties of B335 (R ⋆ , M ⋆ , and f acc = 0.5).
In the 3D continuum modeling explored by Evans et al. (2023), the luminosity depends upon R ⋆ and the temperature, T ⋆ , of the source.Evans et al. (2023) held fixed T ⋆ and varied R ⋆ to match the observed quiescent luminosity, which is the SED modeling result using various photometric data observed before MJD 57000.The fixed value of T ⋆ is 7000 K matches well the near-infrared part of the SED.To match the quiescent luminosity, R ⋆ would be 1.17 R ⊙ .We note that the near-infrared part of the SED may be contaminated due to various causes such as the rotation rate, inclination angle, outflow cavity properties, disk properties, stellar properties, and absorption of various ice molecules (Yang et al. 2017;Chu et al. 2020;Evans et al. 2023).Thus, the fixed T ⋆ of 7000 K has a large uncertainty as does R ⋆ .Nevertheless, here we adopt an R ⋆ of 1.17 R ⊙ as the result of the best-fit model.
This same model provides an infall age, t col = 4 × 10 4 yr, and mass infall rate, 6.26 ×10 −6 M ⊙ yr −1 (Evans et al. 2023).Assuming all the mass that infalls from the envelope into the disk has fallen onto the star, M ⋆ is 0.25 M ⊙ .Other studies (Yen et al. 2015b;Imai et al. 2019), however, estimated the M ⋆ based on the PV diagrams.Yen et al. (2015b) derived the M ⋆ of 0.05 M ⊙ , while Imai et al. (2019) provided a range of mass as 0.02 -0.06 M ⊙ .Therefore, we adopt the M ⋆ to range between 0.02 and 0.25 M ⊙ .
Ṁacc is estimated using equation 5 and the values adopted above and presented in Figure 5.The average Ṁacc during the burst event (the gray shaded region, ∼2.6 yr) that launched the high-velocity ejection clumps is in the range of ∼3.1 ×10 −6 to ∼3.8 ×10 −5 M ⊙ yr −1 , depending on the adopted stellar mass.Therefore, the total mass accreted over the 2.6 yr, M acc becomes 8.0 × 10 −6 to 1.0 × 10 −4 M ⊙ .
The average mass ejection rate, Ṁejection , during the period when the high-velocity ejection components were ejected is estimated by dividing the total mass of the high-velocity ejection components by the timescale over which the ejection occurred, 2.6 yr, resulting in 2.9 × 10 −8 M ⊙ yr −1 .Therefore, at least, ∼0.07% to ∼0.9% of the accreted mass was ejected by the recent burst event.

Stellar properties
As discussed in Section 4.2, the estimation of Ṁacc depends on stellar properties such as L acc , R ⋆ , M ⋆ , and f acc .We thus require accurate stellar properties to obtain an accurate mass accretion rate.The stellar properties adopted in this study are, however, quite uncertain.
First, for L acc , the luminosity enhancement factor is uncertain at the 30% level, ranging between 5 and 7 times the quiescent luminosity in Evans et al. (2023).
The relation, L/L ⊙ = [1.837+0.845(Sν (W 2)×0.67)+0.001(Sν (W 2)×0.67) 2 ], suggests the luminosity is enhanced by almost a factor of 7, and we adopt this value for our calculation.However, if the luminosity increases by only a factor of 5 over the quiescent luminosity, and the other stellar properties are fixed, then the derived Ṁacc would decrease by about 30%.
Second, the best-fit model with T ⋆ = 7000 K has R ⋆ = 1.17 R ⊙ .However, the fixed T ⋆ of 7000 K is quite uncertain; thus, the 1.17 R ⊙ is poorly constrained.The typical protostellar radius is estimated as ∼3 R ⊙ (Stahler et al. 1980), which is 2.6 times larger.The T ⋆ is 4370 K when we adopt the R ⋆ = 3 R ⊙ .Using the typical protostellar radius and fixing all other parameters, the mass accretion rate would increase by a factor of 2.6.
Finally, f acc has a value as a f acc ≤ 1 (Hartmann et al. 2016).Previous studies (Evans et al. 2015;Yen et al. 2015b) have adopted a value of 0.5 and 1, respectively, for calculating the mass accretion rate of B335.This introduces an additional factor of 2 uncertainty in the mass accretion rate.
Therefore, the mass accretion rate varies with the adopted stellar properties, as discussed above.Ṁacc has a maximum value of ∼9.8 × 10 −5 M ⊙ yr −1 if the luminosity enhancement fator is 7, M ⋆ = 0.02 M ⊙ , and R ⋆ = 3 R ⊙ .On the contrary, Ṁacc has a minimum value of ∼2.2 × 10 −6 M ⊙ yr −1 if the luminosity increases by only a factor of 5, M ⋆ = 0.25 M ⊙ , and R ⋆ = 1.17 R ⊙ .The f acc is fixed at 0.5 in both cases.These extremes result in a factor of 45 difference in the derived mass accretion rate.Considering the uncertainties, the ejected mass ranges from 0.03% to 1.3% of the accreted mass during the recent burst event.

Excitation temperature
We adopt T ex =17.5 K to determine the lower limit of the ejection mass.It is not unreasonable to adopt this value since the low-J (J = 2 − 1) 12 CO line is dominated by the coldest gas.However, as mentioned above, the 12 CO J = 2 − 1 line used in this study was induced by the recent outburst, and therefore T ex could be higher than the value we adopted.
If T ex is ∼200 K, the total ejection mass is a factor of 4.6 greater than when T ex is ∼17.5 K (See right yaxis in Figure 4).The mass ratio of ejection to accretion also increases by the same factor.Thus, the ratio, M ejection /M acc ranges from 0.03% to 6.0%.To determine a more accurate mass ratio between accretion and ejection, other 12 CO transition lines, such as 12 CO J = 3 − 2, are required for an accurate calculation of T ex .

Inclination
The typical velocity of shock knots observed in outflows ranges 100 -150 km s −1 (Bachiller 1996;Dutta et al. 2020;Lee 2020;Jhan et al. 2022).The velocity calculated with the inclination of 10 • (Hirano et al. 1988) is consistent with the typical value.However, if we adopt the disk inclination of 87 • from the line of sight, as derived by Stutz et al. (2008), which corresponds to the outflow inclination of 3 • from the plane of the sky, the estimated velocity of knots reaches an extremely high velocity of ∼500 km s −1 .In particular, the velocity and distance of the earliest ejected component are 483.1 km s −1 and 142.4 au, while the velocity and distance of the latest ejected component are 549.8km s −1 and 70.1 au.If the velocity of the ejection components is related to the Keplerian velocity at the outflow launching radius, it can be calculated by r 0 = GM ⋆ /v kep 2 .Thus, the launching radius is ∼0.02R ⊙ for M ⋆ = 0.02 M ⊙ and ∼0.2 R ⊙ for M ⋆ = 0.25 M ⊙ .This radius is smaller than the stellar radius we adopted.Nonetheless, adopting these velocities and distances, the high-velocity ejection components were ejected 1.4 to 0.6 years before the ALMA observation date.The Ṁejection with outflow inclination of 3 • would be 9.5 ×10 −8 M ⊙ yr −1 , which is calculated by dividing the M ejection by the timescale, 0.8 yr.
The inclination of the high-velocity outflow/jet is not necessarily the same as that of the larger slow outflow component.Therefore, it is important to know the correct inclination of the high-velocity ejection component to reduce the uncertainties associated with its physical parameters.For the accurate measurement of the inclination, monitoring observations of the high-velocity ejection component is critical.

CONCLUSIONS
We have revealed the contemporaneous accretion and ejection events that occurred recently in B335, using both ALMA and WISE/NEOWISE data.From the PV diagram of the ALMA 12 CO J = 2 − 1 image, a group of isolated high-velocity ejection components was revealed at ∼ −27.3 km s −1 .The high-velocity CO gas components were likely ejected with the velocity of ∼156 km s −1 , about 4.6 to 2.0 years before the ALMA observation date.The WISE/NEOWISE light curve of B335 presents a brightening event by 2.5 mag at 4.6 µm over ∼8 years, with the maximum brightness in 2018.A recent accretion event could have caused both the brightness increase of B335 and the ejection event, as recognized by the observed high-velocity ejection components.
The ejection properties such as E ejection , P ejection , and F ejection for the high-velocity ejection components are estimated.To investigate the physical causality between contemporaneous ejection and accretion activities, we calculate the total mass of the high-velocity ejection components and the accreted mass for the same period.During the latest accretion burst, the mass ejection rate is 2.9 × 10 −8 M ⊙ yr −1 , while the mass accretion rate ranges from 3.1 × 10 −6 M ⊙ yr −1 to 3.8 × 10 −5 M ⊙ yr −1 , depending on the adopted physical parameters, in the timescale of about 2.6 years, suggesting that about 0.07 -0.9% of the accretion mass was ejected.
Ṁejection , Ṁacc , the ratio between ejected mass and accreted mass, and ejection properties have uncertainty due to the uncertain parameters, particularly the inclination, adopted for the calculations.Therefore, constraining ejection physical parameters accurately is critical to understanding the physical causality between contemporaneous ejection activities and accretion events.

Figure 1 .
Figure 1.Light curve of B335 in the IRAC I2 band of Spitzer, W2 band of WISE, NEOWISE -Post Cryogenic, and NEO-WISE survey.The short solid lines attached to the data points represent the magnitude uncertainties at given epochs.The black dashed vertical line indicates the ALMA observation date (October 2017).The gray-shaded region represents the period when the high-velocity ejection clumps identified in the PV diagram shown in Figure 2 were ejected (March 2013 -October 2015).

Figure 2 .
Figure 2. Left: Integrated intensity map of 12 CO over the velocity range (-49.6 km s −1 < v range < 41.7 km s −1 ).Black contours represent continuum emission.Contour levels are from 10, 20, 40, 80, 160, 320 σ cont (σ cont = 0.015 mJy beam −1 ).The synthesized beams for the 12 CO cube data and continuum data are presented at the lower left corner of the map as a black ellipse and gray ellipse, respectively.The red dashed line represents the direction of the cut for a PV diagram.The purple dashed circle represents the location of the high-velocity ejection component indicated in the right panel.Right: PV diagram for 12 CO along the direction of the red dashed line shown in the left panel.Black contours indicate signals higher than 3 σ (σ = 1.7 mJy beam −1 ≃ 5 K).The high-velocity ejection components are marked with a purple dashed circle.The FWHM of the synthesized beam (∼0.087 ′′ ) and the velocity resolution (0.16 km s −1 ) are indicated with the vertical and horizontal bars, respectively.We note that the mark for the velocity resolution is too small to be visible.

Figure 3 .
Figure 3. Left: Peak intensity (moment 8) map of the high-velocity 12 CO clumps.Right: Velocity of the peak intensity (moment 9) map of the high-velocity 12 CO clumps.The velocity range used for both maps is from −37.3 km s −1 to −17.3 km s −1 , and only emission above 5 σ was used.The long red and short blue arrows indicate the distances from the central continuum peak to the earliest and latest ejected clumps among the high-velocity 12 CO clumps.The black contours represent continuum emission.Contour levels are from 10, 20, 40, 80, 160, 320 σ cont .The synthesized beams for each moment map and the continuum data are presented at the lower left corner of the maps as black and gray ellipses, respectively.The orange dashed box indicates the region over which the total mass of the high-velocity ejection components is calculated.

Table 1 .
Summary of ALMA Band 6 data used in this study