# Which Part of Dense Cores Feeds Material to Protostars? The Case of L1489 IRS

Jinshi Sai, Nagayoshi Ohashi, Anaëlle J. Maury, Sébastien Maret, Hsi-Wei Yen, Yusuke Aso, Mathilde Gaudel

► To cite this version:

Jinshi Sai, Nagayoshi Ohashi, Anaëlle J. Maury, Sébastien Maret, Hsi-Wei Yen, et al. Which Part of Dense Cores Feeds Material to Protostars? The Case of L1489 IRS. The Astrophysical Journal, 2022, 925, 10.3847/1538-4357/ac341d. insu-03705394

HAL Id: insu-03705394

https://insu.hal.science/insu-03705394

Submitted on 27 Jun 2022

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
# Which Part of Dense Cores Feeds Material to Protostars? The Case of L1489 IRS

Jinshi Sai (Insa Choi)<fnref n="1" />,<fnref n="2" />, ID, Nagayoshi Ohashi<fnref n="2" />, ID, Anaëlle J. Maury<fnref n="3" />,<fnref n="4" />, ID, Sébastien Maret<fnref n="5" />, ID, Hsi-Wei Yen<fnref n="2" />, ID, Yusuke Aso<fnref n="6" />, ID, and Mathilde Gaudel<fnref n="7" />, ID

## Abstract

We have conducted mapping observations (~2' × 2') of the Class I protostar L1489 IRS using the 7 m array of the Atacama Compact Array and the IRAM 30 m telescope in C¹⁸O 2–1 emission to investigate the gas kinematics on 1000–10,000 au scales. The C¹⁸O emission shows a velocity gradient across the protostar in a direction almost perpendicular to the outflow. The radial profile of the peak velocity was measured from a C¹⁸O position–velocity diagram cut along the disk major axis. The measured peak velocity decreases with radius at radii of ~1400–2900 au, but increases slightly or is almost constant at radii of $r \gtrsim 2900$ au. Disk-and-envelope models were compared with the observations to understand the nature of the radial profile of the peak velocity. The measured peak velocities are best explained by a model where the specific angular momentum is constant within a radius of 2900 au but increases with radius outside 2900 au. We calculated the radial profile of the specific angular momentum from the measured peak velocities and compared it to analytic models of core collapse. The analytic models reproduce well the observed radial profile of the specific angular momentum and suggest that material within a radius of ~4000–6000 au in the initial dense core has accreted to the central protostar. Because dense cores are typically ~10,000–20,000 au in radius, and as L1489 IRS is close to the end of its mass accretion phase, our result suggests that only a fraction of a dense core eventually forms a star.

*Unified Astronomy Thesaurus concepts: Star formation (1569); Protostars (1302); Low mass stars (2050); Circumstellar disks (235); Radio astronomy (1338); Interstellar medium (847); Submillimeter astronomy (1647); Millimeter astronomy (1061)*

## 1. Introduction

Low-mass stars are formed through the gravitational collapse of dense cores, which typically have a radius of ~0.05–0.1 pc (or 10,000–20,000 au) and mass of ~0.5 $M_\odot$ (di Francesco et al. 2007; André et al. 2014; Könyves et al. 2015; Pokhrel et al. 2018). Comparison between the initial mass function (IMF) and core mass functions (CMFs) of starless cores suggests that only ~40% of the mass in each dense core contributes to stellar mass (André et al. 2010, 2014; Könyves et al. 2015). If correct, the similarity of the IMF and CMFs implies that most of the material contained in dense cores does not end up in the resulting star(s). However, observationally, it remains unclear that which part of the dense core is eventually accreted onto the central protostar(s), preventing the understanding of star formation efficiency in individual cores. In this paper, we define a zone within the initial core, wherein material forming protostars resides, as the “stellar mass feeding zone” (SMFZ) and investigate it.

A way to estimate SMFZs around protostars is to measure the radial dependence of the specific angular momentum; infalling material conserving its specific angular momentum can be traced back to its original location within the initial

dense core as long as the distribution of its specific angular momentum within the initial dense core is known.

The radial dependence of the specific angular momentum within dense cores has been statistically studied based on the relation between the core size and the mean specific angular momentum $J/M$, which is well fitted by a power law $J/M \propto r^{1.6}$ within a radius range of ~0.03–0.3 pc (Goodman et al. 1993; Caselli et al. 2002a; Tatematsu et al. 2016). On the other hand, specific angular momenta $j$ were measured in several infalling envelopes at radii of 200–2000 au, and it was found that these envelopes show similar specific angular momenta of ~$10^{-3}$ km s$^{-1}$ pc. This finding was interpreted as a consequence of angular momentum conservation during the infalling process (e.g., Ohashi et al. 1997). Ohashi et al. (1997) have estimated the typical size scale of the SMFZ to be ~6000 au by comparing a heterogeneous ensemble of measurements of the specific angular momentum on different spatial scales (see also Belloche 2013).

Recent observations at higher spatial resolutions have allowed us to directly measure radial profiles of rotational velocity, i.e., specific angular momentum in individual protostellar systems. It has been reported that the specific angular momentum around protostars is often constant within radii of ~100–1000 au (Lee 2010; Murillo et al. 2013; Yen et al. 2013, 2017; Harsono et al. 2014; Ohashi et al. 2014; Aso et al. 2015, 2017; Maret et al. 2020). These results are basically consistent with the finding made by earlier studies described above, although these results also suggested that older systems

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. 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">¹Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan</footnote>

<footnote n="2">²Academia Sinica Institute of Astronomy and Astrophysics, 11F of Astro-Math Bldg, 1, Section 4, Roosevelt Rd, Taipei 10617, Taiwan; ohashi@asiaa.sinica.edu.tw</footnote>

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

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

<footnote n="5">⁵Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="6">⁶Korea Astronomy and Space Science Institute (KASI), 77 Daedokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea</footnote>

<footnote n="7">⁷LERMA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, 75014 Paris, France</footnote>

<footnote>Received 2021 May 6; revised 2021 October 25; accepted 2021 October 26; published 2022 January 20</footnote>
have larger constant specific angular momenta (Yen et al. 2017). On the other hand, radial profiles of the specific angular momentum were measured at larger radii of 800–10,000 au in three protostellar and prestellar objects and found to follow $j \propto r^{1.8}$ (Pineda et al. 2019). This radial dependence of the specific angular momentum is very similar to the relation between the core size and the mean specific angular momentum of dense cores found by Goodman et al. (1993).

These previous studies have consistently demonstrated that the radial distributions of the specific angular momentum in dense cores and infalling envelopes are divided into two regimes: a $j$-constant regime at inner radii and a $j$-increasing regime at outer radii. Gaudel et al. (2020), who measured specific angular momentum distributions around 12 Class 0 sources over a wide range of the radius (~50–5000 au), detected transitions of the two regimes for first time. After taking the average of the measured specific angular momentum profiles, Gaudel et al. showed that the averaged profile is divided into two regimes at a radius of ~1600 au and the $j$-increasing regime is described by $j \propto r^{1.6}$.

Gaudel et al. (2020) have clearly demonstrated that measurements of the specific angular momentum distribution over a wider range of the radius enable us to find the transition between the two regimes of the specific angular momentum distribution. Note that the specific angular momentum distribution in the $j$-increasing regime helps us to constrain the specific angular momentum distributions of the initial dense cores, which are required to estimate SMFZs. Gaudel et al. (2020) have also shown that specific angular momentum distributions seem to vary from source to source to some degree, suggesting the importance of the measurements for individual sources.

In this paper, we have examined an SMFZ around the Class I protostar L1489 IRS as a case study. L1489 IRS is relatively isolated in the Taurus molecular cloud ($d \sim 140$ pc; Zucker et al. 2019) and located close to the western edge of the dense core L1489 (Benson & Myers 1989; Hogerheijde & Sandell 2000; Motte & André 2001; Wu et al. 2019). The object is considered to be a late Class I protostar according to its bolometric temperature and luminosity of 226 K and 3.5 $L_\odot$, respectively (Green et al. 2013). Hence, it is a good target to examine an SMFZ. Observations at (sub)millimeter wavelengths reported a dusty envelope around the protostar, which has a radius of ~2000 au and mass of ~0.02–0.03 $M_\odot$ (Hogerheijde & Sandell 2000; Motte & André 2001), and rotational motion of the envelope (Hogerheijde 2001; Yen et al. 2013). A bipolar outflow is seen in the north–south direction spanning scales of thousands of au (Tamura et al. 1991; Hogerheijde et al. 1998; Yen et al. 2014). Yen et al. (2014) have identified a Keplerian disk with a radius of ~700 au, which is the largest one identified around a Class I protostar so far, and two infalling flows accreting onto the disk surface based on their Atacama Large Millimeter/submillimeter Array (ALMA) 12 m array observations in $C^{18}\text{O}$ 2–1 emission at angular resolutions of ~1". Sai et al. (2020) have re-examined the gas kinematics at an angular resolution about three times higher than that in Yen et al. (2014), revealing that a Keplerian disk with a radius of ~600 au is surrounded by an infalling envelope with a constant specific angular momentum. The dynamical mass of L1489 IRS was also estimated to be ~1.6 $M_\odot$ (Yen et al. 2014; Sai et al. 2020), which is larger than that of other protostars whose dynamical masses have been estimated (Murillo et al. 2013; Harsono et al. 2014;

Chou et al. 2014; Lee et al. 2014; Aso et al. 2015, 2017). The systemic velocity of L1489 IRS was estimated from the Keplerian rotation to be 7.22 km s⁻¹ (Sai et al. 2020). We have adopted $v_{LSR} = 7.22$ km s⁻¹ for the systemic velocity of L1489 IRS in this paper.

The outline of this paper is as follows. Observations and data reduction are summarized in Section 2. Observational results are presented in Section 3. In Section 4, velocity structures of the $C^{18}\text{O}$ 2–1 emission are examined in detail by measuring the peak velocity as a function of radius and comparing it to that measured with kinematic models. We then derive a radial profile of the specific angular momentum of the envelope, discuss its nature, and examine the SMFZ around L1489 IRS in Section 5. We have also discussed an infalling velocity that is lower than the freefall velocity expected from the stellar mass of L1489 IRS. Finally, all results and discussions are summarized in Section 6.

## 2. Observations and Data Reduction

### 2.1. IRAM 30 m Observations

