Multiple Outflows around a Single Protostar IRAS
15398-3359

Jinshi Sai, Hsi-Wei Yen, Masahiro Machida, Nagayoshi Ohashi, Yusuke Aso,

Anaëlle Maury, Sébastien Maret

► To cite this version:

Jinshi Sai, Hsi-Wei Yen, Masahiro Machida, Nagayoshi Ohashi, Yusuke Aso, et al. Multiple Outflows around a Single Protostar IRAS 15398-3359. The Astrophysical Journal, 2024, 966 (2), pp.192-215. arXiv:2402.04357 [astro-ph.HE] 10.3847/1538-4357/ad34b7. hal-04803311

HAL Id: hal-04803311

https://hal.science/hal-04803311v1

Submitted on 25 Nov 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# Multiple Outflows around a Single Protostar IRAS 15398–3359

Jinshi Sai<fnref n="1" />, Hsi-Wei Yen<fnref n="1" />, Masahiro N. Machida<fnref n="2" />, Nagayoshi Ohashi<fnref n="1" />, Yusuke Aso<fnref n="3" />, Anaëlle J. Maury<fnref n="4" />,<fnref n="5" />, and Sébastien Maret<fnref n="6" />

Received 2024 February 7; revised 2024 March 14; accepted 2024 March 15; published 2024 May 7

## Abstract

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.

*Unified Astronomy Thesaurus concepts: Protostars (1302); Young stellar objects (1834); Star formation (1569); Stellar winds (1636); Millimeter astronomy (1061); Submillimeter astronomy (1647); Radio astronomy (1338); Stellar accretion disks (1579)*

## 1. 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 the 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 (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

the magnetic braking and allow disks to form (e.g., Joos et al. 2012; Seifried et al. 2013; Matsumoto et al. 2017). Nonideal 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 presences 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; 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 as well as disk orientations are all frequently found to be misaligned (Galametz et al. 2018, 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 still remains 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

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Academia Sinica Institute of Astronomy and Astrophysics, 11F of Astro-Math Bldg, 1, Sec. 4, Roosevelt Rd, Taipei 10617, Taiwan</footnote>

<footnote n="2">² Department of Earth and Planetary Sciences, Faculty of Sciences, Kyushu University, Fukuoka, Fukuoka 819-0395, Japan</footnote>

<footnote n="3">³ Korea Astronomy and Space Science Institute, 776 Daedeok-daero, Yuseong-gu, Daejeon 34055, Republic of Korea</footnote>

<footnote n="4">⁴ AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Citè, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="5">⁵ Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138, USA</footnote>

<footnote n="6">⁶ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>
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 a 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_\odot$, 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 by Okoda et al. (2021) through molecular line observations with high sensitivity using the Atacama Large Millimeter/submillimeter Array (ALMA). 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 $\lesssim 0.01 M_\odot$, based on observations at angular resolutions of $\lesssim 0.02''$ (Yen et al. 2017; Okoda et al. 2018). Recent observations at an angular resolution of $\sim 0.02''$ have revealed a possible differential rotation signature of a Keplerian disk and constrained the central protostellar mass range of $\sim 0.02–0.1 M_\odot$ 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 we discuss the possible origins of the misaligned multiple outflows. The finding of the largely misaligned multiple outflows in IRAS 15398 hints at 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. The observations and data reduction are summarized in Section 2. The observational results are presented in Section 3. Analysis of the velocity structures is provided in Section 4, and possible implications are discussed in Section 5. The main results and conclusions are summarized in Section 6.

## 2. Atacama Compact Array Observations

We conducted mosaic observations of IRAS 15398 using the 7 m array of the 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.538000 GHz) results of the observations. The mosaic observations with 28 pointings covered a $\sim 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. 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 $\sim 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 pipeline version 42866M (Pipeline-CASA56-P1-B) of the Common Astronomy Software Applications (CASA) package (McMullin et al. 2007), version 5.6.1.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Moment 0 (contours) and 1 (color) maps of the CO $J=2–1$ emission of IRAS 15398–3359. The contour levels are $3\sigma$, $6\sigma$, $12\sigma$, $24\sigma$, ..., where $1\sigma$ is $0.58$ Jy beam$^{-1}$. The central cross denotes the protostellar position. The ellipse at the bottom left corner indicates the synthesized beam size of $7''$ $7 \times 4''$ $1$ (87°).</figcaption></figure>

We produced images from the visibility data using CASA 5.6.1 with the `tclean` task. We adopted the Briggs-style weighting, with a robust parameter of 0.5 and a velocity resolution of $0.08$ km s$^{-1}$. A multiscale clean was used, with scale parameters of zero and 10, which correspond to the point source and about one beam size, respectively. 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 rms noise level are $7''$ $7 \times 4''$ $1$ (87°) and $0.13$ Jy beam$^{-1}$, respectively. We adopted the source position of $(\alpha_{J2000}, \delta_{J2000}) = (15:43:02.24, -34:09:06.8)$ for the figures and analyses in this paper (Yen et al. 2017).

