Multiple Outflows around a Single Protostar IRAS 15398−3359

We present the results of our mosaic observations of a single Class 0 protostar IRAS 15398−3359 with the Atacama Compact Array in the CO J = 2–1 line. The new observations covering a ∼2′ square region reveal elongated redshifted and blueshifted components, which are located at distances of ∼30″–75″ on the northern and southern sides of the protostar, respectively, in addition to the previously observed primary and secondary outflows. These elongated components exhibit Hubble-law-like velocity structures, i.e., an increase of velocity with increasing distance from the protostar, suggesting that they comprise the third outflow associated with the protostar. Besides, a new redshifted component is detected at radii of ∼40″–75″ on the northwestern side of the protostar. This redshifted component also exhibits a Hubble-law-like velocity profile, which could be a counterpart of the secondary outflow mostly detected at blueshifted velocities in a previous study. The three outflows are all misaligned by ∼20°–90°, and the dynamical timescale of the primary outflow is shorter than those of the other outflows, approximately by an order of magnitude. These facts hint that the outflow launch direction has significantly changed with time. The outflow direction may change if the rotational axis and the magnetic field are misaligned or if the dense core is turbulent. We favor the second scenario as the origin of the multiple outflows in IRAS 15398−3359, based on a comparison between the observational results and numerical simulations.