We carried out observations using the IRAM 30 m telescope from 2019 April 9 to April 14. Two molecular lines, $C^{18}\text{O}$ $J=2-1$ and $N_2H^+$ $J=1-0$, were observed with the heterodyne Eight Mixer Receivers E230 and E090 in the 1.3 mm and 3 mm atmospheric windows, respectively. The rest frequencies are 219.560358 GHz for $C^{18}\text{O}$ 2–1 and 93.173770 GHz for $N_2H^+$ 1–0. The $N_2H^+$ 1–0 line is a multiplet with seven hyperfine structure (HFS) components (see, e.g., Caselli et al. 1995). Here we provide the rest frequency of the brightest HFS component, $N_2H^+$ $J_{F_i,F} = 1_2,3-0_1,2$. The Versatile Spectrometer Arrays backend was connected to the E230 and E090 receivers, providing spectral resolutions of 6.5 kHz (0.0089 km s⁻¹) for $C^{18}\text{O}$ 2–1 and 20 kHz (0.063 km s⁻¹) for $N_2H^+$ 1–0. The velocity resolution for $C^{18}\text{O}$ 2–1 was smoothed to 0.17 km s⁻¹ to increase the signal-to-noise ratio (S/N) and to be combined with Atacama Compact Array (ACA) data. The observations were carried out in the on-the-fly mapping mode with position-switching, with a reference position located at $\Delta\alpha = -10'$, $\Delta\delta = -10'$ with respect to the map center of $\alpha(J2000) = 4^h4^m43s.07$, $\delta(J2000) = +26^\circ18'56.2''$, which is a peak position of the 1.3 mm continuum emission in previous ALMA observations (Yen et al. 2014; Sai et al. 2020). We consider this position to be the protostellar position throughout this paper. The maps cover a ~2' × 2' region around the protostar. The atmospheric opacity at 225 GHz, $\tau_{225}$, was ~0.2 on average and ranged from 0.10 to 0.31. The telescope pointing and focus were corrected using Mars every 1–2 hr and ~4 hr, respectively. The obtained data set was reduced with the GILDAS software package.<fnref n="8" /> The antenna temperature $T_A^*$ was converted to the main beam temperature $T_{mb}$ using the values of $B_{\text{eff}} = 0.60$ and $F_{\text{eff}} = 0.92$ for $C^{18}\text{O}$ 2–1 and $B_{\text{eff}} = 0.81$ and $F_{\text{eff}} = 0.95$ for $N_2H^+$ 1–0, where $B_{\text{eff}}$ is the main beam efficiency, $F_{\text{eff}}$ is the forward efficiency, and $T_{mb} = T_A^* F_{\text{eff}} / B_{\text{eff}}$. The half-power beamwidth (HPBW) is 12'' for the $C^{18}\text{O}$ 2–1 line and 28'' for the $N_2H^+$ 1–0 line. The rms noise levels of the $C^{18}\text{O}$ 2–1 and $N_2H^+$ 1–0 images are 110 mK and 42 mK, respectively in $T_{mb}$. We mainly focus on the $C^{18}\text{O}$ emission in this paper. The map of the $N_2H^+$ emission is presented in Appendix B.

<footnote n="8">⁸ http://www.iram.fr/IRAMFR/GILDAS</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) The moment 0 map of the C¹⁸O 2–1 emission (background color) and the 1.3 mm NIKA2 continuum map (contours). Contour levels start at 3σ and increase in steps of 2σ, where 1σ is 5 mJy beam⁻¹ km s⁻¹. The central and northern crosses denote the protostellar and C¹⁸O second-peak positions, respectively. The filled ellipse and the dashed circle in the bottom-left corner show the beam size of the C¹⁸O and continuum maps, respectively. (b) The moment I map of the C¹⁸O 2–1 emission (background color) overlaid with the moment 0 map of the C¹⁸O 2–1 emission (contours). Contour levels start at 5σ and increase in steps of 5σ between 5σ and 20σ, steps of 10σ between 20σ and 60σ, and steps of 20σ between 60σ and 100σ, where 1σ is 0.12 mJy beam⁻¹ km s⁻¹. Blue and red arrows denote the direction of the blueshifted and redshifted components of the outflow, respectively (Hogherheijde et al. 1998). The filled ellipse in the bottom-left corner shows the beam size of the maps. (c) The same as (b) but the moment II map of the C¹⁸O 2–1 emission shown in the background color.</figcaption></figure>

We have also mapped a 10' square region around the protostar in the 1.3 mm continuum with the millimeter camera NIKA2 (Adam et al. 2018; Perotto et al. 2020) on 2019 February 4 in the pool observations. The observing time was ~0.4 hr and $\tau_{225}$ ranged from 0.24 to 0.31. The obtained data were reduced in the pipeline for NIKA2 data reduction. We then smoothed the obtained map to increase S/N. The final HPBW is 20'', and the rms measured at the off-point in the map is 5 mJy beam⁻¹.

## 2.2. ACA 7 m array Mosaic Observations

Observations using the 7 m array of the ACA were conducted on five nights during 2019 December 1–17 in the ALMA Cycle 7 phase with ten antennas. The baseline length ranged from 8.9 to 45.0 m, and the minimum baseline provided sensitivity to structures extending to the 31'' scale at a level of 10% (Wilner & Welch 1994). The spectral setup of the observations consisted of five spectral windows for the 1.3 mm continuum and four molecular lines of C¹⁸O 2–1, ¹³CO 2–1 (220.398684 GHz), ¹²CO 2–1 (230.538000 GHz), and N₂D⁺ 3–2 (231.321828 GHz). The spectral windows for the C¹⁸O and ¹³CO 2–1 lines, placed within individual basebands, had a bandwidth of 120.0 MHz and a spectral resolution of 61.0 kHz. The spectral windows for the ¹²CO 2–1 and N₂D⁺ 3–2 lines, placed within a single baseband, had a bandwidth of 62.5 MHz and a spectral resolution of 61.0 kHz. The spectral window for the 1.3 mm continuum with a bandwidth of 2 GHz was placed within the other baseband. In this paper, we focus on the C¹⁸O 2–1 line. The elevation of the target source was ~40° on average. The data set was calibrated in the pipeline with the Common Astronomy Software Applications package (CASA, McMullin et al. 2007) version 5.6.1 and its pipeline version 42866M (Pipeline-CASA56-P1-B). The calibration sources were the quasar J0426 + 2327 for the phase calibration and the quasars J0423-0120 and J0725-0054 for the bandpass and flux calibration. Mosaic observations consisted of 27 fields and covered a ~2' × 2' region. The on-source time per pointing was ~2 minutes.

The C¹⁸O 2–1 map was produced with CASA 5.6.1 using the *tclean* task. The frequency channels were smoothed to achieve a higher S/N and the resultant velocity resolution was set at 0.17 km s⁻¹. We adopted Briggs weighting with a robust

value of 0.5, the Hogbom deconvolver, and CLEAN masks interactively drawn. The synthesized beam size is ~7''7 × 6''4 (~85°) and the rms noise level is ~83 mJy beam⁻¹ (~43 mK).

## 2.3. Combination of the ACA and IRAM 30 m Data

The ACA 7 m array and IRAM 30 m images of the C¹⁸O emission, which covered similar regions (~2' × 2'), were combined for further analysis. Before the images were combined, consistency of the flux scales of the ACA 7 m array and IRAM 30 m data was confirmed by comparing the total flux of the two data sets in both the image and the uv domain. For this comparison of the flux scales, pseudo-visibility of the IRAM 30 m data was produced through the Fourier transform of the IRAM 30 m image with the task *uvshort* in CLASS of the GILDAS software.

The two images were combined with the *feather* task in CASA 5.6.1. This task combines interferometer and single-dish images in the uv domain through Fourier transformation, weighting them by the spatial frequency response of the images.<fnref n="9" /> The detailed processes of combining images are as follows: (1) the IRAM 30 m and ACA 7 m array maps were trimmed to exclude masked and noisy regions at the map edge, (2) the IRAM 30 m map was multiplied by the primary-beam response of the ACA 7 m mosaic, (3) the primary-beam-modified IRAM 30 m map and the ACA 7 m array map were combined with the *feather* task, and (4) the combined map was corrected by the primary-beam response of the ACA 7 m mosaic. The synthesized beam size and rms of the combined image are 7''7 × 6''4 (~85°) and ~0.14 Jy beam⁻¹ (~70 mK), respectively. We used the combined image for all the analysis of the C¹⁸O emission in the following sections.

# 3. Results

## 3.1. Overall Structures of Emissions

Figure 1(a) shows the 1.3 mm NIKA2 continuum map and the velocity-integrated intensity (moment 0) map of the C¹⁸O emission. The moment 0 map was calculated by integrating the

<footnote n="9">⁹ See https://casa.nrao.edu/casadocs/casa-5.4.1/image-combination/feather for more detail.</footnote>
C¹⁸O emission over the velocity range 5.09–9.34 km s⁻¹, within which the emission was detected above at least the 3σ level at the protostellar position. The 1.3 mm continuum emission shows a compact structure with a radius of ~15″ (~2100 au), which is resolved with the smoothed NIK2A beam, at the protostellar position. This traces a dense, dusty envelope around the protostar. The continuum emission is extended to the northeast from the center and shows the second peak to the east of the protostar. This second peak originates within a starless core near the protostar, as reported in previous observations (Hogerheijde & Sandell 2000; Motte & André 2001; Wu et al. 2019). The C¹⁸O emission exhibits a peak at the protostellar position and it thus likely traces an envelope associated with the protostar. The second peak appears at a position ~50″ northeast of the protostar, where the 1.3 mm continuum emission does not present a peak.

## 3.2. Velocity Structure of the C¹⁸O 2–1 Emission

Figure 1(b) provides the velocity-integrated intensity and mean velocity (moment 0 and I) maps of the C¹⁸O emission. A velocity gradient is observed from the northeast toward the southwest across the protostellar position. This direction is almost perpendicular to the outflow direction and parallel to that of the velocity gradient due to the rotation of the disk and the envelope seen at radii of ~100–1000 au (Yen et al. 2014; Sai et al. 2020), suggesting that the velocity gradient of ≥1000 au scales is also due to rotational motion. The velocity gradient appears steeper in the vicinity of the protostar and gentler on the outside. To confirm this, the magnitude of the velocity gradient was measured across different scales in the same manner as Goodman et al. (1993) did. We fitted the following function to the mean velocities within a radius of r<sub>fit</sub>:

$$v_{\text{LSR}} = v_0 + a\Delta\alpha + b\Delta\delta, \quad (1)$$

where $v_0$ is the systemic velocity, $\Delta\alpha$ and $\Delta\delta$ are offsets from the protostellar position in R.A. $\alpha$ and decl. $\delta$, and a and b are the magnitude of the velocity gradient in the $\alpha$ and $\delta$ directions, respectively. The magnitude and position angle of the velocity gradient are calculated as follows (Maret et al. 2020):

$$G = (a^2 + b^2)^{1/2}/d, \quad (2)$$

$$\theta = \tan^{-1}(a/b), \quad (3)$$

where G is the magnitude, $\theta$ is the position angle measured from north to east, and d is the distance of the source. Table 1 lists r<sub>fit</sub>, G, and $\theta$, clearly showing that the magnitude of the velocity gradient is larger in inner regions.

The velocity dispersion (or moment II) map presented in Figure 1(c) shows that the velocity dispersion decreases rapidly with radius from ~0.7 km s⁻¹ to ~0.2 km s⁻¹ within a radius of ~10″–20″ around the protostellar position, while it is almost constant at larger radii.

The velocity structure of the C¹⁸O emission is shown in more detail in the velocity channel maps presented in Figure 2. As shown in Figure 1(b), a clear velocity gradient is seen from the northeast to southwest. The emission is compact and within r<~10″ around the protostar at higher blueshifted ($v_{LSR} \le 6.28$ km s⁻¹) and higher redshifted ($v_{LSR} \ge 8.66$ km s⁻¹) velocities, although there is a weak extension to the northeast at 6.28 km s⁻¹. Keplerian rotation was found at these velocities in previous ALMA observations (Yen et al. 2014; Sai et al. 2020).

<table><caption>**Table 1**
The Magnitude and Position Angle of the Velocity Gradient Measured with the C¹⁸O Emission</caption><thead><tr><th>r<sub>fit</sub><br>(arcsec)</th><th>G<br>(km s<sup>-1</sup> pc<sup>-1</sup>)</th><th>θ<br>(deg)</th></tr></thead><tbody><tr><td>10</td><td>89 ± 1</td><td>-121.8 ± 0.7</td></tr><tr><td>20</td><td>27.3 ± 0.2</td><td>-125.2 ± 0.4</td></tr><tr><td>30</td><td>14.79 ± 0.09</td><td>-122.0 ± 0.3</td></tr><tr><td>40</td><td>9.58 ± 0.05</td><td>-115.2 ± 0.3</td></tr><tr><td>50</td><td>6.98 ± 0.03</td><td>-110.9 ± 0.2</td></tr><tr><td>60</td><td>5.46 ± 0.02</td><td>-111.7 ± 0.2</td></tr></tbody></table>