## 3. Results

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 with a length of $\sim 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 the blueshifted and redshifted components have a relatively high velocity of $\pm 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 an LSR velocity of $\sim 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 (position angle or PA $\sim 140°$; Okoda et al. 2021), although its size is larger than the size of the secondary outflow. At distances from the protostar larger than $30''$, blueshifted and redshifted components are
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) High-velocity and (b) low-velocity components of the CO J = 2–1 emission. The red and blue colors show the redshifted and blueshifted components, respectively. The LSR velocity range of each velocity component is shown in the top left corners in units of kilometers per second. The contour levels are 3σ, 6σ, 12σ, 24σ, ..., where 1σ is the rms noise level. The rms of each map is 0.16 Jy beam⁻¹ for the high-velocity redshifted and high-velocity blueshifted components, 0.20 Jy beam⁻¹ for the low-velocity redshifted component, and 0.46 Jy beam⁻¹ 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 (b) show the directions of the PV cuts.</figcaption></figure>

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⁻¹ with respect to the systemic velocity. These features suggest that they are a 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 the velocity channel maps shown in Figure A1). The high-velocity components shown in Figure 2(a) correspond to the elongated emission 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⁻¹ and from 8.3 to 15.0 km s⁻¹, 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⁻¹ and ~1.7–4.9 km s⁻¹, 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.8 km s⁻¹, which are closer to the systemic velocity of 5.4 km s⁻¹ 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⁻¹, than the NR1 and NR2 components.

# 4. Analysis

## 4.1. Velocity Structure

The velocity structures of the low-velocity components are investigated in more detail in this section. Figure 3 presents

position-velocity (PV) diagrams cut toward the NR1 (PA = 355°) and SB (PA = 182°) components. The PA 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⁻¹. This would be because the optically thick extended foreground 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⁻¹ and 6–10 km s⁻¹, and it corresponds to the high-velocity components of the primary outflow. The NR1 and SB components are found at offsets of 20″–75″. They exhibit 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⁻¹ 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 (PA = 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⁻¹ as the offset increases. The Hubble-law-like velocity profile suggests that the NR2 component is also part of an outflow. Although its symmetrical counterpart is not found, the NR2 and SEB components are almost aligned in the same line passing through the central protostar, as shown in Figure 2(b). The relationship of the NR2 component to the other outflows is discussed in more detail in Section 5.1.

Figure 5 shows a PV diagram cut toward the SEB component (PA = 155°). While the emission within an offset of ~10″ is the high-velocity component of the primary outflow, the SEB component appears at offsets of ~10″–20″ and
<figure><img src="image_6.png" /><figcaption>**Figure 3.** PV diagrams of the CO J = 2–1 emission cut toward the (a) SB and (b) NR1 components. The contour levels are from 3σ, 6σ, 12σ, 24σ, ..., where 1σ is 0.13 Jy beam⁻¹. The black dots are data points tracing 5σ contours that are used to fit the wind-driven shell model. Most of the error bars of the data points are smaller than the marker size. The 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⁻¹, respectively. The vertical and horizontal dashed lines denote the systemic velocity of 5.4 km s⁻¹ and the protostellar position, respectively.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** PV diagram of the CO J = 2–1 emission cut toward the NR2 component. The contour levels and symbols are the same as those in Figure 3.</figcaption></figure>

exhibits an almost constant velocity of ~4.3 km s⁻¹. 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⁻¹ are likely part of the 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. 2018, 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 that of the shock remnant caused by the secondary outflow (Okoda et al. 2021). The SEB component and secondary outflow are compared in more detail in Section 5.1.

The WB and WR filamentary structures do not appear to be pointing toward 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 PAs are measured though Gaussian fitting and are presented in Figure 6(a). Figures 6(b) and (c) show that these structures have almost constant LSR velocities around 4.3 and 6.8 km s⁻¹, respectively, along the elongation. No emission is detected from ~4.5 to 6.5 km s⁻¹ 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).

## 4.2. 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 Hubble-law-like velocity profiles. The NR2 component shows a Hubble-law-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:

$$z = CR^2, \quad v_R = v_0 \frac{R}{R_0}, \quad v_z = v_0 \frac{z}{z_0}, \qquad (1)$$