INTRODUCTION
Low-mass star formation occurs in dense cores through gravitational collapse.Throughout the star formation process, circumstellar disks form around protostars as a consequence of conservation of the angular momentum inherent in the initial dense cores (Terebey et al. 1984).Revealing the disk formation process is essential for understanding the diversity of protoplanetary disks and planets forming within them.
Early theoretical works developed analytical models of the disk formation process considering the gravitational collapse of a rotating dense core (e.g., Terebey et al. 1984;Basu 1998).On the other hand, later ideal magnetohydrodynamics (MHD) simulations indicated that the presence of a magnetic field aligned with the rotational axis of the dense core prevents disks from forming, because the magnetic field extracts angular momentum too efficiently (Allen et al. 2003;Mellon & Li 2008).These results contradicted the fact that protostellar/protoplanetary disks are ubiquitous around young stellar objects (YSOs) (e.g., Williams & Cieza 2011;Andrews 2020;Tobin et al. 2020;Ohashi et al. 2023).This issue has been solved by more detailed numerical simulations.Turbulence in dense cores, and/or misalignment between the magnetic field and the rotational axis of the dense core can suppress the efficiency of magnetic braking and allow disks to form (e.g., Joos et al. 2012;Seifried et al. 2013;Matsumoto et al. 2017).Non-ideal MHD effects also enable sizable disks to form even in strongly magnetized dense cores, as the magnetic field can be partially decoupled from neutral gas (Machida & Matsumoto 2011;Tomida et al. 2015;Hennebelle et al. 2020).
The presence of both turbulence and misalignment between the rotational axis and the associated magnetic field in dense cores have been observed.Molecular line observations of dense cores have revealed complex velocity structures that are not explained by a simple rotational motion, suggesting the presence of turbulence in these dense cores (Caselli et al. 2002;Chen et al. 2019;Sai et al. Gaudel et al. 2020;Sai et al. 2023).Polarization observations in continuum emission have reported that magnetic fields are randomly aligned with outflows associated with protostars, which may represent the rotational axes of disks (Hull et al. 2014;Yen et al. 2021b).Velocity gradients and magnetic fields of dense cores, and disk orientations are all frequently found to be misaligned (Galametz et al. 2018(Galametz et al. , 2020;;Yen et al. 2021a;Gupta et al. 2022), and these misalignments might affect the efficiency of the angular momentum transfer from large to small scales (Galametz et al. 2020;Gupta et al. 2022).
Although dense cores may be commonly turbulent or associated with a misaligned magnetic field, the dynamical role of these conditions in the star formation process remains still unclear.Numerical simulations of disk formation in turbulent dense cores or those associated with misaligned magnetic fields predict variations in the direction of the disk normal over time, because the angular momentum vector of infalling material changes with time (Seifried et al. 2013;Matsumoto et al. 2017;Bate 2018;Hirano & Machida 2019;Machida et al. 2020).In these simulations, the outflow orientation is also expected to vary with time as consequence, given that the outflow is typically launched in the direction of the disk normal (Seifried et al. 2013;Machida et al. 2020).
IRAS 15398−3359 (hereafter IRAS 15398) is a Class 0 protostar in the Lupus I star forming cloud at a distance of 155 pc (Zucker et al. 2020;Santamaría-Miranda et al. 2021;Thieme et al. 2023).Its bolometric temperature and luminosity are 50 K and 1.4 L ⊙ , respectively (Ohashi et al. 2023).The protostar is embedded in an infalling and rotating envelope and associated with a primary bipolar outflow along the northeast-southwest direction (Oya et al. 2014;Bjerkeli et al. 2016;Yen et al. 2017).In addition to the primary outflow, a monopolar, blueshifted secondary outflow launched in a direction orthogonal to the primary outflow has been reported through molecular line observations with a high sensitivity using Atacama Large Millimeter/submillimeter Array (ALMA) (Okoda et al. 2021).Molecular line observations with ALMA have also suggested the presence of a Keplerian disk with a radius of ≳ 30 au (Okoda et al. 2018;Thieme et al. 2023).Previous works have inferred a very small protostellar mass of ≲ 0.01 M ⊙ based on observations at angular resolutions of ≲ 0. ′′ 2 (Yen et al. 2017;Okoda et al. 2018).Recent observations at an angular resolution of ∼ 0. ′′ 15 have revealed a possible differential rotation signature of a Keplerian disk, and constrained the central protostellar mass range of ∼ 0.02-0.1 M ⊙ through fitting to the differential rotation signature (Thieme et al. 2023).
In this paper, we report the discovery of a third outflow around IRAS 15398, which is misaligned with the known primary and secondary outflows, and discuss the possible origins of the misaligned, multiple outflows.The finding of the largely misaligned multiple outflows in IRAS 15398 hints the dynamical importance of turbulence in the dense core and/or misalignment between the magnetic field and the rotational axis of the dense core.The outline of this paper is as follows.Observations and data reduction are summarized in Section 2. Observational results are presented in Section 3. Analysis on velocity structures is provided in Section 4, and possible implications are discussed in Section 5. Main results and conclusions are summarized in Section 6.

ACA 7-M ARRAY OBSERVATIONS
We have conducted mosaic observations of IRAS 15398 using the 7-m array of Atacama Compact Array (ACA) on 2019 December 19 during ALMA Cycle 7 with 10 antennae in Band 6.In this paper, we present the CO J =2-1 (230.538000GHz) results of the observations.The mosaic observations with 28 pointings covered ∼ 2 ′ square region centered at the protostellar position.The projected baseline lengths ranged from 6.5 to 37 kλ at 230 GHz.The observations were carried out in the frequency division mode (FDM).A spectral window with a bandwidth of 62.5 MHz and a spectral resolution of 61.0 kHz was assigned to the CO J =2-1 line, which provided a velocity resolution of ∼0.080 km s −1 .J1534−3526 and J1337−1257 were observed for the phase calibration and for the bandpass and flux calibrations, respectively.The data were calibrated with the pipeline version 42866M (Pipeline-CASA56-P1-B) of the Common Astronomy Software Applications package (CASA;McMullin et al. 2007) version 5.6.1.
We produced images from the visibility data using CASA 5.6.1 with the tclean task.We adopted the Briggs style weighting with the robust parameter of 0.5 and the velocity resolution of 0.08 km s −1 .Multiscale clean was used with the scale parameters of zero and 10, which correspond to the point source and about one beam size, respectively.The primary beam correction was applied to the cleaned image.We produced a final image cube with a velocity resolution of 0.16 km s −1 for analysis by binning two velocity channels.The angular resolution and root-mean-square (rms) noise level are 7. ′′ 7×4.′′ 1 (87 • ) and 0.13 Jy beam −1 , respectively.We adopted the source position of (α J2000 , δ J2000 )=(15:43:02.24,−34:09:06.8)for the figures and analyses in this paper (Yen et al. 2017).Moment 0 and 1 maps of the CO J = 2-1 emission are presented in Figure 1.The emission reveals several elongated structures at various spatial scales.Within a radius of 30 ′′ around the protostellar position, we see an elongated structure from southwest to northeast in a length of ∼30 ′′ .The elongated emission shows a velocity gradient in the same direction: a blueshifted component on the southwestern side and a redshifted component on the northeastern side of the protostar.Both blueshifted and redshifted components have a relatively high velocity of ±3 km s −1 with respect to the systemic velocity of 5.4 km s −1 (Thieme et al. 2023) in the moment 1 map.This elongated structure likely traces the primary outflow that has been reported in previous observations (Oya et al. 2014;Yen et al. 2017;Okoda et al. 2021) based on its morphology and velocity structure.Another component extends from the protostellar position to the southeast.This component has a LSR velocity of ∼ 4-5 km s −1 .Its elongation direction is close to those of the major axis of the disk and envelope (Yen et al. 2017;Okoda et al. 2018;Thieme et al. 2023) and the monopolar, blueshifted secondary outflow (P.A.∼ 140 • ; Okoda et al. 2021), although its size is larger than the size of the secondary outflow.At distances from the pro-tostar larger than 30 ′′ , blueshifted and redshifted components are observed on the southern and northern sides of the protostar, respectively.These components appear symmetric with respect to the protostellar position.They also exhibit velocity gradients in the north-south direction with relative velocities of ∼1-2 km s −1 with respect to the systemic velocity.These features suggest that they are an bipolar outflow associated with IRAS 15398.Another elongated, redshifted emission is found at (∆α, ∆δ) ∼ (−30 ′′ , 50 ′′ ), which is approximately aligned with a line passing through the protostellar position and the blueshifted component elongated in the southeast-northwest direction.Another elongated component is located at ∼40 ′′ west from the protostellar position.
In order to probe further details of the velocity structures of the CO emission, moment 0 maps of different velocity components are presented in Figure 2 (see also velocity channel maps shown in Figure A1).The highvelocity components shown in Figure 2(a) correspond to the emission elongated from the southwest to northeast near the protostar in Figure 1.These components are detected over a wide velocity range: the high-velocity blueshifted and redshifted components range from −6.0  to 1.6 km s −1 and from 8.3 to 15.0 km s −1 , respectively.More components appear at lower velocities, as shown in Figure 2(b).The northern redshifted (NR1 and NR2) components and the southern blueshifted (SB) component have velocities of ∼ 5.9-8.1 km s −1 and ∼ 1.7-4.9km s −1 , respectively.The blueshifted elongation from the protostellar position toward the southeast (SEB) and the blueshifted filamentary structure on the western side of the protostar (WB) have LSR velocities of ∼ 3.5-4.8km s −1 , which is closer to the systemic velocity of 5.4 km s −1 than the velocity of the SB component (see Figure A1).The redshifted filamentary structure on the western side (WR) also has a smaller velocity of ∼ 6.4-7.1 km s −1 than the the NR1 and NR2 components.

Velocity Structure
The velocity structures of low-velocity components are investigated in more detail in this section.Figure 3 presents position-velocity (PV) diagrams cut toward the NR1 (P.A.= 355 • ) and SB (P.A.= 182 • ) components.The position angle of each component is measured as an angle from north to the line passing through the protostellar position and the peak position of the emission, which is determined by a Gaussian fit.No clear emission is found around the systemic velocity of 5.4 km s −1 .This would be because the optically thick, extended fore-ground CO emission is resolved out.Two distinct components are found in each PV diagram.Within an offset of 10 ′′ , the emission extends over velocity ranges of 0-5 km s −1 and 6-10 km s −1 , which is the high-velocity components of the primary outflow.The NR1 and SB components are found at offsets of 20-75 ′′ .They exhibit the so-called Hubble-law like velocity profiles, i.e., an increase in the relative velocity with respect to the systemic velocity with increasing distance from the protostar.The velocity ranges of the NR1 and SB components are ∼ 1-3.5 km s −1 relative to the systemic velocity.These velocity structures and symmetric geometry with respect to the protostellar position suggest that the NR1 and SB components trace an outflow launched from IRAS 15398.This third outflow is highly misaligned with the primary and secondary outflows, as discussed in more detail later.
A PV diagram cut toward the NR2 component (P.A.= 332 • ) is presented in Figure 4.The emission of the NR2 component is located at offsets of 40 ′′ -75 ′′ , while the emission within an offset of 10 ′′ is associated with the primary outflow.The velocity structure of the NR2 component is similar to that of the NR1 component; the velocity increases from ∼ 6 to 7.5 km s −1 as the offset increases.The Hubble-law like velocity profile suggests that the NR2 component is also part of an outflow.Although its symmetric counterpart is not found, the NR2 and the SEB components are almost aligned in the same  offsets of ∼ 10 ′′ -20 ′′ and exhibits an almost constant velocity of ∼ 4.3 km s −1 .Its velocity structures are not well resolved with our limited spatial resolution.The other components at further larger offsets of ∼25 ′′ -70 ′′ at velocities of ∼4.5 and 5 km s −1 are likely part of ambient gas, as the CO emission is highly complex or almost resolved out at these velocities (see Figure A1).The SEB component is approximately along the disk/envelope major axis, as well as the secondary outflow detected on a smaller scale in previous observations (Yen et al. 2017;Okoda et al. 2018Okoda et al. , 2021;;Thieme et al. 2023).The previous studies have indicated that the infalling and rotating envelope exhibits redshifted velocity on the southeastern side, which does not match the blueshifted velocity of the SEB component.Thus, the SEB component is likely unrelated to the envelope.On the other hand, the velocity of the SEB component is similar to those of shock remnant cause by the secondary outflow (Okoda et al. 2021).The SEB component and the secondary outflow are compared in more detail in Section 5.1.
The WB and WR filamentary structures do not appear to be pointing towards the protostellar position, which suggests that they may not be part of the outflows launched from IRAS 15398.The PV diagrams cut along the WB and WR filamentary structures are presented in Figure 6.Their centers and position angles are measured though Gaussian fitting and are presented in Figure 6(a).Figure 6(b) and (c) show that these structures have almost constant LSR velocities around 4.3 and 6.8 km s −1 , respectively, along the elongation.No emission is detected from ∼ 4.5 to 6.5 km s −1 in the PV diagrams, which would be due to the resolved-out, optically thick foreground gas.The constant velocities close to the cloud velocity suggest that the WB and WR filamentary structures are part of the ambient cloud.Similar filamentary structures tracing large-scale ambient material have also been reported in other protostellar systems (Dutta et al. 2022).

Physical Properties of the Outflows
The NR1 and SB components are likely a pair of lobes of the third outflow because of their symmetric morphology and velocity structures with the Hubble-law like velocity profiles.The NR2 component shows the Hubblelaw like velocity profile and thus would also trace an outflow lobe.To characterize their geometry and velocity structures, we perform fitting with a wind-driven shell model (Shu et al. 1991;Lee et al. 2000).In the wind-driven shell model, the morphology and velocity structure of the outflow cavity are expressed as follows: where R is the cylindrical radius in a direction perpendicular to the outflow axis, z is the distance from the protostar along the outflow axis, and v R and v z are the expansion velocities.The constants R 0 and z 0 are fixed at 1 ′′ , and the unit of C is arcsec −1 .The outflow cavity is then projected onto the plane of the sky with the inclination angle i.Here, we define the inclination angle as an angle between the line-of-sight and the outflow axis; in our definition, i = 90 • corresponds to the edge-on disk configuration.We search for the parameter set of (C, v 0 , i) that best explains the observations.The fitting is performed in two steps.First, we fit the outflow shell morphology to constrain C sin i using the moment 0 maps presented in Figure 2.Then, adopting the constrained C sin i, we fit the velocity structure using the PV diagrams.In both steps, 5σ contours of the emission are used to define morphologies and velocity structures of the observed outflows.Details of the fitting process are presented in Appendix B.
The fitting results are presented in Figures 3, 4 and B1 and summarized in Table 1.The inclination angles of the three lobes are all ∼ 70 • .While the expansion velocities (v 0 ) of the SB and NR1 lobes are very similar, the expansion velocity of the NR2 lobe is smaller, approximately by a factor of two, than those of the SB and NR1 lobes.This suggests that the NR2 lobe is tracing a different ejection event than the third outflow.
We also estimate the dynamical timescales of these lobes to probe their relationship with the known primary and secondary outflows.We use their lengths and the maximum velocities to calculate the dynamical timescales, following the method used in previous works (Yıldız et al. 2015;Okoda et al. 2021).The lengths of the outflows, their maximum velocities, and the derived dynamical timescales are summarized in Table 2.Note that all these lobes are extended to the map edge, and thus, the outflow lengths represent lower limits.Consequently, the calculated dynamical timescales are also lower limits.The dynamical timescales of the NR1, NR2  and SB lobes are calculated to be ≳ 0.49-1.1 × 10 4 yr with correction of the inclination angle.A previous work estimated the dynamical timescale of the primary outflow to be ∼ 5 × 10 2 yr with correction of the inclination angle in the same manner (Bjerkeli et al. 2016), suggesting that the SB, NR1 and NR2 lobes are launched first, and then followed by the primary outflow.

Configurations of the Multiple Outflows
The SEB lobe is well aligned with the secondary outflow on a smaller scale of ∼ 200-1200 au (Okoda et al. 2021).The previous work has also revealed an arc-like structure in SO and SiO lines at the tip of the secondary outflow at a distance of 1200 au from the protostar.Based on the relatively narrow velocity range of the detected SO and SiO emission from ∼ 4 to 5.5 km s −1 , the authors suggested that the arc-like structure traces an old shock caused by the secondary outflow.The velocity range of the arc-like structure is almost consistent with the velocity range of the SEB lobe.The lack of the emission of the SEB lobe at velocities of 5-5.5 km s −1 would be due to the resolved-out, optically thick foreground gas (Figure A1).Faint, blueshifted CO emission was also detected near the SO/SiO arc-like structure in other interferometric observations (Bjerkeli et al. 2016;Yen et al. 2017).Okoda et al. (2021) have pointed out that the faint CO emission would be a remnant of part of the arc-like structure.The location and velocity of the CO emission detected in the previous works are mostly coincident with those of the SEB lobe, although the SEB lobe also traces emission on larger scales than the previous observations.Thus, the SEB lobe could be a dissipating remnant of the secondary outflow.
The NR2 lobe appears aligned with the direction of the SEB lobe and the secondary outflow, suggesting that the NR2 lobe is the counterpart of the secondary outflow.The similar velocities of the NR2 and SEB lobes of ∼ 1-2 km s −1 with respect to the systemic velocity also support this possibility.Numerical simulations suggest that the outflow can be asymmetric in the presence of turbulence because of the interaction between the outflow and an inhomogeneous envelope (Offner et al. 2011).Offner et al. (2011) show that the lengths of two lobes of an outflow with velocities of a few km s −1 can differ by more than a few thousands au.One might wonder whether the NR1 and NR2 lobes originate from a cavity wall of a single outflow.However, their distinct expansion velocities suggest that they trace different ejection events, as was mentioned in Section 4.1.Hence, the NR1 and NR2 lobes are most likely two different outflows but not different sides of the same outflow cavity.The dynamical timescale of the secondary outflow is estimated to be ∼ 5 × 10 3 yr without the correction of the inclination angle (Okoda et al. 2021), which is shorter than that of the NR2 lobe.This would be because the blueshifted lobe of the secondary outflow seems to have experienced shock (Okoda et al. 2021) and has not propagated to radii larger than ∼ 3000 au.Since the dynamical timescales for the SB, NR1 and NR2 are only lower limits, it is challenging to further discuss which outflow, the second or the third, is ejected earlier or later with the current datasets.
The configurations of the primary, secondary and third outflows on the plane of the sky are summarized in Figure 7 and Table 3. Considering the position angles of the redshifted and blueshifted lobes of the outflows, the misalignment angles between these three outflows on the plane of the sky are ∼ 20-90 • .

Origin of the Multiple Outflows
As discussed in Okoda et al. (2021), IRAS 15398 has been known as a single protostar so far.Ohashi et al. (2023) and Thieme et al. (2023) have recently presented 1.3 mm continuum images at an angular resolution of 0. ′′ 04 (∼6 au), indicating that the source is not a binary with a separation wider than ∼ 6 au.The presence of a Keplerian disk with a radius of 30 au around IRAS 15398 have been suggested (Okoda et al. 2018;Thieme et al. 2023).Observations of close binary systems show that most of them having a circumbinary disk with such a small binary separation exhibit a single outflow or parallel outflows (Lee et al. 2015;Tobin et al. 2016;Lim et al. 2016;Alves et al. 2017), as pointed out by Okoda et al. (2021).Thus, these multiple outflows would not be likely driven by unresolved very close binary protostars.
Another possible explanation of the multiple outflows around IRAS 15398 is time variation of the outflow orientation.The orientation of the outflow is expected to vary with time when the angular momentum axis of the infalling material is randomly aligned.Indeed, the dynamical timescales of the outflows are different by an order of magnitude, as discussed above.Theoretical works predict that such time variation of the outflow orientation occurs when dense cores are turbulent or the associated magnetic field is misaligned with the initial rotational axis of the dense core (Matsumoto & Tomisaka 2004;Matsumoto et al. 2017;Machida et al. 2020).In these simulations, the orientation of the disk normal changes with time because angular momentum vectors of infalling gas elements are not aligned with a single direction.Consequently, the outflow orientation also changes with time, since the outflow is generally launched in the direction of the disk normal.In the following subsections, whether the turbulence and misalignment scenarios fit IRAS 15398 will be discussed one by one.

Turbulence Scenario
Numerical simulations show that the disk orientation, which would determine the outflow direction, varies over time significantly, sometimes by more than 180 • , in turbulent dense cores (Bate et al. 2010;Seifried et al. 2013;Offner & Chaban 2017;Matsumoto et al. 2017;Bate 2018).In turbulent dense cores, material infalling onto the disk has a random angular momentum vector and infalling material having larger angular momentum alters the net angular momentum vector of the disk.Some of these simulations yield that the disk orientation can change with time in dense cores when the ratio of the turbulent energy to the gravitational energy (β turb ) is ∼ 0.06-0.09(Seifried et al. 2013;Matsumoto et al. 2017).
For a quantitative comparison to those simulations, we calculate β turb in the dense core of IRAS 15398 assuming the uniform density distribution: where Here, E turb is the turbulent energy, E grav is the gravitational energy, M core is the core mass, R core is the core radius, G is the gravitational constant, and δv turb is the turbulent velocity.The core mass is estimated from the JCMT 850 µm flux, which is presented by Yen & the eDisk team (2024), with the following equation: where F ν is the flux density, d is the distance, κ ν is the dust opacity, B ν is the Planck function, and T dust is the dust temperature.F ν is measured to be 2.9 Jy with a typical flux uncertainty of 10%.We adopt the dust opacity κ ν = 0.02 g −1 cm 2 (Ossenkopf & Henning 1994), and the dust temperature of 15 ± 5 K.The core radius is estimated to be 6700 ± 800 au from the 3σ contour of the 850 µm continuum emission (Yen & the eDisk team 2024).The uncertainty of the core radius is the standard deviation of the radii of the 3σ contour, which vary with the azimuthal angle.We used the 2 ′ square map of the C 18 O J = 2-1 emission presented by Sai et al. (2023) to measure the turbulent velocity in two ways.One way is to measure the full width at half maximum (FWHM) of a spectrum and subtract the thermal line width.The FWHM of the spectrum is measured to be 0.50 ± 0.02 km s −1 within a circle with a radius of 8 ′′ at an off position, 40 ′′ north of the protostellar position, where systemic motions such as infall or rotation are not observed.By assuming the gas temperature same as the dust temperature above, the non-thermal velocity (σ nth ) is estimated to be 0.20 ± 0.02 km s −1 .This would give an upper limit of the turbulent velocity of the dense core, as the C 18 O J = 2-1 emission may include the surrounding cloud component.These parameters yield β turb ∼ 1.4 ± 0.9, where the uncertainty is calculated through the error propagation of the uncertainty of each parameter.The calculated turbulent energy against the gravitational energy is much larger than those in numerical simulations reporting the time variation of the orientation of the disk normal.
The above calculation of the turbulent velocity may overestimate the turbulence energy in the dense core.The other way to estimate the turbulent velocity is to calculate the velocity deviation in the moment 1 map across the dense core diameter, since it is suggested that the velocity gradient on scales of ≳ 1000 au originates from turbulent motion (Sai et al. 2023).Sai et al. (2023) provide a scaling relation between the velocity deviation (δv) and the spatial scale (τ ) of δv = 10 −1.69±0.08 τ 1000 au 0.6±0.2(km s −1 ).( 6) This relation could provide a better estimation on the turbulent velocity on a specific spatial scale, because the line-of-sight components that could be larger than the map area are averaged out in the moment 1 map.By adopting this scaling relation and τ = 2R core , the turbulent velocity across the dense core is estimated to be 0.10 ± 0.01 km s −1 .This estimate results in β turb ∼ 0.33 ± 0.20, which is also larger than those of the simulations.Hence, considering the strength of the turbulence in the dense core, turbulence could cause the dynamical time variation of the outflow direction in the dense core of IRAS 15398.

Misalignment Scenario
A misalignment between the rotational axis of the initial dense core and the associated magnetic field could result in a significant change in the outflow direction over time.Machida et al. (2020) show that the orientations of the disk normal and the low-velocity outflow (< 10 km s −1 ) change by ∼ 90 • when the initial misalignment angle between the magnetic field and the rotational axis of the dense core is > 80 • .
Previous observations with SOFIA reveal orientations of magnetic fields on scales of ∼ 9000 au around IRAS 15398, providing a mean position angle of the magnetic field is ∼ 225 • (Redaelli et al. 2019).On the other hand, the JCMT/POL-2 map of 850 µm continuum emission suggests that the overall direction of the magnetic field on scales of ∼ 6000 au is north-south (P.A. of 180 • ; Yen & the eDisk team 2024).The angular momentum vector of the initial dense core could be estimated from the velocity gradient on scales of ∼ 6000-9000 au around IRAS 15398.Sai et al. (2023) have revealed that the P.A. of the velocity gradient varies from ∼ 107 to 166 • at radii of 6000-9000 au, using the C 18 O J = 2-1 emission.These gradient directions imply that P.A. of the angular momentum vector of the dense core is ∼ 200-260 • , which is orthogonal to the gradient directions.In Figure 7, we summarize these possible orientations of the magnetic field and the angular momentum vector of the dense core on the plane of the sky with the configurations of the multiple outflows.
The misalignment angle between the magnetic field and the angular momentum vector can be 20-35 • on the plane of the sky, when comparing the magnetic field orientation and velocity gradient measured on the same spatial scale.The probability of a significant misalignment between the magnetic field and the rotational axis of the dense core in 3D space, given a projected misalignment angle, can be estimated a statistical approach.Assuming that the magnetic field and the rotational axis of the dense core are randomly oriented in 3D space, the possibility that the misalignment angle is > 80 • in 3D space and appears to be 20-35 • on the plane of the sky after projection is less than 10% (Galametz et al. 2018;Gupta et al. 2022).Thus, the magnetic field is less likely largely misaligned with the rotational axis of the dense core of IRAS 15398, although the possibility of a large misalignment of > 80 • cannot be completely ruled out.

How Common are Multiple Outflows around
Single Protostars?

Is IRAS 15398 a Peculiar Source?
Multiple, highly-misaligned outflows possibly launched by single protostars, such as those found in IRAS 15398, have not been observed except for very few recent reports (Okoda et al. 2021;Sato et al. 2023), although it should be noted that multiple outflows have not been searched for in a systemic way.As discussed in the previous section, the multiple outflows in IRAS 15398 could be attributed to turbulence in its dense core.A question worth discussing is whether the turbulence in the dense core of IRAS 15398 is notably stronger compared to that in other dense cores.
We compare the turbulent velocity and energy of the dense core of IRAS 15398 with those of the dense cores in the Perseus star-forming regions, which have been well characterized with comprehensive surveys in both (sub)mm line and continuum emission (Kirk et al. 2006(Kirk et al. , 2007)).To calculate the turbulent energy and the gravitational energy for the Perseus dense cores, we compiled the JCMT 850 µm fluxes, the core radii and the total line widths of the C 18 O J = 2-1 and N 2 H + J = 1-0 lines from the catalog provided by Kirk et al. (2007).Then, we calculated the core masses and non-thermal velocities assuming the same dust opacity and temperature as for IRAS 15398.We adopted the distance of 282 ± 28 pc for the Perseus dense cores, derived from the mean and standard deviation of Gaia distances to different regions of the Perseus star-froming region (Zucker et al. 2022).Note that Kirk et al. (2007) identified the dense cores using the Clumpfind algorithm and adopted the 3σ contour edge as the definition of the core radius.Our definition of the core radius for IRAS 15398 is almost identical with theirs, and the rms of the 850 µm map for IRAS 15398 falls within a typical range of rms noises of their 850 µm maps.First, the C 18 O non-thermal velocities (σ nth ) and β turb calculated with them are compared between IRAS 15398 and the Perseus dense cores in Figure 8.The core mass, turbulent velocity and β turb of IRAS 15398 are comparable to or less than the typical values in both the starless and protostellar dense cores in the Perseus star-forming region.
The C 18 O emission observed by single-dish telescopes may include the surrounding cloud components, and the turbulent velocity could be overestimated.Thus, we performed similar comparisons using the turbulent velocities of the Perseus dense cores measured with the N 2 H + emission, which typically traces denser regions than the C 18 O emission (Kirk et al. 2007).Because there is no available N 2 H + data taken by single-dish telescopes for IRAS 15398, we used the turbulent velocity measured with the scaling relation of Equation ( 6).Line widths of the N 2 H + emission of dense cores are typically smaller than those of the C 18 O emission by a factor of two (Johnstone et al. 2010).The turbulent velocity measured with Equation ( 6) is also about twice smaller than that measured with the C 18 O emission, suggesting that our method is tracing similar spatial scales to those of the N 2 H + emission.Figure 9 presents histograms of these parameters.The core mass, turbulent velocity, and β turb of IRAS 15398 are similar to or less than the typical values in the Perseus dense cores.These comparisons suggest that turbulence of the dense core of IRAS 15398 is not notably stronger compared to other dense cores.Furthermore, most of the protostellar and starless cores in the Perseus star-forming region exhibit β turb ≳ 0.1, which is larger than those of numerical simulations reporting time evolution of the disk and/or outflow orientations (Seifried et al. 2013;Matsumoto et al. 2017).This may suggest that protostellar systems in these dense cores also experience changes in disk and outflow directions over time, which contradicts the few reports of multiple outflows.Based on the results of the numerical simulations and the observed distributions of the turbulence energy of the dense cores, we would expect that the disk/outflow orientations should evolve with time and multiple outflows form in many sources.A possible explanation for few reports of multiple outflows is that the significant change in the outflow direction only occurs within a short period at an early phase of the star formation process.While the disk would be small and its orientation would be easily altered at the beginning of the star formation process, it could be difficult to change the orientation of a larger, more massive disk at later evolutionary stages.If this is the case and the outflow launched along different directions in an early phase also dissipates within a similar timescale, then multiple outflows would be only observed in young sources.

Age of the System and Outflow Dissipation
IRAS 15398 has been considered to be at an early stage of the protostellar phase because of its small protostellar mass (Okoda et al. 2018;Thieme et al. 2023).Thieme et al. (2023) estimated the protostellar age (t * ) from the protostellar mass and mass accretion rate ( Ṁacc ) by t * = M * / Ṁacc , assuming a constant accretion rate over time.They obtained the mass accretion rate of (1.3-6.1)× 10 −6 M ⊙ yr −1 adopting the accretion luminosity of L acc ∼ L bol = 1.4 L ⊙ , M * ∼ 0.02-0.1 M ⊙ , and the protostellar radius of R * = 3 R ⊙ through the equation of Ṁacc = L acc R * /(GM * ), where G is the gravitational constant.This mass accretion rate and the protostellar mass yield t * ∼ (0.4-7.5) × 10 4 yr.The actual mass accretion rate would be variable over time.Theoretical works predict that the mass accretion rate for low-mass protostars is typically a few ×10 −6 -10 −5 M ⊙ yr −1 and beyond 10 −5 M ⊙ yr −1 when strong accretion outbursts occur via the disk gravitational instability (Dunham & Vorobyov 2012).The mass accre-tion rate estimated for IRAS 15398 with the current luminosity is closer to the accretion rate at a steady state in the numerical simulations.The actual protostellar age could be further smaller than the above estimate, if the mass accretion rate of IRAS 15398 is higher during some periods due to accretion outbursts.Thieme et al. (2023) also estimated the protostellar age to be t * ∼ (0.6-1.9) × 10 4 yr from comparison between the specific angular momentum measured for the envelope of IRAS 15398 and a core collapse model provided by Takahashi et al. (2016), which assumes the angular momentum conservation during the core collapse.These protostellar ages of ∼ 10 4 yr are smaller than the typical lifetime of the Class 0 protostars of ∼ 10 5 yr (Kenyon et al. 1990;Evans et al. 2009;Dunham et al. 2015), suggesting that IRAS 15398 is at an early stage of the protostellar phase.
The dynamical timescales of the secondary and third outflows are τ dyn ≳ 0.5-1 × 10 4 yr, which is comparable to the protostellar age of IRAS 15398.This suggests that the second and third outflows have been ejected at a very early phase just after formation of the protostar.The misaligned outflows ejected at such an early phase would dissipate within a short timescale, because there is no more momentum injection in the same direction.The dissipation timescale of the observed outflows can be estimated by equating their momentum and the drag force that the outflows experience (Codella et al. 1999).The drag force can be written as , where C D is a drag coefficient, ρ core is the density of the surrounding dense core, v out is the outflow velocity and S is the cross section area.The momentum of the outflow is P out = M out v out , where P out and M out are the momentum and mass of the outflow, respectively.By relating these two equations, the dissipation timescale of the outflow is obtained as: where ρ out is the mean density of the outflow, l out is the outflow length in 3D space, l ′ out and v ′ out are the outflow length and velocity on the plane of the sky, and i is the inclination angle.Here, we assume C D = 1 (Codella et al. 1999) and M out ∼ ρ out l out S. The typical core density is n core ∼ 1.0 × 10 4 cm −3 (e.g., Ward-Thompson et al. 2007), corresponding to ρ core ∼ 4.7 × 10 −20 g cm −3 assuming the mean molecular weight of 2.8 (Kauffmann et al. 2008).An averaged mass of the outflows are estimated to be ∼ (0.32-1.1) × 10 −3 M ⊙ from the mean integrated flux of ∼ 96 Jy km s −1 (Figure 2), assuming T ex = 50-200 K (van Kempen et al. 2009;Yıldız et al. 2012Yıldız et al. , 2013;;Dunham et al. 2014), a CO molecular abundance of 10 −4 , the mean molecular weight of 2.8 and the local thermal equilibrium (LTE).Note that the mass of the outflows is likely underestimated, because we do not take into account low-velocity components obscured by the optically thick, extended foreground gas.Offner et al. (2011) found that the outflow mass is underestimated by a factor of ∼ 5-10 by neglecting low-velocity emission with |v| ≲ 2 km s −1 based on synthetic observations of numerical simulations.Hence, we adopt outflow mass of ∼ 1 × 10 −2 M ⊙ for an order of magnitude estimate.By adopting outflow lengths of ∼ 8500 au (or ∼ 55 ′′ ) and widths of ∼ 1600 au (or ∼ 10 ′′ ) based on Figure 2, the outflow density is calculated to be ∼ 9 × 10 −20 g cm −3 .The ratio of the outflow density to the core density (ρ out /ρ core ) of ∼ 2 calculated from the above densities is roughly consistent with the density ratio found in numerical simulations (Machida & Matsumoto 2012).The density ratio, the outflow length, the outflow velocity of ∼ 3 km s −1 and the inclination angle of i ∼ 70 • yield the dissipation timescale of ∼ 2 × 10 4 yr.The calculated dissipation timescale is comparable with the dynamical timescale and the age of the system.Thus, these misaligned, multiple outflows that may form at an early stage of the star formation process would dissipate on a similar timescale to the current age of IRAS 15398.
As discussed above, dense cores including IRAS 15398 commonly exhibit turbulent energy large enough to form multiple outflows, and a possible reason why the multiple outflows are rarely observed is that they form only at an early stage of the star formation process and dissipate in a short timescale.The ratio of the age of IRAS 15398 to the typical lifetime of protostars is about ∼ 1/10, suggesting that one of ten protostellar systems may have similar multiple outflows.Systemic investigations on presence/absence of multiple outflows around single protostars are required to further test this hypothesis.

CONCLUSIONS
We reported our finding of the third outflow, as well as the possible redshifted lobe of the secondary outflow, around a single protostar IRAS15398−3359 based on mosaic observations in the CO J = 2-1 emission using the ACA 7-m array.Main results and conclusions are summarized below: 1. Our ACA 7-m array observations have newly revealed low-velocity lobes or filamentary structures at distances from the protostar larger than ∼ 30 ′′ in addition to the known high-velocity primary and blueshifted secondary outflows.
2. Two lobes, named as NR1 and SB, appear to be symmetric with respect to the protostar and exhibit Hubble-law like velocity profiles, suggesting that the NR1 and SB lobes likely trace the third outflow.The SEB lobe corresponds to the known secondary outflow.The NR2 lobe, which aligned in a line passing through the SEB lobe and the protostar, also exhibits a Hubble-law like velocity profile, suggesting that this is the counterpart of the blueshifted secondary outflow.
3. The dynamical timescales of the third outflow is estimated to be ≳ 5 × 10 3 yr.It suggests that the third outflow would be the oldest one and followed by the primary outflows with the dynamical timescale of ∼ 500 yr.These multiple outflows are misaligned with each other by ∼ 20-90 • on the plane of the sky.IRAS 15398 has been found to be a single protostar with the observation at ∼ 6 au resolution.These observations suggest that the highly-misaligned, multiple outflows are a result by the evolution of the disk and outflow orientations.
4. The ratio of the turbulent energy to the gravitational energy in the dense core of IRAS 15398 is calculated to be β ∼ 0.33-1.4based on previous observations.The turbulent energy in IRAS 15398 is sufficiently large compared to those in numerical simulations reporting the evolution of the disk/outflow orientations, while the magnetic field is unlikely to be largely misaligned (≳ 80 • ) with the rotational axis of the dense core.Hence, the multiple outflows in IRAS 15398 could be attributed to the turbulence in its dense core. 5.The turbulent energy of the dense core of IRAS 15398 is within a typical range of those of other dense cores in the Perseus star-forming region.These comparisons hint that the significant change in the disk and outflow directions over time can be a common event at the early stage of the star formation process.The few identifications of multiple outflows could be due to the short timescale of the formation and dissipation of the misaligned outflows.
This paper used the following ALMA data: ADS/JAO.The velocity channel maps of the CO J = 2-1 emission are presented in Figure A1.

B. WIND-DRIVEN SHELL MODEL FITTING
We present the fitting method of the wind-driven shell model (Shu et al. 1991;Lee et al. 2000) in this section.The method consists of two steps: fitting of the morphology and fitting of the velocity structure.In the  first step, we fit the outflow shell morphology using the moment 0 maps presented in Figure 2. We use the 5σ contours of the observed emission to trace outflow shell morphologies.First, the moment 0 maps are rotated by P.A. of the outflow lobes, which is measured as lines passing through the protostellar position and the emission peak of the outflow lobes, so that the outflow axis corresponds to the y axis of the plane of the sky (Figure B1).Then, we extract intensity profiles along the x axis at y offsets sampled by a half of the beam size, and measure positions of the intensity at the 5σ level.Uncertainties of the measured positions are calculated as positional offsets at 5 ± 1σ contours from the 5σ contour.The NR1 and NR2 lobes spatially overlap.Therefore, we separate them at local intensity minima of the extracted intensity profiles.At a given y offset, two x values corresponding to the outflow shell edges on different sides are obtained.To better trace the outflow shell morphology, the inner data point is removed when both two points are on the same side (i.e., when both are x < 0 or x > 0).Similarly, data points at artificial edges, i.e., the local intensity minima where the lobes are separated, are removed.The projected outflow shell morphology is written as y = C sin ix 2 − cos 2 i 4C sin i , where x and y are plane-of-sky coordinates.The inclination angle is defined as an angle between the line-ofsight axis and the outflow axis.We search for the parameter set of (C sin i, cos i) to best explain the observations by fitting the data points with the Markov-Chain Monte Carlo (MCMC) code emcee (Foreman-Mackey et al. 2013).Uncertainties of the fitting parameters are calculated as the standard deviation of the posterior distributions.The obtained data points and fitting results are shown in Figure B1.The posterior distributions are also presented in Figure B2.The cavity shape, C sin i, is well constrained, while the cosine of the inclination angle, cos i, does not converge.Thus, we leave the inclination angle as a free parameter in the second fitting step.
In the second step, we fit the velocity structure using PV diagrams and C sin i constrained in the first step.The 5 ± 1σ contours are used to measure velocities and their uncertainties at offsets sampled by a half of the beam size.Around the systemic velocity, the emission is highly affected by the cloud contamination and its resolve-out.Thus, data points at velocities less than 4.4 km s −1 are removed from the fitting for the SB lobe to better trace the velocity structure of the outflow.Similarly, data points at velocities larger than 6.34 km s −1 and 6.25 km s −1 are removed from fits for the NR1 and NR2 lobes, respectively.We perform the fitting of the velocity structure in the same manner as in the first step.The posterior distributions of the parameters are shown in Figure B3, and fitting results are presented in Figures 3 and 4, and in Table 1.Note that these fitting results do not change significantly even if we include the data points that are likely affected by the cloud contamination.The fitting results using all the data points are presented in Figure B4 and Table B1.The uncertainties of v 0 and i shown in Table 1 and B1 do not take into account systemic errors of the fitting method itself.Based on the differences of v 0 and i values between Table 1 and B1, systemic errors of v 0 and i associated with the selection of the data point around the cloud velocity are about 10% and 1%, respectively.

Figure 2 .
Figure 2. (a) High-and (b) low-velocity components of the CO J =2-1 emission.Red and blue colors show redshifted and blueshifted components, respectively.The LSR velocity range of each velocity component is shown in the top left corners in units of km s −1 .Contour levels are 3, 6, 12, 24, ... ×σ, where 1σ is the rms noise level.The rms of each map is 0.16 Jy beam −1 for the high-velocity redshifted and high-velocity blueshifted components, 0.20 Jy beam −1 for the low-velocity redshifted component, and 0.46 Jy beam −1 for the low-velocity blueshifted component.The central crosses denote the protostellar position.The ellipses at the bottom left corners indicate the synthesized beam size of 7. ′′ 7 × 4. ′′ 1 (87 • ).The dashed lines in panel b show directions of the position-velocity (PV) cuts.

SBFigure 3 .NR2Figure 4 .Figure 5 .
Figure 3. PV diagrams of the CO J =2-1 emission cut toward the (a) SB and (b) NR1 components.Contour levels are from 3, 6, 12, 24, ... ×σ, where 1σ is 0.13 Jy beam −1 .Black dots are data points tracing 5σ contours that are used to fit the wind-driven shell model.Most of error bars of the data points are smaller than the marker size.Black curves show velocity structures of the best-fit wind-driven shell models.The vertical and horizontal bars at the bottom left corners indicate the synthesized beam size of 7. ′′ 7 × 4. ′′ 1 (87 • ) and the velocity resolution of 0.16 km s −1 .The vertical and horizontal dashed lines denote the systemic velocity of 5.4 km s −1 and the protostellar position, respectively.

Figure 5
Figure5shows a PV diagram cut toward the SEB component (P.A.= 155 • ).While the emission within an offsets of ∼ 10 ′′ is the high-velocity component of the primary outflow, the SEB component appears at

Figure 6 .
Figure 6.(a) Same as Figure 2(b) but zoomed into the WB and WR filamentary structures.Blue/red dashed lines and dots denote the directions and centers of the PV cuts, respectively.(b) PV diagrams of the CO 2-1 emission cut along the WB and (c) WR filamentary structures centered their emission peaks.The vertical and horizontal bars at the bottom left corners and contour levels are the same as those in Figure 3.The vertical and horizontal dashed lines denote the systemic velocity of 5.4 km s −1 and the emission peak of the filamentary structures.

Figure 7 .
Figure 7. Summary of the configurations of the multiple outflows and the possible orientations of the magnetic field and the angular momentum vector of the dense core on the plane of the sky.Lobes drawn with solid and dashed curves indicate the redshifted and blueshifted components of the outflows, respectively.

Figure 8 .
Figure8.Histograms of the mass, non-thermal velocity, and the ratio of the turbulence energy to the gravitational energy of the dense cores in the Perseus star forming region from left to right, which are referred fromKirk et al. (2007).Vertical lines show the parameters for the dense core of IRAS 15398 and grey shaded area indicate their uncertainties.Non-thermal velocities are measured with the FWHM of the C 18 O J = 2-1 spectrum for both Perseus dense cores and the dense core of IRAS 15398.

Figure 9 .
Figure9.Same as Figure8but with measurements using the N2H + J =1-0 emission for the dense cores in the Perseus star-forming region, referred fromKirk et al. (2007), and with the non-thermal velocity measured using the scaling relation of Equation (6) for IRAS 15398.

Figure A1 .
Figure A1.Velocity channel maps of (a) blueshifted and (b) redshifted components of the CO J = 2-1 emission.Contour levels are 3, 6, 12, 24, ... ×σ, where 1σ is 0.13 Jy beam −1 .Maps are shown in steps of two channels.The label at the top left corner in each panel denotes the LSR velocity of the channel in km s −1 .The central crosses denote the protostellar position.The ellipse at the bottom left corner indicates the synthesized beam size of 7. ′′ 7 × 4. ′′ 1 (87 • ).

Figure B1 .Figure B2 .
Figure B1.Results of the fitting of the outflow shell morphology.Black points are measured data points corresponding to 5σ contours.Black curves indicate the best-fit outflow shell morphologies.
Figure B3.Same as Figure B2 but for the fitting of the velocity structure.

Figure B4 .
Figure B4.Same as Figures 3 and 4 but with fitting results for all data points including those likely affected by the cloud contamination.

Table 1 .
Results of Wind-driven Shell Model Fitting

Table 2 .
Lengths, Velocities and Dynamical Timescales of the Outflow Lobes

Table B1 .
Results of Wind-driven Shell Model Fits for All Data Points