Thus, the Keplerian disk component is dominant in the compact emission at these velocities. At lower redshifted velocities (7.64 km s⁻¹ ≤ $v_{LSR}$ ≤ 8.49 km s⁻¹), the emission structures are also relatively compact but with a slight extension of r ~ 10″–30″ in the southwest direction, even though there is also a weak extension from the northwest to southeast at 7.64 km s⁻¹. Our previous work has found rotational motion of the infalling envelope around L1489 IRS at these velocities on a scale of r ~ 1000 au ~ 7″ (Sai et al. 2020). These results show that the redshifted emissions trace either the disk or infalling envelope.

On the other hand, the emission structures at lower blueshifted velocities (6.45 km s⁻¹ ≤ $v_{LSR}$ ≤ 7.47 km s⁻¹) are extended across the entire maps. The extended emissions at LSR velocities of 6.62–6.96 km s⁻¹ appear to be associated with the C¹⁸O secondary peak seen in the moment 0 map because they show stronger peaks at the C¹⁸O secondary peak position. The molecular emissions associated with the L1489 starless core show their peaks at LSR velocities of ~6.6–6.8 km s⁻¹ (Caselli et al. 2002a; Wu et al. 2019), suggesting that the extended emission is also a part of the starless core.

Figure 3 shows a position–velocity (PV) diagram cut along the direction of the disk major axis (position angle P.A. = 54°, Sai et al. 2020) to investigate the velocity gradient of the C¹⁸O emission in more detail. It is clear that the velocity of the redshifted emission increases as the position approaches the protostellar position, indicative of differential rotation. The velocity structure of the redshifted emission is mostly consistent with the power-law profile $r^{-1}$ of the C¹⁸O infalling envelope measured at radii of ~600–1000 au by Sai et al. (2020), shown by the black curves in Figure 3, although a part of the emission at $v_{LSR} > 8.66$ km s⁻¹ arises from the Keplerian disk, as was mentioned above. In order to understand how Keplerian and infalling motions contribute to the redshifted emission in the PV diagram, further careful analysis of the velocity structures is required (see Section 4.2).

The velocity structures of the blueshifted emission, on the other hand, appear different from those of the redshifted emission. They are consistent with the $r^{-1}$ profile of the infalling envelope measured in Sai et al. (2020) at velocities of $v_{LSR} < 6$ km s⁻¹, showing a feature of differential rotation. At velocities of $v_{LSR} > 6$ km s⁻¹, however, the velocity structures do not follow the $r^{-1}$ profile of the infalling envelope, and the peak velocities are ~0.5 km s⁻¹ smaller than those expected from the $r^{-1}$ profile. This inconsistency is probably due to the extended emission seen at $v_{LSR} = 6.45–7.47$ km s⁻¹ in the channel maps.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Velocity channel maps of the C$^{18}$O 2–1 emission. Contour levels are 3, 6, 12, 24, and 48 × σ, where 1σ = 0.14 Jy beam$^{-1}$. The labels in the top-right corner show the LSR velocity of each channel in km s$^{-1}$. The central and northern crosses denote the protostellar position and the C$^{18}$O second-peak position measured in the moment 0 map, respectively. The black filled ellipse in the bottom-left corner denotes the beam size.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Position–velocity diagram of the C$^{18}$O 2–1 emission cut along the direction of the disk major axis (P.A. = 54°). Contour levels are 3, 6, 12, and 24 × σ, where 1σ = 0.14 Jy beam$^{-1}$. The black curves indicate the power-law profile r$^{-1}$ of the C$^{18}$O infalling envelope measured at radii of ~600–1000 au in a previous work (Sai et al. 2020). The vertical and horizontal bars in the bottom-left corner denote the FWHM of the beam major axis and the velocity resolution, respectively.</figcaption></figure>

# 4. Analysis

## 4.1. Measurement of the Peak Velocity of the C$^{18}$O 2–1 Emission

In order to characterize the velocity gradient seen at the redshifted velocity in the PV diagram in Figure 3 more quantitatively, the peak velocity was measured as a function of radius in the diagram in the same manner as performed by Sai et al. (2020). Similar methods were also used in other works (Yen et al. 2013; Ohashi et al. 2014; Aso et al. 2015, 2017; Maret et al. 2020). The measurement was only made for the redshifted component because the blueshifted component is highly contaminated by the extended gas associated with the starless core.

A Gaussian function was fitted to the spectrum at each position in the PV diagram to measure peak velocities. The fitting starts from the offset of −10" because velocity structures are not spatially well-resolved near the protostar due to the limited angular resolution. The sampling step is 3"/2, which corresponds to half of the beam size. The fitting was performed at positions where the emission was detected at levels above 6σ, resulting in the offset of the outermost point = −35". Velocity channels within ±2 channels around the channel with the maximum intensity were used for the fitting. This channel selection method traces the peak velocity more accurately than using all the velocity channels when a spectrum is not of a Gaussian shape (see Appendix B in Sai et al. 2020). The uncertainties of the representative data points along position and velocity axes were assumed to be the position accuracy, which is given as angular resolution/(S/N), and fitting error, respectively. Figure 4(a) shows representative data points on the PV diagram.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) Same as Figure 3 but with zoomed-in offset and velocity ranges of −40"−0" and 7−10 km s⁻¹, respectively. Black points represent peak velocities measured through Gaussian fittings to the PV diagram. The vertical and horizontal bars in the bottom-right corner denote the FWHM of the beam major axis and the velocity resolution, respectively. (b) Plots of the data points measured in the PV diagram on a log r−log v plane. The solid line shows the power-law function with the best-fit parameters.</figcaption></figure>

These data points are plotted on a log r−log v plane, as presented in Figure 4(b). The measured peak velocity decreases with radius within ~3000 au, while it slightly increases or is almost constant outside ~3000 au. To characterize the two different behaviors, we fitted the following double power-law function to the data points through χ² fitting:

$$v_{\text{peak}} = \begin{cases} v_{\text{break}} \left( \frac{r}{r_{\text{break}}} \right)^{p_{\text{in}}} & (r \le r_{\text{break}}) \\ v_{\text{break}} \left( \frac{r}{r_{\text{break}}} \right)^{p_{\text{out}}} & (r > r_{\text{break}}) \end{cases}, \quad (4)$$

where $v_{\text{peak}} = |v_{\text{LSR}} - v_{\text{sys}}|$. We evaluated uncertainties of the fitting parameters with the Monte Carlo method by iterating the fitting procedure 3000 times, changing the offset and velocity values of the data points within their errors. The uncertainties of the fitting parameters were obtained as the standard deviation of the Gaussian distribution fitted to the posterior distributions. In this procedure, uniform and Gaussian prior distributions were assumed for the data errors along the velocity and offset axes, respectively.

A solid line in Figure 4(b) shows the function with the best-fit parameters of ($v_{\text{break}}$, $r_{\text{break}}$, $p_{\text{in}}$, $p_{\text{out}}$) = (0.29 ± 0.02 km s⁻¹, 2900 ± 200 au, −1.2 ± 0.1, 0.3 ± 0.2). The best-fit power-law index at radii smaller than the break radius is close to −1. This suggests that the peak velocity mostly traces rotating motion with a constant specific angular momentum, although further careful analysis is required because part of the infalling motion could cause contamination (see the next subsection). On the other hand, the derived power-law index at radii larger than 2900 au is 0.3 ± 0.2, indicating that the velocity slightly increases or is almost constant with radius.

## 4.2. Disk-and-envelope Model

In order for us to interpret the origin of the velocity gradient seen in the PV diagram, disk-and-envelope models are constructed and compared with observations in this section.

Although several disk-and-envelope models are examined to find the one explaining observations best in the following discussions, our base model is based on the envelope model proposed by Ulrich (1976) and the disk model discussed in Sai et al. (2020). The envelope model proposed by Ulrich (1976) is spherical and consists of trajectories of material freefalling from infinite radii and conserving angular momentum. The assumption of the spherical envelope is plausible because the integrated intensity of the C¹⁸O emission shows a spherical shape at the protostellar position, as seen in Figure 1. The envelope model is symmetric about the z-axis. The velocity and density profiles are given as follows in spherical coordinates:

$$\rho(r, \theta) = \rho_0 \left(\frac{r}{r_c}\right)^{-1.5} \left(1 + \frac{\cos \theta}{\cos \theta_0}\right)^{-0.5} \times \left(\frac{\cos \theta}{2 \cos \theta_0} + \frac{r_c}{r} \cos^2 \theta_0\right)^{-1}, \quad (5)$$

$$v_r(r, \theta) = - \left( \frac{GM_*}{r} \right)^{0.5} \left( 1 + \frac{\cos \theta}{\cos \theta_0} \right)^{0.5}, \quad (6)$$

$$v_{\theta}(r, \theta) = \left(\frac{GM_*}{r}\right)^{0.5} (\cos \theta_0 - \cos \theta) \times \left(\frac{\cos \theta_0 + \cos \theta}{\cos \theta_0 \sin^2 \theta}\right)^{0.5}, \quad (7)$$

$$v_{\phi}(r, \theta) = \left(\frac{GM_*}{r}\right)^{0.5} \frac{\sin \theta_0}{\sin \theta} \left(1 - \frac{\cos \theta}{\cos \theta_0}\right)^{0.5}. \quad (8)$$

Here G is the gravitational constant, $M_*$ is the stellar mass, $\theta_0$ is the initial location of the infalling material in $\theta$ at an infinite radius, and $r_c = j^2/GM_*$ is the centrifugal radius, where j is the specific angular momentum. The reference density $\rho_0$ is defined as $\rho = 2\rho_0/3\sqrt{2}$ at $(r, \theta) = (r_c, 0^\circ)$.

Ulrich's envelope model has five parameters—$M_*$, $r_c$, $r_{\text{env}}$, i, and $\rho_0$—where $r_{\text{env}}$ is the envelope radius and i is the inclination angle of the system. First, we adopted the following values for the parameters based on observations for a base
<table><caption>**Table 2**
Parameters of the Base Disk-and-envelope Model</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td colspan="2"><em>Protostellar and global parameters</em></td></tr><tr><td>Protostellar mass (M<sub>*</sub>)</td><td>1.6 M<sub>☉</sub></td></tr><tr><td>C<sup>18</sup>O abundance (X<sub>C<sup>18</sup>O</sub>)</td><td>1.7 × 10<sup>-7</sup></td></tr><tr><td colspan="2"><em>Disk parameters</em></td></tr><tr><td>Disk mass (M<sub>disk</sub>)</td><td>0.0071 M<sub>☉</sub></td></tr><tr><td>Power of surface density profile (p)</td><td>0.5</td></tr><tr><td>Disk radius (R<sub>disk</sub>)</td><td>600 au</td></tr><tr><td>Position angle (P.A.)</td><td>54°</td></tr><tr><td colspan="2"><em>Envelope parameters</em></td></tr><tr><td>Centrifugal radius (r<sub>c</sub>)</td><td>600 au</td></tr><tr><td>Envelope radius (r<sub>env</sub>)</td><td>5000 au</td></tr><tr><td>Inclination angle (i)</td><td>73°</td></tr><tr><td>Density (ρ<sub>0</sub>)</td><td>1.4 × 10<sup>-18</sup> g cm<sup>-3</sup></td></tr></tbody></table>

model. For the protostellar mass, $M_* = 1.6M_☉$ derived from the Keplerian rotation of the disk was adopted (Yen et al. 2014; Sai et al. 2020). $r_c = R_{disk} = 600$ au was assumed for $r_c$, as measured by Sai et al. (2020), and $r_{env}$ was fixed at 5000 au based on the extent of the C<sup>18</sup>O emission. The inclination angle of the system was assumed to be 73°, according to the inclination angle of the disk estimated from the aspect ratio of the dust continuum emission (Sai et al. 2020). Although the inclination angle is uncertain, the model results do not strongly depend on it if it is greater than 60°, as is the case with L1489 IRS (Stark et al. 2006; Brinch et al. 2007; Eisner 2012; Sheehan & Eisner 2017). The volume density $\rho_0$ was initially fixed at $1.4 \times 10^{-18} \text{ g cm}^{-3}$, which provides an intensity of C<sup>18</sup>O almost consistent with the observations, while it was adjusted when our models were modified to explain the observations better.