where R is the cylindrical radius in a direction perpendicular to the outflow axis, z is the distance from the protostar along the
<figure><img src="image_8.png" /><figcaption>**Figure 5.** PV diagram of the CO *J* = 2–1 emission cut along the SEB component. The contour levels and symbols are the same as those in Figure 3.</figcaption></figure>

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^\circ$ 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\sigma$ contours of the emission are used to define the 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 2, 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 $\ge 0.49-1.1 \times 10^4$ yr with correction of the inclination angle. A previous work estimated the dynamical timescale of the primary outflow to be $\sim 5 \times 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.

<table><caption>**Table 1**

Results of Wind-driven Shell Model Fitting</caption><thead><tr><th>Lobe</th><th>C sin i</th><th>v<sub>0</sub><br>(km s<sup>-1</sup>)</th><th>i<br>(deg)</th></tr></thead><tbody><tr><td>SB</td><td>0.81 ± 0.02</td><td>0.1129 ± 0.0004</td><td>71.83 ± 0.05</td></tr><tr><td>NR1</td><td>1.24 ± 0.02</td><td>0.1017 ± 0.0006</td><td>74.24 ± 0.12</td></tr><tr><td>NR2</td><td>0.89 ± 0.02</td><td>0.0542 ± 0.0016</td><td>69.75 ± 0.71</td></tr></tbody></table>

<table><caption>**Table 2**

Lengths, Velocities, and Dynamical Timescales of the Outflow Lobes</caption><thead><tr><th>Lobe</th><th>r<sub>CO</sub><br>(au)</th><th>v<sub>max,LSR</sub><br>(km s<sup>-1</sup>)</th><th>v<sub>max,rel</sub><br>(km s<sup>-1</sup>)</th><th>τ'<sub>dyn</sub><br>(10<sup>4</sup> yr)</th><th>τ<sub>dyn</sub><br>(10<sup>4</sup> yr)</th></tr></thead><tbody><tr><td>SB</td><td>≥12,000</td><td>1.72</td><td>3.68</td><td>≥1.5</td><td>≥0.49</td></tr><tr><td>NR1</td><td>≥12,000</td><td>8.12</td><td>2.72</td><td>≥2.1</td><td>≥0.59</td></tr><tr><td>NR2</td><td>≥12,000</td><td>7.32</td><td>1.92</td><td>≥3.0</td><td>≥1.1</td></tr></tbody></table>

Note. $\tau_{dyn}$ and $\tau'_{dyn}$ are the dynamical timescales with and without inclination angle correction, respectively, and $\tau_{dyn} = \tau'_{dyn}/\tan i$.

# 5. Discussion

## 5.1. 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 arclike 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 arclike structure traces an old shock caused by the secondary outflow. The velocity range of the arclike structure is almost consistent with the velocity range of the SEB lobe. The lack of 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 arclike 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 arclike structure. The location and velocity of the CO emission detected in 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 to be aligned with the direction of the SEB lobe and the secondary outflow, suggesting that the NR2 lobe is a 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 kilometers per second can differ by more than a few thousand astronomical units. 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 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
<figure><img src="image_9.png" /><figcaption>**Figure 6.** (a) The same as Figure 2(b), but zoomed in to the WB and WR filamentary structures. The 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 at their emission peaks. The vertical and horizontal bars at the bottom left corners and the 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⁻¹ and the emission peak of the filamentary structures, respectively.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**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.</figcaption></figure>

to be ~5 × 10³ 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

<table><caption>**Table 3**
PAs of the Outflows</caption><thead><tr><th colspan="2">Primary</th><th colspan="2">Secondary</th><th colspan="2">Third</th></tr><tr><th>Red<sup>1</sup></th><th>Blue<sup>1</sup></th><th>Red<sup>3</sup></th><th>Blue<sup>2</sup></th><th>Red<sup>3</sup></th><th>Blue<sup>3</sup></th></tr></thead><tbody><tr><td>60</td><td>230</td><td>332</td><td>140</td><td>355</td><td>182</td></tr></tbody></table>

**Notes.** (1) Yen et al. (2017). (2) Okoda et al. (2021). (3) This work. All values are in degrees.

(Okoda et al. 2021) and has not propagated to radii larger than ~3000 au. Since the dynamical timescales for the SB, NR1, and NR2 lobes 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 data sets.

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 PAs 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°.

## 5.2. Origin of the Multiple Outflows

As discussed in Okoda et al. (2021), IRAS 15398 has so far been known as a single protostar. 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 has 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; Lim et al. 2016; Tobin et al. 2016; Alves et al. 2017), as pointed out by Okoda et al. (2021). Thus, these multiple outflows would not likely be driven by unresolved very close binary protostars.

Another possible explanation of the multiple outflows around IRAS 15398 is the 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.

### 5.2.1. 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; Matsumoto et al. 2017; Offner & Chaban 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 the disk orientation changing with time in dense cores when the ratio of the turbulent energy to the gravitational energy ($\beta_{turb}$) is ~0.06–0.09 (Seifried et al. 2013; Matsumoto et al. 2017).