Regarding the disk in our models, the disk model discussed in Sai et al. (2020) was adopted except for its warped and gap structures. The warped and gap structures affect only the morphology at $r < 600$ au and are thus negligible in comparison with the current observations with an angular resolution of ~1000 au. Note that the same inclination and position angles were assumed for both the disk and the envelope, although they do not necessarily have the same rotational axis as suggested in the past observations (Brinch et al. 2007; Sai et al. 2020). All the parameters for the base model are summarized in Table 2. Because the infall velocity is freefall in the base model, the base model is called the disk-and-envelope model with freefall in the following discussions.

In order to compare our models with the observations, synthetic observations were carried out based on our models. The intensity scales of the synthetic observations were based on the radiative transfer calculations using RADMC-3D<sup>10</sup> with the assumption of the LTE condition. The opacity table derived by Semenov et al. (2003)<sup>11</sup> was adopted. The protostellar temperature ~4000 K calculated from $3.5 L_{bol}$ was used for the protostellar spectral energy distribution to calculate the temperature profile within the envelope and the disk. The molecular abundance of C<sup>18</sup>O relative to H<sub>2</sub> gas and the gas-to-dust mass ratio were assumed to be $1.7 \times 10^{-7}$ and 100, respectively (Frerking et al. 1982). We added Gaussian noise to

the model images and convolved them with a Gaussian beam so that they have the same resolution and S/N as the observations. PV diagrams were generated from the calculated model images and compared to the observed PV diagram.

The first row in Figure 5 shows the comparison between the disk-and-envelope model with freefall and the observations. As discussed in Section 3, the blueshifted emission in the observations is not likely associated with L1489 IRS. Hence, comparisons were made only for the redshifted components. PV diagrams zoomed in to a velocity range of 7 to 10 km s<sup>-1</sup> and an offset range of -40" to 7" are presented in the second column in Figure 5 to compare the redshifted components between the observations and models. As seen in the second panel in the first row, the model and observed PV diagrams are not consistent. The model PV diagram shows intensity peaks at higher velocities than the observations, and is wider along the offset axis than in the observations. In order to make further comparisons, the peak velocity was measured in the model PV diagram in the same way as that used for the observed PV diagram in Section 4.1. The third column in Figure 5 shows the radial profiles of the measured peak velocity on the log r-log v plane. The red solid and dashed lines show the rotational and infalling velocity along the disk midplane in the models, i.e., $v_φ(r, θ = 90°)$ and $v_r(r, θ = 90°)$, respectively, described with Equations (6), (8), and (9). The third panel in the first row shows that the measured velocity for the freefall model is much higher than that from the observations, suggesting that the infalling or rotational velocity in the model is too large to explain the observations.

To examine a case where the rotational velocity is suppressed, we produced a disk-and-envelope model with slow rotation with $\rho_0 = 1.1 \times 10^{-17} \text{ g cm}^{-3}$ and $r_c = 150$ au, which correspond to a rotational velocity half that in the model with freefall. Note that the reference density was also modified because it is defined as a value at $r_c$. The second row in Figure 5 compares the model with slow rotation model and the observations. The first panel shows that the overall shape of the PV diagram of the model with slow rotation is very similar to that of the model with freefall but the former is more axisymmetric about the offset axis. Regardless of the lower rotational velocity, the intensity peaks appear at higher velocities than the observations, as seen in the second panel. The measured peak velocity was also much higher than the observed one, as shown in the third panel. These results suggest that the peak velocity is almost determined by the infalling velocity in the case where the infalling velocity is freefall.

The infalling velocities in the two models are too high to explain the observations. Although one might consider adopting a smaller $M_*$ to reduce the infalling velocity in the model, $M_*$ is well constrained from Keplerian rotation of the disk. Thus, we introduced a new parameter $\alpha$ as $v_r = \alpha v_{r,ff}$ to suppress the infalling velocity, maintaining the same stellar mass. Here, $v_{r,ff}$ is the radial velocity in the case of the freefall expressed by Equation (6). Strictly speaking, Ulrich's envelope model does not show self-consistency with the additional parameter $\alpha$; nonetheless, this approach would still be acceptable to observe a tendency of the velocity structures with suppressed infalling velocity. As a case where the infalling velocity was lower than the freefall velocity, a model with $\alpha = 0.4$ was produced and compared to the observations. In the second panel in the third row in Figure 5, the model with slow infall reproduces the overall shape of the observed PV

<footnote n="10"><sup>10</sup> http://www.ita.uni-heidelberg.de/~dullemond/software/radmc-3d/</footnote>

<footnote n="11"><sup>11</sup> https://www2.mpi-ahd.mpg.de/home/henning/Dust_opacities/Opacities/opacities.html</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** The first and second columns compare PV diagrams cut along the disk major axis between the observations (black contours and gray scale) and models (red contours). Contour levels of both observed and model PV diagrams are the same as those in Figure 3. Red data points are the peak velocities measured with the model PV diagrams. The third column shows radial profiles of the measured peak velocity. Black and red data points are peak velocities measured with the observed and model PV diagrams, respectively. Black solid lines represent the best-fit power-law function for the data points from the observations. The red solid and dashed lines show the azimuthal and radial velocity components along the disk midplane in the models, i.e., $v_\phi(r, \theta = 90^\circ)$ and $v_r(r, \theta = 90^\circ)$, respectively, described with Equations (6), (8), and (9). From the top row to the bottom, results of the model with freefall, slow rotation, slow infall, and modified $v_\phi$ are presented.</figcaption></figure>

diagram well. The radial profile of the measured peak velocity, indicated in the third panel, is also consistent with the profile measured from the observations at radii smaller than the break

radius of ~2900 au. Although the measured peak velocity is a blend of the rotational velocity and infalling velocity, its absolute values are very similar to the given rotational velocity.
Among the three disk-and-envelope models discussed above, the one with slow infall ($\alpha = 0.4$) reproduces the observations best, suggesting that the infall motion in the envelope is slower than the freefall motion. Nevertheless, we should stress that none of them can reproduce the radial profile of the observed peak velocities at radii larger than 2900 au. The nature of the radial profile of the measured peak velocity, including the increase in velocity outside the break radius, will be further discussed in the next section.

## 5. Discussion

### 5.1. Nature of the Radial Profile of the Measured Peak Velocity

The radial profile of the measured peak velocity shows a break at a radius of ~2900 au, as shown in Section 4.1. The profile inside the break radius is consistent with an $r^{-1}$ power law, and this is indeed expected if the motions are dominated by rotation in an infalling and rotating envelope with a constant $j$, as was demonstrated using the disk-and-envelope model with slow infall. One might wonder whether the observed velocity structures could originate from non-axisymmetric infalling flows suggested in a previous work (Yen et al. 2014) rather than rotational motion of a more spherical envelope. However, the extent of the suggested infalling flows is only ~10" on the plane of the sky, while the emission explained by the spherical envelope model extends up to a radius of ~20" from the protostellar position. The previous ALMA observations in Yen et al. (2014) would have detected only the density-enhanced non-axisymmetric infalling flows within a more spherical envelope due to missing flux. The limited field of view of the previous ALMA observations, whose FWHM of the primary beam is ~15", could also be a reason why the emission extending up to a radius of ~20" from the protostellar position was not detected in Yen et al. (2014). The infalling flows reported in the previous work are almost along the north-south direction, which also suggests that they are not likely responsible for the velocity structure in the PV diagram cut along the disk major axis (P.A. = 54°) in the current work.

The measured peak velocity outside the break radius, on the other hand, increases or is almost flat. A simple explanation of the increase of the measured peak velocity outside the break radius may be given by a change in the radial dependence of the rotational velocity at the break radius. Such a change is possible if the infalling material carries a greater specific angular momentum at outer radii in the envelope (Yen et al. 2011; Takahashi et al. 2016). In order to demonstrate this possibility, the disk-and-envelope model with slow infall was modified to have the following rotational velocity:

$$v_{\phi}(r, \theta) = \begin{cases} v_{\text{break}} \left( \frac{r}{r_{\text{break}}} \right)^{-1} \frac{\sin \theta_0}{\sin \theta} & (r \le r_{\text{break}}) \\ v_{\text{break}} \left( \frac{r}{r_{\text{break}}} \right)^{p_{\text{out}}} \frac{\sin \theta_0}{\sin \theta} & (r > r_{\text{break}}) \end{cases}, \quad (9)$$

where $v_{\text{break}} = \sqrt{GM_*r_c \sin \theta_0^2}/r_{\text{break}}$, $r_{\text{break}} = 2900$ au, and $p_{\text{out}} = 0.3^{12}$. The PV diagram of the model with the slow infall and modified $v_{\phi}$, presented in the fourth row in Figure 5, is very similar to that obtained with the original model with slow infall (before the modification), shown in the third row, while the

radial profile of the peak velocity obtained from the model with slow infall and modified $v_{\phi}$ reproduces the observed radial profile well even at radii larger than 2900 au. The increase in the measured peak velocity could be caused by a change in the radial dependence of the infalling velocity instead of the rotational velocity. Disk-and-envelope models with such changes, however, fail to explain the observations, as demonstrated in Appendix A.1.

Note that the break in the profile could also occur if the dominant mechanism creating the velocity gradient changes from rotation to another, such as infalling flow and turbulence. Although non-axisymmetric infalling flows were suggested in a previous work (Yen et al. 2014), these flows extend only ~10" (~1400 au) and are almost along the north-south direction as was mentioned above. Hence, we examined whether an infalling flow coming from another direction could reproduce the increase in the peak velocity outside the break radius in Appendix A.2. As shown in Figure A3, the infalling flow model barely reproduces the increase in the peak velocity outside 2900 au, although the direction of the infalling flow has to be very specific, as explained in Appendix A.2. Another possible motion to explain the increase in the peak velocity would be cloud-scale turbulence characterized by Larson's law ($\Delta v \propto r^{-0.5}$; Larson 1981; McKee and Ostriker 2007). Cloud-scale turbulence can produce a velocity gradient on subparsec scales in projected maps (Burkert & Bodenheimer 2000), and the velocity deviation is expected to be proportional to $r^{0.5}$. Thus, the increase in the peak velocity outside the break radius could be explained as the tail of the scaling law of the cloud-scale turbulence, as also discussed by Gaudel et al. (2020). On the other hand, the velocity gradient in L1489 IRS appears systematic and coherent. Table 1 shows that the directions of the velocity gradient have only a small dispersion of ~±5° over a spatial scale of 10"–60" (~1400–8400 au). Numerical simulations suggest that the angular momentum axis often varies by ~20°–90° with a radius from 1000 to 10,000 au in turbulent dense cores (Joos et al. 2013; Matsumoto et al. 2017). The velocity gradient in the projected map is also expected to change the directions with radius for such a turbulent velocity field. Hence, it would be more natural to interpret the coherent velocity structure of L1489 IRS as rotational motion rather than turbulence.

In summary, a simple and probably more natural explanation of the increase in the measured peak velocity outside the break radius is that the radial dependence of the rotational velocity changes at the break radius. An infalling flow or turbulence could also explain the increasing peak velocity if such motions were more dominant than rotation outside the break radius. In the following sections, we further discuss the simple and more natural case where the radial dependence of the rotational velocity changes at the break radius.

### 5.2. Specific Angular Momentum Profile in L1489 IRS

As the radial profile of the measured peak velocity is well explained by the rotational motion of the envelope, the specific angular momentum can be calculated from the measured peak velocity with the relation $j = r \times v$ after a correction of the inclination angle of 73°. Again, the measured peak velocity is not exactly, but almost, comparable to the rotational velocity, as shown in the comparison to the models. The difference between the measured peak velocity and the expected rotational

<footnote n="12">12 Equation (9) for $r \le r_{\text{break}}$ is another expression of Equation (8), and derived from Equation (10) in Ulrich (1976) and Equation (8).</footnote>
<figure><img src="image_8.png" /><figcaption>Figure 6. Radial profile of the specific angular momentum measured in L1489 IRS (data points) and those profiles at certain times in a calculation of the inside-out collapse model (solid and dashed lines).</figcaption></figure>

velocity, with a correction of the inclination angle in the slow-infall model, is about 10%. The specific angular momentum of the disk and envelope within 1000 au was also derived from the rotation curve measured in the previous work of Sai et al. (2020). The derived specific angular momentum is plotted as a function of the radius in Figure 6. The uncertainty of the derived specific angular momentum due to the uncertainty of the inclination angle between 60° and 90° is roughly the size of the symbol in the plot. Figure 6 shows that the specific angular momentum profile consists of three different regimes: a Keplerian disk inside ~600 au, a j-constant regime at ~600–2900 au, and a j-increasing regime outside ~2900 au.

A kinematic transition between the j-constant and j-increasing regimes has been reported in Class 0 sources (Gaudel et al. 2020). They calculated apparent specific angular momenta from LSR velocities measured at radii of 50–5000 au for 12 Class 0 protostellar systems. Although the obtained radial dependences of the apparent specific angular momentum were largely scattered from one source to another, the mean profile showed a break at a radius of ~1600 au. The mean apparent specific angular momentum in the j-constant regime is ~6 × 10⁻⁴ km s⁻¹ pc. These values are smaller than those in L1489 IRS, i.e., r<sub>break</sub> ~ 2900 au and j ~ 4 × 10⁻³ km s⁻¹ pc. This is consistent with simple analytical calculations of angular momentum transfer from an initial dense core where specific angular momentum is larger at a larger radius (Yen et al. 2011, 2017; Takahashi et al. 2016).

In order to further investigate the nature of the specific angular momentum profile and to discuss the size of the SMFZ, we performed simple analytical model calculations of the angular momentum transfer from an initial dense core to an infalling envelope and compared the calculated radial profiles of the specific angular momentum to that measured in the observations.

First, we calculate the angular momentum transfer along an equatorial plane of an initial dense core, based on the inside-out collapse of a singular isothermal sphere with a finite angular momentum (Shu 1977; Yen et al. 2011, 2017). We assume that

the specific angular momentum profile in the initial dense core is given by

$$j = j_0 \left( \frac{r}{10^4 \text{ au}} \right)^{1.6}, \qquad (10)$$

as suggested by previous observations of dense cores (Goodman et al. 1993; Caselli et al. 2002a; Tatematsu et al. 2016; Pineda et al. 2019; Gaudel et al. 2020). The density distribution of the initial dense core is assumed as follows (Shu 1977):

$$\rho = \frac{c_s^2}{2\pi G} r^{-2}, \qquad (11)$$

where c_s is the isothermal sound speed. The enclosed mass including the central object mass is calculated from the density distribution as

$$M(r) = \frac{2c_s^2}{G} r. \qquad (12)$$

Assuming that infalling material follows freefall, infalling motion of material is calculated from the equation of motion of the gas element:

$$\frac{dv_r}{dt} = -\frac{GM}{r^2}. \qquad (13)$$

Here M is the enclosed mass within the initial radius r<sub>ini</sub>, at which the infalling material is initially located. The left-hand side of the equation can be written as

$$\frac{dv_r}{dt} = \frac{dr}{dt} \frac{dv_r}{dr} = \frac{1}{2} \frac{dv_r^2}{dr}. \qquad (14)$$

By substituting this in Equation (13) and integrating, the following equations are obtained:

$$\begin{gathered} \frac{1}{2} v_r^2 = \int_{r_{\text{ini}}}^{r} -\frac{GM}{r'^2} dr', \\ \therefore v_r = \frac{dr}{dt} = -\left\{2GM\left(\frac{1}{r} - \frac{1}{r_{\text{ini}}}\right)\right\}^{1/2}, \end{gathered} \qquad (15)$$

and thus

$$dt = - \left\{ 2GM \left( \frac{1}{r} - \frac{1}{r_{\text{ini}}} \right) \right\}^{-1/2} dr. \qquad (16)$$

Note that M = M(r<sub>ini</sub>) is constant on the frame of the gas element and v<sub>r</sub> is negative because infalling motion is considered. By replacing r with the normalized radius x ≡ r/r<sub>ini</sub> and integrating Equation (16), the relation between the radial location and time is obtained:

$$t = - \sqrt{\frac{r_{\text{ini}}^3}{2GM}} \int_{1}^{x} \frac{dx'}{\sqrt{x'^{-1} - 1}}. \qquad (17)$$

In the inside-out collapse model, the expansion wave propagates from the inside to the outside at the isothermal sound speed c_s, and gas located at a certain radius starts to infall once the expansion wave has reached there. Therefore, we set t = 0 at the time when the innermost part of the core collapses, and a gas element at r = r<sub>ini</sub> starts infalling when the expansion wave reaches r<sub>ini</sub>, i.e., at t = t<sub>wave</sub> = r<sub>ini</sub>/c_s. The radial
distribution of the specific angular momentum at *t* is then derived assuming conservation of angular momentum, *j*(*r*(*t*)) = *j*(*r*<sub>ini</sub>). The protostellar mass at *t* is also estimated from *M*(*r*<sub>ini</sub>) with *r*<sub>ini</sub> satisfying *r*(*t*, *r*<sub>ini</sub>) < 0.

Figure 6 shows the result of a calculation with $c_s = 0.42 \text{ km s}^{-1}$ and the initial specific angular momentum profile with $j_0 = 1.8 \times 10^{-2} \text{ km s}^{-1} \text{ pc}$, which is roughly consistent with recent measurements of specific angular momentum in two Class 0 sources and a candidate for the first hydrostatic core (Pineda et al. 2019). The expected angular momentum profile at the age of $8 \times 10^4 \text{ yr}$, shown in Figure 6, matches the measured profile including that of the $j$-increasing regime at radii larger than 2900 au. The central stellar mass at the discussed age is $\sim 1.5 M_\odot$, which is comparable to the protostellar mass of L1489 IRS. Note that the assumed $c_s$ corresponds to a temperature of $\sim 50$ K, whereas a typical temperature of dense cores is $\sim 10$ K (Benson & Myers 1989). These results suggest that the initial core of L1489 IRS was about five times more massive or denser than one in hydrostatic equilibrium. Such a massive initial dense core is actually expected for L1489 IRS because an infalling velocity lower than the freefall velocity suggests that an additional force supports infalling material against gravity, as discussed in the next section. Note that computed specific angular momentum profiles basically do not change except for the timescale even if slow infall is considered.

Another analytical model proposed by Takahashi et al. (2016) was also compared to the case of L1489 IRS. The model calculates the collapse of a supercritical Bonnor–Ebert sphere assuming that the entire dense core collapses simultaneously without following the inside-out collapse. Gas pressure is considered in the model in addition to gravity:

$$ \frac{dv_r}{dt} = \frac{GM}{fr_{ini}r} - \frac{GM}{r^2}. \quad (18) $$

$f(r_{ini})$ is proportional to the ratio between the gravitational potential and the thermal energy density at $r_{ini}$, $-\frac{GM}{r_{ini}\beta c_s^2}$, where

$\beta(r_{ini}) \equiv \frac{\partial \ln \rho}{\partial \ln r} = \frac{\partial \ln \rho_{ini}}{\partial \ln r_{ini}}$ assuming that $\frac{\partial \ln \rho}{\partial \ln r}$ is constant in time. From consideration similar to that for the inside-out collapse model, the relation between the location and time is derived as follows:

$$ t = -\sqrt{\frac{r_{ini}^3}{2GM}} \int_1^x \frac{dx'}{\sqrt{f^{-1} \ln x' + x'^{-1} - 1}}. \quad (19) $$

In the calculation, we adopt a typical core radius $R_{core} = 10,000$ au, the central density $n_c = 2.3 \times 10^5 \text{ cm}^{-3}$, and $c_s = 0.2 \text{ km s}^{-1}$ for a Bonnor–Ebert sphere and set $f=3$ to have the Bonnor–Ebert sphere unstable and explain the protostellar mass of L1489 IRS.

Figure 7 compares the observations and the result of the calculation assuming $j_0 = 9 \times 10^{-3} \text{ km s}^{-1} \text{ pc}$. The angular momentum profile at $7.6 \times 10^4 \text{ yr}$ after the start of the collapse roughly matches the measured profile including that of the $j$-increasing regime at $\sim 2900-5000$ au. The central stellar mass was estimated to be $\sim 1.6 M_\odot$, which is consistent with that of L1489 IRS. The age of the protostar was calculated as $t_{sys} - t_{prop}$, where $t_{sys}$ is the time after the collapse starts and $t_{prop}$ is the epoch of protostar formation (see Takahashi et al. 2016 for more detail). In the calculations, the protostellar age was estimated to be $\sim 2.8 \times 10^4 \text{ yr}$.

<figure><img src="image_9.png" /><figcaption>Figure 7. Same as Figure 6 but with the result for the model based on Takahashi et al. (2016).</figcaption></figure>

Our model calculations yield a protostellar age of $\sim (3-8) \times 10^4 \text{ yr}$. The typical lifetime of Class I sources has been estimated to be $\sim (0.9-5) \times 10^5 \text{ yr}$ from the relative numbers of sources in different evolutionary stages (Wilking et al. 1989; Kenyon et al. 1990; Evans et al. 2009; Dunham et al. 2015; Kristensen & Dunham 2018). Thus, the minimum value of the typical lifetime of Class I sources is comparable to, or $\sim 3$ times larger than, the protostellar ages calculated in our models. Note that the lifetime derived from the previous observations mentioned above represents median lifetime or half-life of a sample of Class I protostars, and true lifetime has some diversity (Evans et al. 2009; Kristensen & Dunham 2018). In fact, the lifetime of individual dense cores depends on the core density and can vary by an order of magnitude (Kirk et al. 2005). The protostellar ages estimated in our calculations can also be longer by a factor of $\sim 2-3$ if the infall velocity is lower than the freefall velocity, as suggested in Section 4.2, although current calculations assume the infall timescale is determined by the freefall time.

It should be noted that the analytical calculations adopted here were simplified, and did not consider the possible effects of magnetic field. Nonideal MHD simulations of disk formation showed that the radial profile of the specific angular momentum in the collapsing dense core could have a shallow slope at radii of $\sim 100-6000$ au when ambipolar diffusion was efficient in the simulations (Zhao et al. 2018). This shallow profile from the simulations could be similar to the green dashed curve at 1000–5000 au in Figure 7. However, the simulations stopped when the stellar mass was less than 0.1 $M_\odot$ and the disk radius was $\sim 20$ au, which are more than ten times smaller than those of L1489 IRS. Thus, it is difficult to make a direct comparison. Further simulations covering the later evolutionary stages or various initial conditions to produce larger stellar masses and disk radii are required to understand the observational results including the effects of magnetic field.

## 5.3. Stellar Mass Feeding Zone

Comparisons between simple analytical models and the observations demonstrate that the observed profile of the
specific angular momentum can be explained by the gravita-
tional collapse of dense cores with greater specific angular
momenta than the value measured by Goodman et al. (1993).
Based on these models, we estimate the radius of the SMFZ.
The size of the SMFZ can be estimated by considering from
which radius of the initial core the infalling material showing
the *j*-constant regime comes, because material having a specific
angular momentum smaller than that of the *j*-constant regime
would have already accreted into the disk or central protostar.
In the case of the model based on the inside-out collapse,
shown in Figure 6, the radius of the SMFZ was estimated to be
~4000 au by extrapolating the *j*-constant profile to the initial
distribution of the specific angular momentum shown by the
black line. Similarly, in the case of the collapsing supercritical
Bonnor–Ebert sphere, shown in Figure 7, the radius of the
SMFZ was estimated to be ~6000 au.

The estimated radius of the SMFZ may not be the final value for L1489 IRS because it infalls from the envelope still continues. Regardless, we expect little material to accrete to the disk and the protostar because L1489 IRS is close to the end of the Class I stage. The disk radius of L1489 IRS, ~600 au, which is comparable to the maximum disk radius of single Class II sources in the Taurus star-forming region (~640 au; Guilloteau et al. 2014; Simon et al. 2017), also suggests that little material accretes to the disk and the protostar in L1489 IRS. Even if the disk radius of L1489 IRS increases from ~600 au to 640 au, the size of the SMFZ increases by only 10%, which is less significant.

It would be interesting to estimate sizes of SMFZs for Class II objects in the post-accretion phase because these sizes should be the final ones. However, it is difficult to estimate them for lack of information constraining specific angular momentum distributions in initial dense cores for Class II objects without infalling envelopes. We should stress that Class I objects embedded in infalling envelopes, such as L1489 IRS, allow us to measure the specific angular momentum in the outermost envelopes, thus constraining specific angular momentum distributions in initial dense cores, as shown in Figures 6 and 7.

It is interesting to note that the estimated radii of the SMFZ
are smaller than those of typical dense cores. In the model based
on a supercritical Bonnor–Ebert sphere, R<sub>core</sub> = 10, 000 au was
assumed and the estimated radius of the SMFZ was 60% of
R<sub>core</sub>. The total mass enclosed within 10,000 au was ~3.4 M<sub>☉</sub>,
while the protostellar mass was ~1.6 M<sub>☉</sub>, suggesting that ~50%
of the total mass in the initial dense core forms a star. This value
is close to the star formation efficiency suggested from
comparisons between the CMF and IMF, ~40% (André et al.
2014; Könyves et al. 2015), and our result based on the model
calculations suggests that only a limited area of the dense core
forms the protostar(s) within the timescale of the Class 0 and I
stages. The physical mechanism limiting the area of dense cores
to form stars has to be studied further in the future.

5.4. Regulated Infalling Velocity

Our modeling, presented in Section 4.2, suggests that the
infalling velocity in the envelope of L1489 IRS is lower than
the freefall velocity yielded by the stellar mass of 1.6 M☉ by a
factor of ~2.5. This implies that some force supports infalling
material against the gravitational force of the central protostar.

Magnetic field possibly supplies an additional force against
gravity. Although neither the strength nor morphology of
magnetic field in L1489 IRS is known, we can roughly estimate

the possible strength of the magnetic field based on assump-
tions in the same way as discussed in Aso et al. (2015). We
consider infalling motion along the equatorial plane of the
envelope in cylindrical coordinates. The cylindrical radius R is
taken to be the radius in the equatorial plane. It is assumed that
the magnetic field is perpendicular to the equatorial plane, and
pinched toward the center. Such pinched magnetic fields
showing an hourglass shape have been found in infalling
envelopes mostly around Class 0 protostars (e.g., Girart et al.
2006). Although L1489 IRS is a more evolved source, we
adopted it here since very few measurements of magnetic fields
exist at disk to envelope scales for Class I protostars. The
equation of motion can be described as follows under a
symmetry condition, where ∂/∂θ = 0 and ∂/∂z = 0 (Aso et al.
2015):

$$
\rho(R) v_r(R) \frac{dv_r}{dR} = - \frac{GM_* \rho(R)}{R^2} - \frac{1}{2\mu_0} \frac{dB^2}{dR} + \frac{B(R)^2}{\mu_0 R_{\text{curv}}(R)}, \quad (20)
$$

where $R_{\text{curv}}$ is the curvature radius of the magnetic field. We rewrite the equation in a more convenient form assuming the radial dependence of $B$ as $B_0(R/R_0)^{-q}$ (see Appendix C for more detail):

$$
B(R) = \left[ \frac{GM_* \rho(R)}{R^2} (1 - \alpha^2) \mu_0 \left\{ \frac{q}{R} + \frac{1}{R_{\text{curv}}(R)} \right\}^{-1/2} \right]^{1/2} . \quad (21)
$$

For an estimate of order, R<sub>curv</sub> = R was assumed because the
magnetic field morphology is unknown. Then, from the best
model parameters of α = 0.4 and ρ<sub>0</sub> = 1.4 × 10<sup>-18</sup> g cm<sup>-3</sup>, B
was estimated to be ~0.25 mG at a radius of 1000 au from
Equation (21), assuming B ~ ρ<sup>2/3</sup>, which is the case for
spherically symmetric collapse.

Zeeman observations on cloud and core scales have reported
that the strength of the magnetic fields ranges from ~1 µG to
1 mG (Crutcher et al. 2010). Observations in polarized
continuum emission, which traces magnetic field directions,
also estimated the magnetic field strength on 100–1000 au
scales to be of the order of magnitude of milligauss through the
method of Chandrasekhar & Fermi (1953) or the alternative
method proposed by Koch et al. (2012) (Girart et al. 2006; Rao
et al. 2009; Stephens et al. 2013). Recent Zeeman observations
of the Class II source TW Hya using ALMA in CN lines have
reported nondetection of Zeeman splitting in the CN line from
the protoplanetary disk around TW Hya, providing an upper
limit of |B<sub>z</sub>| < 0.8 mG (Vlemmings et al. 2019). The magnetic
field strength estimated for L1489 IRS is within the range of
observationally suggested values. Hence, the slow infall in
L1489 IRS could be explained by the magnetic field if the
magnetic field strength in L1489 IRS is similar to that in other
sources.

Similarly, a slow infalling velocity with α = 0.3 has been
reported in the Class I protostar TMC-1A (Aso et al. 2015).
Aso et al. (2015) estimated the strength of the magnetic field in
TMC-A to be ~2 mG at 200 au, which corresponds to 0.4 mG
at 1000 au if extrapolated using the relation B ~ ρ²/³. Although
L1489 IRS has a larger stellar mass than TMC-1A (0.68 M☉),
the values of α and the estimated strengths of the magnetic
fields in both sources are similar. It is also reported that the
infalling velocities in the Class 0 protostar L1527 IRS and the
Class I protostar L1551 IRS 5 are lower than the expected freefall velocity by factors of ~2 and ~3, respectively (Ohashi et al. 2014; Chou et al. 2014). Interestingly, all of these four sources suggest similar α values of 0.3–0.5 regardless of their different evolutionary stages. Lower infalling velocity could be common in the protostellar phase.

## 6. Summary

We have conducted mapping observations covering a ~2' × 2' region around the protostar L1489 IRS using ACA and IRAM 30 m in C¹⁸O 2–1 emission in order to kinematically investigate a zone feeding mass into the central star, the “stellar mass feeding zone”. The main results and conclusions are summarized as follows.

1. We detected intensity peaks of both the 1.3 mm continuum and C¹⁸O emission at the protostellar position, which traced an envelope associated with the protostar. The C¹⁸O emission shows a velocity gradient at $r \gtrsim 1000$ au in almost the same direction as that of the velocity gradient due to its disk rotation. The 1.3 mm continuum emission exhibits an intensity peak to the east of the protostar, which is associated with a starless core. The C¹⁸O emission shows the second intensity peak to the northeast of the protostar, suggesting another associated component.

2. The peak velocity was measured as a function of radius at radii of ~1000–5000 au to the southwest of the protostar with the PV diagram cut along the major axis of the disk rotation to characterize the observed velocity gradient. The measured peak velocity decreases with radius inside 2900 au, suggesting differential rotation of the envelope, but is constant or increases slightly outside 2900 au.

3. The measured peak velocity is considered to trace the rotational velocity of the envelope based on comparison of the disk-and-envelope models with the observations. The model with slow infall, wherein the infalling velocity is lower than the freefall velocity by a factor of 2.5 and the specific angular momentum is constant, best reproduced the observed PV diagram and the radial profile of the measured peak velocity within 2900 au. The model also suggested that the measured peak velocity was almost comparable to the expected rotational velocity of the envelope with a constant specific angular momentum. In order to explain the increase in the peak velocity outside 2900 au in radius, the specific angular momentum has to increase outside that radius in the model.

4. We calculated the radial profile of the specific angular momentum at ~300–5000 au from the peak velocity measured in current and previous works. The specific angular momentum profile consists of three parts: a Keplerian disk inside ~600 au, a *j*-constant regime at ~600–2900 au, and a *j*-increasing regime outside ~2900 au. We compared it to analytic models of collapsing dense cores having a finite angular momentum, assuming that infalling materials conserve their specific angular momenta. The measured specific angular momentum profile is explained by the models. The analytic models suggest that an initial core more massive or denser than a singular isothermal sphere or a critical Bonnor–Ebert sphere is preferred to explain the protostellar mass of L1489 IRS of 1.6 $M_\odot$.

5. Based on comparison of the radial profiles of the specific angular momentum derived from the observations and the analytical models, we estimated the size of the SMFZ, a zone where material forming the protostar resides in the initial dense core, to be ~4000–6000 au in radius in L1489 IRS. Such radii are significantly smaller than typical radii of dense cores of ~10,000–20,000 au, suggesting that only a part of the initial core feeds the protostar.

6. The infalling velocity being lower than the freefall velocity by a factor of ~2.5 can be explained by a magnetic field with a strength of $B \sim 0.25$ mG, assuming the morphology of the magnetic field. The estimated strength of magnetic field is consistent with that measured in molecular clouds and protostellar envelopes in previous observations.

This paper 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. This work is based on observations carried out under project number 136-18 with the IRAM 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). We would like to thank the IRAM staff for their support during the campaigns, and the NIKA2 collaboration for enabling tools to reduce the data. We also thank Stéphane Guilloteau for fruitful discussions on the merging of the ACA and IRAM 30 m data. We are grateful to Kengo Tomida, Doug Johnstone, Yuri Aikawa, Patrick M. Kock, Ya-Wen Tang, and Masayuki Yamaguchi for productive discussions. J.S. thanks Shigehisa Takakuwa for the financial support for his trip to the IRAM 30 m telescope to make observations. N.O. acknowledges a grant from the Ministry of Science and Technology (MOST) of Taiwan (MOST 109-2112-M-001-051 and MOST 110-2112-M-001-031). A.M. is supported by the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (Magnetic YSOs project, grant agreement No. 679937). J.S. is supported by Academia Sinica Institute of Astronomy and Astrophysics.

**Facilities:** ALMA, IRAM:30m