For a quantitative comparison to those simulations, we calculate $\beta_{turb}$ in the dense core of IRAS 15398, assuming the uniform density distribution:

$$ \beta_{\text{turb}} = \frac{E_{\text{turb}}}{E_{\text{grav}}}, \qquad (2) $$

where

$$ E_{\text{turb}} = \frac{3}{2} M_{\text{core}} \delta v_{\text{turb}}^2, \qquad (3) $$

$$ E_{\text{grav}} = \frac{3}{5} \frac{GM_{\text{core}}^2}{R_{\text{core}}}. \qquad (4) $$

Here, $E_{\text{turb}}$ is the turbulent energy, $E_{\text{grav}}$ is the gravitational energy, $M_{\text{core}}$ is the core mass, $R_{\text{core}}$ is the core radius, $G$ is the gravitational constant, and $\delta v_{\text{turb}}$ is the turbulent velocity. The core mass is estimated from the James Clerk Maxwell Telescope (JCMT) 850 µm flux, which is presented by Yen (2024), with the following equation:

$$ M_{850} = \frac{F_{\nu} d^2}{\kappa_{\nu} B_{\nu} (T_{\text{dust}})}, \qquad (5) $$

where $F_{\nu}$ is the flux density, $d$ is the distance, $\kappa_{\nu}$ is the dust opacity, $B_{\nu}$ is the Planck function, and $T_{\text{dust}}$ is the dust temperature. $F_{\nu}$ is measured to be 2.9 Jy with a typical flux uncertainty of 10%. We adopt the dust opacity $\kappa_{\nu} = 0.02 \text{ g}^{-1}\text{cm}^2$ (Ossenkopf & Henning 1994) and a dust temperature of $15 \pm 5$ K. The core radius is estimated to be 6700 ± 800 au from the 3σ contour of the 850 µm continuum emission (Yen and 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 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 to be the same as the dust temperature above, the nonthermal velocity ($\sigma_{\text{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 $\beta_{\text{turb}} \sim 1.4 \pm 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 $\ge 1000$ au originates from turbulent motion (Sai et al. 2023). Sai et al. (2023) provide a scaling relation between the velocity deviation ($\delta v$) and the spatial scale ($\tau$) of

$$ \delta v = 10^{-1.69 \pm 0.08} \left( \frac{\tau}{1000 \text{ au}} \right)^{0.6 \pm 0.2} (\text{km s}^{-1}). \qquad (6) $$

This relation could provide a better estimation of 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 $\tau = 2R_{\text{core}}$, the turbulent velocity across the dense core is estimated to be 0.10 ± 0.01 km s$^{-1}$. This estimate results in $\beta_{\text{turb}} \sim 0.33 \pm 0.20$, which is also larger than those of the simulations. Hence, considering the strength of the turbulence in the dense core, the turbulence could cause the dynamical time variation of the outflow direction in the dense core of IRAS 15398.

### 5.2.2. 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 PA of the magnetic field as ~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 (PA of 180°; Yen and 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 PA of the velocity gradient varies from ~107° to 166° at radii of
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Histograms of the mass, nonthermal velocity, and ratio of the turbulence energy to the gravitational energy of the dense cores in the Perseus star-forming region, from left to right, in reference to Kirk et al. (2007). The vertical lines show the parameters for the dense core of IRAS 15398 and the gray shaded areas indicate their uncertainties. Nonthermal velocities are measured with the FWHM of the C$^{18}$O J = 2–1 spectrum for both the Perseus dense cores and the dense core of IRAS 15398.</figcaption></figure>

6000–9000 au, using the C$^{18}$O J = 2–1 emission. These gradient directions imply that the PA 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 through 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 to be 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.

## 5.3. How Common Are Multiple Outflows around Single Protostars?

### 5.3.1. 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, excepting a 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) millimeter line and continuum emission (Kirk et al. 2006, 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 nonthermal velocities, assuming the same dust opacity and temperature as for IRAS 15398. We adopted a distance of 282 ± 28 pc for the Perseus dense cores, derived from the mean and standard deviation of the Gaia distances to different regions of the Perseus star-forming region (Zucker et al. 2020). 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 to theirs, and the rms of the 850 µm map for IRAS 15398 falls within the typical range of the rms noises of their 850 µm maps. First, the C$^{18}$O nonthermal velocities ($\sigma_{nth}$) and $\beta_{turb}$ calculated with them are compared between IRAS 15398 and the Perseus dense cores in Figure 8. The core mass, turbulent velocity, and $\beta_{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). The 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 2 (Johnstone et al. 2010). The turbulent velocity measured with Equation (6) is also about 2 times 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 $\beta_{turb}$ of IRAS 15398 are similar to or
<figure><img src="image_12.png" /><figcaption>**Figure 9.** The same as Figure 8, but with measurements using the $N_2H^+$ J = 1–0 emission for the dense cores in the Perseus star-forming region, in reference to Kirk et al. (2007), and with the nonthermal velocity measured using the scaling relation of Equation (6) for IRAS 15398.</figcaption></figure>

less than the typical values in the Perseus dense cores. These comparisons suggest that the 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 $\beta_{\text{turb}} \gtrsim 0.1$, which is larger than those of numerical simulations reporting the time evolution of the disk and/or outflow orientations (Seifried et al. 2013; Matsumoto et al. 2017). This may suggest that the protostellar systems in these dense cores also experience changes in disk and outflow directions over time, contradicting the few reports of multiple outflows.

### 5.3.2. Age of the System and Outflow Dissipation Timescale

Based on the results of the numerical simulations and the observed distributions of the turbulence energy of the dense cores, we would expect the disk/outflow orientations to evolve with time and multiple outflows to form in many sources. A possible explanation for the 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 only be observed in young sources.

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 ($\dot{M}_{\text{acc}}$) by $t_* = M_*/\dot{M}_{\text{acc}}$, assuming a constant accretion rate over time. They obtained a mass accretion rate of (1.3–6.1) $\times 10^{-6} M_{\odot} \text{yr}^{-1}$, adopting an accretion luminosity of $L_{\text{acc}} \sim L_{\text{bol}} = 1.4 L_{\odot}$, $M_* \sim 0.02–0.1 M_{\odot}$ and a protostellar radius of $R_* = 3 R_{\odot}$ through the equation $\dot{M}_{\text{acc}} = L_{\text{acc}}R_*/(GM_*)$, where G is the gravitational constant. This mass accretion rate and the protostellar mass yield $t_* \sim (0.4–7.5) \times 10^4 \text{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 times $10^{-6}–10^{-5} M_{\odot} \text{yr}^{-1}$ and beyond $10^{-5} M_{\odot} \text{yr}^{-1}$ when strong

accretion outbursts occur via the disk gravitational instability (Dunham & Vorobyov 2012). The mass accretion 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_* \sim (0.6–1.9) \times 10^4 \text{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 $\sim 10^4 \text{yr}$ are smaller than the typical lifetime of the Class 0 protostars of $\sim 10^5 \text{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 $\tau_{\text{dyn}} \gtrsim 0.5–1 \times 10^4 \text{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 $F_d = \frac{1}{2}C_D\rho_{\text{core}}v_{\text{out}}^2 S$, where $C_D$ is a drag coefficient, $\rho_{\text{core}}$ is the density of the surrounding dense core, $v_{\text{out}}$ is the outflow velocity, and S is the cross-sectional area. The momentum of the outflow is $P_{\text{out}} = M_{\text{out}}v_{\text{out}}$, where $P_{\text{out}}$ and $M_{\text{out}}$ are the momentum and mass of the outflow, respectively. By relating these two equations, the dissipation timescale of the outflow is obtained as:

$$t_{\text{diss}} = \frac{P_{\text{out}}}{F_D} \sim 2 \frac{\rho_{\text{out}} l_{\text{out}}}{\rho_{\text{core}} v_{\text{out}}}, \quad (7)$$

$$= 2 \frac{\rho_{\text{out}} l'_{\text{out}}}{\rho_{\text{core}} v'_{\text{out}}} \frac{1}{\tan i}, \quad (8)$$

where $\rho_{\text{out}}$ is the mean density of the outflow, $l_{\text{out}}$ is the outflow length in 3D space, $l'_{\text{out}}$ and $v'_{\text{out}}$ are the outflow length and velocity on the plane of the sky, respectively, and i is the
inclination angle. Here, we assume $C_D = 1$ (Codella et al. 1999) and $M_{\text{out}} \sim \rho_{\text{out}} l_{\text{out}} S$. The typical core density is $n_{\text{core}} \sim 1.0 \times 10^4 \text{ cm}^{-3}$ (e.g., Ward-Thompson et al. 2007), corresponding to $\rho_{\text{core}} \sim 4.7 \times 10^{-20} \text{ g cm}^{-3}$, assuming a mean molecular weight of 2.8 (Kauffmann et al. 2008). The averaged mass of the outflows is estimated to be $\sim(0.32-1.1) \times 10^{-3} M_{\odot}$ from the mean integrated flux of $\sim 96$ Jy km s$^{-1}$ (Figure 2), assuming $T_{\text{ex}} = 50-200$ K (van Kempen et al. 2009; Yıldız et al. 2012, 2013; Dunham et al. 2014), a CO molecular abundance of $10^{-4}$, a mean molecular weight of 2.8, and local thermal equilibrium. 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 $\sim 5-10$ by neglecting low-velocity emission with $|v| \lesssim 2 \text{ km s}^{-1}$, based on synthetic observations of numerical simulations. Hence, we adopt an outflow mass of $\sim 1 \times 10^{-2} M_{\odot}$ for an order-of-magnitude estimate. By adopting outflow lengths of $\sim 8500$ au (or $\sim 55''$) and widths of $\sim 1600$ au (or $\sim 10''$) based on Figure 2, the outflow density is calculated to be $\sim 9 \times 10^{-20} \text{ g cm}^{-3}$. The ratio of the outflow density to the core density ($\rho_{\text{out}}/\rho_{\text{core}}$) of $\sim 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 $\sim 3 \text{ km s}^{-1}$, and the inclination angle of $i \sim 70^\circ$ yield a dissipation timescale of $\sim 2 \times 10^4 \text{ 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 over a short timescale. The ratio of the age of IRAS 15398 to the typical lifetime of protostars is about $\sim 1:10$, suggesting that one of ten protostellar systems may have similar multiple outflows. Systematic investigations on the presence/absence of multiple outflows around single protostars are required to further test this hypothesis.

## 6. Conclusions

We have reported our finding of the third outflow, as well as the possible redshifted lobe of the secondary outflow, around a single protostar IRAS 15398–3359, based on mosaic observations in the CO $J=2-1$ emission using the ACA 7 m array. The 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 $\sim 30''$, in addition to the known high-velocity primary and blueshifted secondary outflows.

2. Two lobes, named as NR1 and SB, appear to be symmetrical 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 aligns in a line passing through the SEB lobe and the protostar, also exhibits a Hubble-law-like velocity profile, suggesting that this is a counterpart of the blueshifted secondary outflow.

3. The dynamical timescale of the third outflow is estimated to be $\gtrsim 5 \times 10^3 \text{ yr}$. This suggests that the third outflow would be the oldest one, followed by the primary outflows with a dynamical timescale of $\sim 500$ yr. These multiple outflows are misaligned with each other by $\sim 20^\circ-90^\circ$ on the plane of the sky. IRAS 15398 has been found to be a single protostar with an observation at $\sim 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 $\beta_{\text{turb}} \sim 0.33-1.4$, based 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 ($\gtrsim 80^\circ$) 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 the typical range of those of other dense cores in the Perseus star-forming region. These comparisons hint that significant changes 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.

## Acknowledgments

This paper has used the following ALMA data: ADS/JAO. ALMA \#2019.1.01063.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. We thank all ALMA staff for conducting the observations. We also thank Yuki Okoda for the discussion on the interpretation of the observational data. H.-W.Y. acknowledges support from the National Science and Technology Council (NSTC) in Taiwan through the grant NSTC 1102628-M-001-003-MY3 and from the Academia Sinica Career Development Award (AS-CDA-111-M03). N.O. acknowledges support from the NSTC through the grants NSTC 109-2112-M-001-051, 110-2112-M-001-031.

Facility: ALMA

Software: CASA (McMullin et al. 2007), Numpy (Oli-phant 2006; van der Walt et al. 2011), Scipy (Virtanen et al. 2020), Astropy (Astropy Collaboration et al. 2013, 2018), Matplotlib (Hunter 2007), emcee (Foreman-Mackey et al. 2013).

## Appendix A
Supplemental Figures

The velocity channel maps of the CO $J=2-1$ emission are presented in Figure A1.
<figure><img src="image_13.png" /><figcaption>**Figure A1.** Velocity channel maps of (a) blueshifted and (b) redshifted components of the CO J = 2–1 emission. The contour levels are 3σ, 6σ, 12σ, 24σ, ..., where 1σ is 0.13 Jy beam⁻¹. Maps are shown in steps of two channels. The label at the top left corner of each panel denotes the LSR velocity of the channel in kilometers per second. 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°).</figcaption></figure>

# Appendix 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: the fitting of the morphology and the 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 the outflow shell morphologies. First, the moment 0 maps

are rotated by the PA 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 half of the beam size, and measure positions of the intensity at the 5σ level. The uncertainties of the measured positions are calculated as positional offsets at 5 ± 1σ contours from the 5σ contour. The NR1 and NR2 lobes spatially
<figure><img src="image_14.png" /><figcaption>**Figure B1.** Results of the fitting of the outflow shell morphology. The black points are measured data points corresponding to $5\sigma$ contours. The black curves indicate the best-fit outflow shell morphologies.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure B2.** Posterior distributions of the parameters for the fitting of the shell morphology.</figcaption></figure>

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 - \frac{\cos^2 i}{4C \sin i}, \quad \therefore x = \pm \sqrt{\frac{1}{C \sin i} y + \frac{\cos^2 i}{4C^2 \sin^2 i}}, \quad (B1)$$

where x and y are plane-of-sky coordinates. The inclination angle is defined as an angle between the line-of-sight 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 code `emcee` (Foreman-Mackey et al. 2013). The 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 the $C \sin i$ constrained in the first step. The $5 \pm 1\sigma$ contours are used to measure velocities and their uncertainties at offsets sampled by 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 \text{ 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 \text{ 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
<figure><img src="image_16.png" /><figcaption>Figure B3. The same as Figure B2, but for the fitting of the velocity structure.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure B4. The same as Figures 3 and 4, but with fitting results for all data points, including those likely affected by the cloud contamination.</figcaption></figure>

<table><caption>Table B1
Results of Wind-driven Shell Model Fits for All Data Points</caption><thead><tr><th>Lobe</th><th>C sin i</th><th>v<sub>0</sub><br>(km s<sup>-1</sup>)</th><th>i<br>(deg)</th></tr></thead><tbody><tr><td>SB</td><td>0.81 ± 0.02</td><td>0.1007 ± 0.0003</td><td>69.35 ± 0.05</td></tr><tr><td>NR1</td><td>1.24 ± 0.02</td><td>0.0963 ± 0.0004</td><td>73.22 ± 0.06</td></tr><tr><td>NR2</td><td>0.89 ± 0.02</td><td>0.0487 ± 0.0005</td><td>66.66 ± 0.20</td></tr></tbody></table>

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 Tables 1 and B1 do not take into account the systemic errors of the fitting method itself. Based on the differences of the $v_0$ and $i$ values between Tables 1 and B1, the 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.

## ORCID iDs

Jinshi Sai https://orcid.org/0000-0003-4361-5577

Hsi-Wei Yen https://orcid.org/0000-0003-1412-893X

Masahiro N. Machida https://orcid.org/0000-0002-0963-0872

Nagayoshi Ohashi https://orcid.org/0000-0003-0998-5064

Yusuke Aso https://orcid.org/0000-0002-8238-7709

Anaëlle J. Maury https://orcid.org/0000-0002-3801-8754

Sébastien Maret https://orcid.org/0000-0003-1104-4554

## References

* Allen, A., Li, Z.-Y., & Shu, F. H. 2003, ApJ, 599, 363

* Alves, F. O., Girart, J. M., Caselli, P., et al. 2017, A&A, 603, L3

* Andrews, S. M. 2020, ARA&A, 58, 483

* Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

* Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

* Basu, S. 1998, ApJ, 509, 229

* Bate, M. R. 2018, MNRAS, 475, 5618

* Bate, M. R., Lodato, G., & Pringle, J. E. 2010, MNRAS, 401, 1505

* Bjerkeli, P., Jørgensen, J. K., & Brinch, C. 2016, A&A, 587, A145

* Caselli, P., Benson, P. J., Myers, P. C., & Tafalla, M. 2002, ApJ, 572, 238

* Chen, C.-Y., Storm, S., Li, Z.-Y., et al. 2019, MNRAS, 490, 527

* Codella, C., Bachiller, R., & Reipurth, B. 1999, A&A, 343, 585

* Dunham, M. M., Allen, L. E., Evans, N. J., II., et al. 2015, ApJS, 220, 11

* Dunham, M. M., Stutz, A. M., Allen, L. E., et al. 2014, Protostars and Planets VI (Tucson, AZ: Univ. Arizona Press), 195

* Dunham, M. M., & Vorobyov, E. I. 2012, ApJ, 747, 52

* Dutta, S., Lee, C.-F., Hirano, N., et al. 2022, ApJ, 931, 130

* Evans, N. J., II., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321

* Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

* Galametz, M., Maury, A., Girart, J. M., et al. 2018, A&A, 616, A139

* Galametz, M., Maury, A., Girart, J. M., et al. 2020, A&A, 644, A47

* Gaudel, M., Maury, A. J., Belloche, A., et al. 2020, A&A, 637, A92
Gupta, A., Yen, H.-W., Koch, P., et al. 2022, ApJ, 930, 67

Hennebelle, P., Commerçon, B., Lee, Y.-N., & Charnoz, S. 2020, A&A, 635, A67

Hirano, S., & Machida, M. N. 2019, MNRAS, 485, 4667

Hull, C. L. H., Plambeck, R. L., Kwon, W., et al. 2014, ApJS, 213, 13

Hunter, J. D. 2007, CSE, 9, 90

Johnstone, D., Rosolowsky, E., Tafalla, M., & Kirk, H. 2010, ApJ, 711, 655

Joos, M., Hennebelle, P., & Ciardi, A. 2012, A&A, 543, A128

Kauffmann, J., Bertoldi, F., Bourke, T. L., Evans, N. J., II., & Lee, C. W. 2008, A&A, 487, 993

Kenyon, S. J., Hartmann, L. W., Strom, K. M., & Strom, S. E. 1990, AJ, 99, 869

Kirk, H., Johnstone, D., & Di Francesco, J. 2006, ApJ, 646, 1009

Kirk, H., Johnstone, D., & Tafalla, M. 2007, ApJ, 668, 1042

Lee, C.-F., Mundy, L. G., Reipurth, B., Ostriker, E. C., & Stone, J. M. 2000, ApJ, 542, 925

Lee, K. I., Dunham, M. M., Myers, P. C., et al. 2015, ApJ, 814, 114

Lim, J., Hanawa, T., Yeung, P. K. H., et al. 2016, ApJ, 831, 90

Machida, M. N., Hirano, S., & Kitta, H. 2020, MNRAS, 491, 2180

Machida, M. N., & Matsumoto, T. 2011, MNRAS, 413, 2767

Machida, M. N., & Matsumoto, T. 2012, MNRAS, 421, 588

Matsumoto, T., Machida, M. N., & Inutsuka, S.-i. 2017, ApJ, 839, 69

Matsumoto, T., & Tomisaka, K. 2004, ApJ, 616, 266

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Mellon, R. R., & Li, Z.-Y. 2008, ApJ, 681, 1356

Offner, S. S. R., & Chaban, J. 2017, ApJ, 847, 104

Offner, S. S. R., Lee, E. J., Goodman, A. A., & Arce, H. 2011, ApJ, 743, 91

Ohashi, N., Tobin, J. J., Jørgensen, J. K., et al. 2023, ApJ, 951, 8

Okoda, Y., Oya, Y., Francis, L., et al. 2021, ApJ, 910, 11

Okoda, Y., Oya, Y., Sakai, N., et al. 2018, ApJL, 864, L25

Oliphant, T. E. 2006, A Guide to NumPy (USA: Trelgol Publishing)

Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943

Oya, Y., Sakai, N., Sakai, T., et al. 2014, ApJ, 795, 152

Redaelli, E., Alves, F. O., Santos, F. P., & Caselli, P. 2019, A&A, 631, A154

Sai, J., Ohashi, N., Yen, H.-W., Maury, A. J., & Maret, S. 2023, ApJ, 944, 222

Santamaría-Miranda, A., de Gregorio-Monsalvo, I., Plunkett, A. L., et al. 2021, A&A, 646, A10

Sato, A., Tokuda, K., Machida, M. N., et al. 2023, ApJ, 958, 102

Seifried, D., Banerjee, R., Pudritz, R. E., & Klessen, R. S. 2013, MNRAS, 432, 3320

Shu, F. H., Ruden, S. P., Lada, C. J., & Lizano, S. 1991, ApJL, 370, L31

Takahashi, S. Z., Tomida, K., Machida, M. N., & Inutsuka, S.-i. 2016, MNRAS, 463, 1390

Terebey, S., Shu, F. H., & Cassen, P. 1984, ApJ, 286, 529

Thieme, T. J., Lai, S.-P., Ohashi, N., et al. 2023, ApJ, 958, 60

Tobin, J. J., Kratter, K. M., Persson, M. V., et al. 2016, Natur, 538, 483

Tobin, J. J., Sheehan, P. D., Megeath, S. T., et al. 2020, ApJ, 890, 130

Tomida, K., Okuzumi, S., & Machida, M. N. 2015, ApJ, 801, 117

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

van Kempen, T. A., van Dishoeck, E. F., Hogerheijde, M. R., & Güsten, R. 2009, A&A, 508, 259

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Ward-Thompson, D., André, P., Crutcher, R., et al. 2007, Protostars and Planets V (Tucson, AZ: Univ. Arizona Press), 33

Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67

Yen, H.-W., Koch, P. M., Hull, C. L. H., et al. 2021b, ApJ, 907, 33

Yen, H.-W., Koch, P. M., Takakuwa, S., et al. 2017, ApJ, 834, 178

Yen, H.-W. & the eDisk team. 2024, ApJ, submitted

Yen, H.-W., Zhao, B., Koch, P. M., & Gupta, A. 2021a, ApJ, 916, 97

Yıldız, U. A., Kristensen, L. E., van Dishoeck, E. F., et al. 2012, A&A, 542, A86

Yıldız, U. A., Kristensen, L. E., van Dishoeck, E. F., et al. 2013, A&A, 556, A89

Yıldız, U. A., Kristensen, L. E., van Dishoeck, E. F., et al. 2015, A&A, 576, A109

Zucker, C., Speagle, J. S., Schlafly, E. F., et al. 2020, A&A, 633, A51