**Software:** CASA (McMullin et al. 2007), GILDAS (http://www.iram.fr/IRAMFR/GILDAS), RADMC-3D (Dullemond et al. 2012), Numpy (Oliphant 2006; van der Walt et al. 2011), Scipy (Virtanen et al. 2020), Astropy (Astropy Collaboration et al. 2013, 2018), Matplotlib (Hunter 2007).

## Appendix A
### Supplemental Models

Two supplemental models are provided to explore further possibilities explaining the increase in the peak velocity outside the radius of 2900 au.

#### A.1. Disk-and-envelope Model with Modified $v_r$

Disk-and-envelope models with modified $v_r$, where infalling velocity increases outside the break radius, were compared with the observations to determine whether such models could
<figure><img src="image_10.png" /><figcaption>Figure A1. Same as Figure 5 except for the two supplemental models, which have modified radial velocity distributions.</figcaption></figure>

explain the increase in peak velocity outside 2900 au. Two models following Equations (5) to (8) except for the radial velocity distribution were constructed with the parameters in Table 2 and $\alpha = 0.4$. One model, referred to as $v_r$-model 1, has the following radial velocity distribution:

$$v_r(r, \theta) = \begin{cases} -\alpha \left(\frac{GM_*}{r}\right)^{0.5} \left(1 + \frac{\cos \theta}{\cos \theta_0}\right)^{0.5} & (r \le r_{\text{break}}) \\ -v_{r, \text{break}} \left(\frac{r}{r_{\text{break}}}\right)^{p_{\text{inf}}} & (r > r_{\text{break}}) \end{cases} \quad (A1)$$

where $v_{r, \text{break}} = 0.4 \text{ km s}^{-1}$ and $p_{\text{inf}}=1.2$ so that the radial velocity increases with radius outside the break radius and reaches the freefall velocity at the edge of the envelope. The other model, $v_r$-model 2, has the radial velocity distribution

$$v_r(r, \theta) = \begin{cases} -\alpha \left(\frac{GM_*}{r}\right)^{0.5} \left(1 + \frac{\cos \theta}{\cos \theta_0}\right)^{0.5} & (r \le r_{\text{break}}) \\ -\left(\frac{GM_*}{r}\right)^{0.5} \left(1 + \frac{\cos \theta}{\cos \theta_0}\right)^{0.5} & (r > r_{\text{break}}) \end{cases} \quad (A2)$$

so that the radial velocity is lower than the freefall velocity inside the break radius but is the freefall velocity outside.

Figure A1 shows the model results. The solid and dashed lines in the third column in Figure A1 show the given radial velocity for the models. Although both models have radial velocities increasing outside the break radius, the observed peak velocity cannot be reproduced.

## A.2. Disk-and-flow Model

A disk-and-flow model was built to investigate whether a flow-like infall such as those reported by Yen et al. (2014) and Pineda et al. (2020) can explain the observed velocity structure outside the break radius or not. We modified the envelope structure of the disk-and-envelope model described in Section 4.2 to build a disk-and-flow model that has a flow falling from certain directions. The direction of the flow is determined by the parameters $\theta_0$ and $\phi_0$, which are the initial location of the infalling material in $\theta$ and $\phi$ at an infinite radius, respectively. We adopted $\theta_0 = 95^{\circ} \pm 5^{\circ}$, where $\theta = 90^{\circ}$ is the disk midplane, so that the direction of the flow is along the direction of the disk major axis (i.e., the direction of the PV cut). A case with $\phi_0$ of $326^{\circ} \pm 2^{\circ}$ was compared with the observations. Other parameters describing the envelope structures are the same as those of the base disk-and-envelope model summarized in Table 2, except for the volume density $\rho_0$, which was fixed at $4 \times 10^{-17} \text{ g cm}^{-3}$. The emission from the flow component of the disk-and-flow model is presented in Figure A2, and shows the trajectory and velocity of the flow projected on the plane of the sky.

The PV diagram of the disk-and-flow model cut along the disk major axis shown in Figure A3 demonstrates that the model does not reproduce the observed PV diagram enlarged in the middle panel of Figure A3. The radial profile of the peak velocity measured in the disk-and-flow model, presented in the third panel of Figure A3, shows that the measured peak velocity slightly increases with radius at $\sim 1400-4000$ au regardless of the decrease in $v_r$ and $v_\phi$. This is because the infalling flow is more parallel to the plane of the sky at inner radii, resulting in a smaller line-of-sight velocity at smaller radii. Although the measured peak velocity of the model cannot explain the observed profile of the peak velocity inside the
<figure><img src="image_11.png" /><figcaption>**Figure A2.** Velocity channel maps of the disk-and-flow model (red contours) overlaid on the channel maps of the observations (background color and black contours). Contour levels and symbols are the same as those in Figure 2. Dashed lines show the direction of the PV cut (P.A. = 54°).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure A3.** Same as Figure 5 except for the disk-and-flow model.</figcaption></figure>

break radius, it is roughly consistent with the observed profile outside the break radius.

It should be noted that the direction of the flow has to be very specific to explain the observed velocity structure outside the break radius. The adopted initial angles of ($\theta_0$, $\phi_0$) = (95° ± 5°, 326° ± 2°) were analytically estimated so that the flow reproduces the observed velocity structure outside the break radius best. Although we have investigated different sets of the initial angles and flow widths, the increase in the peak velocities outside the break radius was not reproduced better than the current model.

## Appendix B
$N_2H^+ J = 1-0$

We present the data of $N_2H^+$ 1-0 and a brief discussion about the spatial distributions of the $N_2H^+$ 1-0 and $C^{18}O$ 2-1 emissions. Figure B1 shows the moment 0 map of the $N_2H^+$ emission with the 1.3 mm continuum map. The $N_2H^+$ emission was integrated over the velocity range from -1.8 to 13.9 km s⁻¹ to include all hyperfine components of the emission. The $N_2H^+$ emission shows a peak at the position of the 1.3 mm continuum second peak, and is likely associated with the starless core to the east of the protostar. The overall distribution of the $N_2H^+$ emission is consistent with previous observations (Caselli et al. 2002a). The spatial distributions of the $C^{18}O$ and $N_2H^+$ emissions are apparently anticorrelated (the primary and secondary peaks of $C^{18}O$ 2-1 are denoted by a green triangle and a yellow square in Figure B1).

To investigate the anticorrelation of these two molecular emissions quantitatively, the abundances of $C^{18}O$ and $N_2H^+$ molecules were estimated from their column densities and the

<figure><img src="image_13.png" /><figcaption>**Figure B1.** The moment 0 map of the $N_2H^+$ 1-0 (background color) and the 1.3 mm continuum map (contours). Contour levels start at $3\sigma$ and increases in steps of $2\sigma$, where $1\sigma$ corresponds to 5 mJy beam⁻¹. The dashed circle and the filled ellipse in the bottom-left corner denote the beam sizes for the 1.3 mm continuum map and the $N_2H^+$ map, respectively. The central cross denotes the protostellar position. Colored markers show positions from where column densities are derived.</figcaption></figure>

$H_2$ column densities. The measurements were made at three different points shown in Figure B1: the $C^{18}O$ primary and secondary peaks, and the $N_2H^+$ peak. We used $C^{18}O$ and continuum images smoothed with the beam of the $N_2H^+$ map for calculations.
<table><caption>**Table B1**
Results of the Fitting to Hyperfine Components of the N₂H⁺ 1–0 Emissions</caption><thead><tr><th>Position</th><th>R.A.<br>(J2000)</th><th>Decl.<br>(J2000)</th><th>T<sub>ex</sub><br>(K)</th><th>τ<sub>tot</sub></th><th>ν<sub>cent</sub><br>(km s<sup>-1</sup>)</th><th>σ<sub>v</sub><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>C<sup>18</sup>O primary peak</td><td>04:04:43.07</td><td>+26:18:55.28</td><td>3.1 ± 0.1</td><td>4.6 ± 1.8</td><td>6.96 ± 0.02</td><td>0.23 ± 0.01</td></tr><tr><td>C<sup>18</sup>O secondary peak</td><td>04:04:44.17</td><td>+26:19:39.59</td><td>4.1 ± 0.2</td><td>3.7 ± 0.6</td><td>6.758 ± 0.003</td><td>0.148 ± 0.003</td></tr><tr><td>N<sub>2</sub>H<sup>+</sup> peak</td><td>04:04:47.40</td><td>+26:19:17.85</td><td>6.58 ± 0.04</td><td>8.9 ± 0.2</td><td>6.7779 ± 0.0005</td><td>0.1115 ± 0.0005</td></tr></tbody></table>

<table><caption>**Table B2**
Measured Column Densities and Molecular Abundances</caption><thead><tr><th>Position</th><th>N<sub>H2</sub><br>(×10<sup>21</sup> cm<sup>-2</sup>)</th><th>N<sub>C<sup>18</sup>O</sub><br>(×10<sup>15</sup> cm<sup>-2</sup>)</th><th>N<sub>N<sub>2</sub>H<sup>+</sup></sub><br>(×10<sup>12</sup> cm<sup>-2</sup>)</th><th>X<sub>C<sup>18</sup>O</sub><br>(×10<sup>-7</sup>)</th><th>X<sub>N<sub>2</sub>H<sup>+</sup></sub><br>(×10<sup>-10</sup>)</th><th>X<sub>C<sup>18</sup>O</sub>/X<sub>N<sub>2</sub>H<sup>+</sup></sub><br>(×10<sup>2</sup>)</th></tr></thead><tbody><tr><td>C<sup>18</sup>O primary peak</td><td>8.7 ± 0.2</td><td>1.211 ± 0.009</td><td>4.3 ± 1.7</td><td>1.39 ± 0.03</td><td>4.9 ± 2.0</td><td>2.8 ± 1.1</td></tr><tr><td>C<sup>18</sup>O secondary peak</td><td>12.4 ± 0.6</td><td>1.53 ± 0.01</td><td>3.2 ± 0.5</td><td>1.23 ± 0.06</td><td>2.6 ± 0.4</td><td>4.8 ± 0.8</td></tr><tr><td>N<sub>2</sub>H<sup>+</sup> peak</td><td>24.1 ± 0.6</td><td>0.92 ± 0.01</td><td>11.6 ± 0.3</td><td>0.38 ± 0.01</td><td>4.8 ± 0.2</td><td>0.79 ± 0.04</td></tr></tbody></table>

Dust column densities are derived from the following equation assuming an optically thin emission:

$$ \Sigma_{\text{dust}} = \frac{I_{\nu}}{B_{\nu}(T_{\text{dust}}) \kappa_{\nu}}, \qquad (\text{B1}) $$

where $I_ν$ is the intensity, $B_ν$ is the Planck function, $T$ is the dust temperature, and $\kappa_ν$ is the dust-mass opacity. For the dust temperature, it was assumed that $T_{\text{dust}} = 20$ K at the C¹⁸O primary peak position and $T_{\text{dust}} = 10$ K at the other two positions. We applied $\kappa_ν = 0.90$ g cm⁻² derived by Ossenkopf & Henning (1994), the gas-to-dust mass ratio of 100, and the mean molecular weight $\mu = 2.8$ (Kauffmann et al. 2008) to derive the dust column density and convert it to the column density of H₂ gas.

The method described in Caselli et al. (2002b) was adopted to calculate the column densities of C¹⁸O and N₂H⁺ molecules. For the C¹⁸O emission, $T_{\text{ex}} = 20$ K at the C¹⁸O primary peak position and $T_{\text{ex}} = 10$ K at the other two positions were assumed. For the N₂H⁺ emission, we performed fitting to its hyperfine components with the Python package PySpecKit and derived the N₂H⁺ column density from the total opacity $\tau_{\text{tot}}$, the intrinsic velocity dispersion $\sigma_v$, and the excitation temperature $T_{\text{ex}}$. Table B1 lists the results of the fitting to the hyperfine components of the N₂H⁺ emission at three positions. Molecular abundances of C¹⁸O and N₂H⁺ were derived from the ratio of the column density of each molecule to that of the H₂ gas. The derived column densities and molecular abundances are summarized in Table B2.

The derived N₂H⁺ abundances are similar within a factor of 2 among the three measurement points, while those of C¹⁸O are different between the C¹⁸O primary/secondary peaks and the N₂H⁺ peak by a factor of ~3–4. Because of the difference in C¹⁸O abundances, the abundance ratios $X_{C^{18}O}/X_{N_2H^+}$ at the C¹⁸O primary/secondary peaks are ~4–6 times larger than that at the N₂H⁺ peak.

C¹⁸O abundances at the C¹⁸O primary/secondary peaks are comparable to the molecular cloud values in the Taurus region (Frerking et al. 1982), and to the values measured in other Class I objects (Jørgensen et al. 2002). A possible reason for the lower C¹⁸O abundance at the N₂H⁺ peak would be that CO molecules including C¹⁸O freeze-out onto dust grains, which is often seen in starless cores (e.g., Bergin et al. 2002; Tafalla et al. 2004). The N₂H⁺ emission is likely associated with the

L1489 starless core (Motte & André 2001; Wu et al. 2019), supporting our idea of a CO freeze-out at the N₂H⁺ peak. Aikawa et al. (2015) shows a relation between CO abundance and the H₂ gas temperature based on calculation of a chemical network. Their calculation suggests that CO abundance can be reduced by an order of magnitude due to a CO freeze-out. Thus, the lower abundance near the center of the starless core could be reproduced by CO freeze-out. The molecular abundance of C¹⁸O at the C¹⁸O second peak is similar to that measured at the protostellar position. However, the C¹⁸O second peak is at the edge of the starless core (Motte & André 2001; Wu et al. 2019) and the volume density is expected to be lower than at the center of the starless core, resulting in a longer timescale for CO freeze-out. Photodesorption would also be more effective at lower density, suppressing the CO depletion at the C¹⁸O second peak. Hence, the anticorrelation between the N₂H⁺ and C¹⁸O emissions would be the CO freeze-out at the starless core.

## Appendix C
### Derivation of Equation (21)

Here, we present the derivation of Equation (21). The radial component of the equation of motion is as follows (Aso et al. 2015):

$$ \begin{aligned} \rho(R)v_r(R)\frac{dv_r}{dR} &= -\frac{GM_*\rho(R)}{R^2} - \frac{1}{2\mu_0}\frac{dB^2}{dR} + \frac{B(R)^2}{\mu_0 R_{\text{curv}}(R)}, \\ &= F_G + F_B, \end{aligned} \qquad (\text{C1}) $$

where

$$ F_G = -\frac{GM_*\rho(R)}{R^2}, \qquad (\text{C2}) $$

$$ \begin{aligned} F_B &= F_{\text{mp}} + F_{\text{int}} \\ &= -\frac{1}{2\mu_0}\frac{dB^2}{dR} + \frac{B(R)^2}{\mu_0 R_{\text{curv}}(R)}. \end{aligned} \qquad (\text{C3}) $$

Here, $F_{\text{mp}}$ is the magnetic pressure and $F_{\text{int}}$ is the magnetic tension. $F_G$ is always negative from Equation (C2).

Now, we assume $v_r(R)$ to be $v_r(R) = \alpha v_{\text{ff}}(R)$, where $v_{\text{ff}}(r) = -\sqrt{2GM_*/r}$ is the freefall velocity and $\alpha$ is constant. Note that this assumption means that the potential is negative at
any radius and zero at an infinite radius, i.e., $|F_G| > |F_B|$. The derivative of $v_r$ is

$$ \begin{aligned} \frac{dv_r(R)}{dR} &= \alpha \frac{dv_{ff}}{dR} \\ &= \alpha \frac{1}{2} \sqrt{\frac{2GM_*}{R^3}}. \end{aligned} \quad (C4) $$

Therefore, the left-hand side of Equation (C1) can be written as

$$ \begin{aligned} \rho(R)v_r(R) \frac{dv_r}{dR} &= -\rho(R)\alpha^2 \frac{GM_*}{R^2} \\ &= \alpha^2 F_G. \end{aligned} \quad (C5) $$

Thus, from Equations (C1) and (C5),

$$ \alpha^2 = 1 + \frac{F_B}{F_G}, \quad (C6) $$

$$ \therefore \alpha = \sqrt{1 + \frac{F_B}{F_G}}. \quad (C7) $$

Assuming $B = B_0(R/R_0)^{-q}$, the term $F_{mp}$ can be written as follows:

$$ F_{mp} = -\frac{1}{2\mu_0} \frac{dB^2}{dR} \quad (C8) $$

$$ \begin{aligned} &= \frac{q}{\mu_0} \frac{B_0^2}{R_0} \left(\frac{R}{R_0}\right)^{-2q-1} \\ &= \frac{1}{\mu_0} B_0^2 \left(\frac{R}{R_0}\right)^{-2q} \frac{q}{R} \\ &= \frac{1}{\mu_0} B(R)^2 \frac{q}{R}. \end{aligned} \quad (C9) $$

Therefore, from Equation (C3), we obtain

$$ F_B = \frac{B(R)^2}{\mu_0} \left\{ \frac{q}{R} + \frac{1}{R_{\text{curv}}(R)} \right\}. \quad (C10) $$

From this equation and Equation (C6), we obtain the magnetic field B as a function of R and α:

$$ \begin{aligned} F_B &= \frac{B(R)^2}{\mu_0} \left\{ \frac{q}{R} + \frac{1}{R_{\text{curv}}(R)} \right\} = F_G(\alpha^2 - 1), \\ \therefore B(R) &= \left[ F_G(\alpha^2 - 1) \mu_0 \left\{ \frac{q}{R} + \frac{1}{R_{\text{curv}}(R)} \right\}^{-1} \right]^{1/2}. \end{aligned} \quad (C11) $$

Note that $F_B/F_G$ has no R-dependence according to Equation (C6) because α is assumed to be constant here.

## ORCID iDs

Jinshi Sai, (Insa Choi) https://orcid.org/0000-0003-4361-5577

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

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

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

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

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

Mathilde Gaudel https://orcid.org/0000-0002-5482-5206

## References

Adam, R., Adane, A., Ade, P. A. R., et al. 2018, A&A, 609, A115

Aikawa, Y., Furuya, K., Nomura, H., & Qi, C. 2015, ApJ, 807, 120

André, P., Di Francesco, J., Ward-Thompson, D., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 27

André, P., Men'shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102

Aso, Y., Ohashi, N., Aikawa, Y., et al. 2017, ApJ, 849, 56

Aso, Y., Ohashi, N., Saigo, K., et al. 2015, ApJ, 812, 27

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

Belloche, A. 2013, in EAS Publications Series, Vol. 62, ed. P. Hennebelle & C. Charbonnel (Paris: EDP Sciences), 25

Benson, P. J., & Myers, P. C. 1989, ApJS, 71, 89

Bergin, E. A., Alves, J., Huard, T., & Lada, C. J. 2002, ApJL, 570, L101

Brinch, C., Crapsi, A., Jørgensen, J. K., Hogerheijde, M. R., & Hill, T. 2007, A&A, 475, 915

Burkert, A., & Bodenheimer, P. 2000, ApJ, 543, 822

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

Caselli, P., Myers, P. C., & Thaddeus, P. 1995, ApJL, 455, L77

Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002b, ApJ, 565, 344

Chandrasekhar, S., & Fermi, E. 1953, ApJ, 118, 113

Chou, T.-L., Takakuwa, S., Yen, H.-W., Ohashi, N., & Ho, P. T. P. 2014, ApJ, 796, 70

Crutcher, R. M., Wandelt, B., Heiles, C., Falgarone, E., & Troland, T. H. 2010, ApJ, 725, 466

di Francesco, J., Evans, N. J. I., Caselli, P., et al. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil (Tucson, AZ: Univ. of Arizona Press), 17

Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A Multi-purpose Radiative Transfer Tool, Astrophysics Source Code Library, ascl:1202.015

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

Eisner, J. A. 2012, ApJ, 755, 23

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

Frerking, M. A., Langer, W. D., & Wilson, R. W. 1982, ApJ, 262, 590

Gaudel, M., Maury, A. J., Belloche, A., et al. 2020, A&A, 637, A92

Girart, J. M., Rao, R., & Marrone, D. P. 2006, Sci, 313, 812

Goodman, A. A., Benson, P. J., Fuller, G. A., & Myers, P. C. 1993, ApJ, 406, 528

Green, J. D., Evans, N. J., II, Jørgensen, J. K., et al. 2013, ApJ, 770, 123

Guilloteau, S., Simon, M., Piétu, V., et al. 2014, A&A, 567, A117

Harsono, D., Jørgensen, J. K., van Dishoeck, E. F., et al. 2014, A&A, 562, A77

Hogerheijde, M. R. 2001, ApJ, 553, 618

Hogerheijde, M. R., & Sandell, G. 2000, ApJ, 534, 880

Hogerheijde, M. R., van Dishoeck, E. F., Blake, G. A., & van Langevelde, H. J. 1998, ApJ, 502, 315

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

Joos, M., Hennebelle, P., Ciardi, A., & Fromang, S. 2013, A&A, 554, A17

Jørgensen, J. K., Schöier, F. L., & van Dishoeck, E. F. 2002, A&A, 389, 908

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, J. M., Ward-Thompson, D., & André, P. 2005, MNRAS, 360, 1506

Koch, P. M., Tang, Y.-W., & Ho, P. T. P. 2012, ApJ, 747, 79

Könyves, V., André, P., Men'shchikov, A., et al. 2015, A&A, 584, A91

Kristensen, L. E., & Dunham, M. M. 2018, A&A, 618, A158

Larson, R. B. 1981, MNRAS, 194, 809

Lee, C.-F. 2010, ApJ, 725, 712

Lee, C.-F., Hirano, N., Zhang, Q., et al. 2014, ApJ, 786, 114

Maret, S., Maury, A. J., Belloche, A., et al. 2020, A&A, 635, A15

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

McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565

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

Motte, F., & André, P. 2001, A&A, 365, 440

Murillo, N. M., Lai, S.-P., Bruderer, S., Harsono, D., & van Dishoeck, E. F. 2013, A&A, 560, A103

Ohashi, N., Hayashi, M., Ho, P. T. P., et al. 1997, ApJ, 488, 317

Ohashi, N., Saigo, K., Aso, Y., et al. 2014, ApJ, 796, 131

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

Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943
Perotto, L., Ponthieu, N., Macías-Pérez, J. F., et al. 2020, A&A, 637, A71

Pineda, J. E., Segura-Cox, D., Caselli, P., et al. 2020, NatAs, 4, 1158

Pineda, J. E., Zhao, B., Schmiedeke, A., et al. 2019, ApJ, 882, 103

Pokhrel, R., Myers, P. C., Dunham, M. M., et al. 2018, ApJ, 853, 5

Rao, R., Girart, J. M., Marrone, D. P., Lai, S.-P., & Schnee, S. 2009, ApJ, 707, 921

Sai, J., Ohashi, N., Saigo, K., et al. 2020, ApJ, 893, 51

Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611

Sheehan, P. D., & Eisner, J. A. 2017, ApJ, 851, 45

Shu, F. H. 1977, ApJ, 214, 488

Simon, M., Guilloteau, S., Di Folco, E., et al. 2017, ApJ, 844, 158

Stark, D. P., Whitney, B. A., Stassun, K., & Wood, K. 2006, ApJ, 649, 900

Stephens, I. W., Looney, L. W., Kwon, W., et al. 2013, ApJL, 769, L15

Tafalla, M., Santiago, J., Johnstone, D., & Bachiller, R. 2004, A&A, 423, L21

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

Tamura, M., Gatley, I., Waller, W., & Werner, M. W. 1991, ApJL, 374, L25

Tatematsu, K., Ohashi, S., Sanhueza, P., et al. 2016, PASJ, 68, 24

Ulrich, R. K. 1976, ApJ, 210, 377

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

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

Vlemmings, W. H. T., Lankhaar, B., Cazzoletti, P., et al. 2019, A&A, 624, L7

Wilking, B. A., Lada, C. J., & Young, E. T. 1989, ApJ, 340, 823

Wilner, D. J., & Welch, W. J. 1994, ApJ, 427, 898

Wu, Y., Lin, L., Liu, X., et al. 2019, A&A, 627, A162

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

Yen, H.-W., Takakuwa, S., & Ohashi, N. 2011, ApJ, 742, 57

Yen, H.-W., Takakuwa, S., Ohashi, N., & Ho, P. T. P. 2013, ApJ, 772, 22

Yen, H.-W., Takakuwa, S., Ohashi, N., et al. 2014, ApJ, 793, 1

Zhao, B., Caselli, P., Li, Z.-Y., & Krasnopolsky, R. 2018, MNRAS, 473, 4868

Zucker, C., Speagle, J. S., Schlafly, E. F., et al. 2019, ApJ, 879, 125