OPEN ACCESS

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

CrossMark

# Probing Velocity Structures of Protostellar Envelopes: Infalling and Rotating Envelopes within Turbulent Dense Cores

Jinshi Sai (Insa Choi)<fnref n="1" />,<fnref n="2" />, ∗, Nagayoshi Ohashi<fnref n="2" />, †, Hsi-Wei Yen<fnref n="2" />, ‡, Anaëlle J. Maury<fnref n="3" />,<fnref n="4" />, §, and Sébastien Maret<fnref n="5" />, ¶

Received 2022 September 22; revised 2023 January 14; accepted 2023 January 15; published 2023 February 28

## Abstract

We have observed the three low-mass protostars, IRAS 15398–3359, L1527 IRS, and TMC-1A, with the ALMA 12 m array, the ACA 7 m array, and the IRAM-30 m and APEX telescopes in the C¹⁸O J = 2–1 emission. Overall, the C¹⁸O emission shows clear velocity gradients at radii of ~100–1000 au, which likely originate from the rotation of envelopes, while velocity gradients are less clear and velocity structures are more perturbed on scales of ~1000–10,000 au. IRAS 15398–3359 and L1527 IRS show a break at radii of ~1200 and ~1700 au in the radial profile of the peak velocity, respectively. The peak velocity is proportional to $r^{-1.38}$ or $r^{-1.7}$ within the break radius, which can be interpreted as indicating the rotational motion of the envelope with a degree of contamination by gas motions on larger spatial scales. The peak velocity follows $v_{peak} \propto r^{0.68}$ or $v_{peak} \propto r^{0.46}$ outside the break radius, which is similar to the J/M–R relation of dense cores. TMC-1A exhibits a radial profile of the peak velocity that is not consistent with the rotational motion of the envelope nor the J/M–R relation. The origin of the relation of $v_{peak} \propto r^{0.46}-r^{0.68}$ is investigated by examining correlations of the velocity deviation ($\delta v$) and the spatial scale ($\tau$) in the two sources. The obtained spatial correlations, $\delta v \propto \tau^{-0.6}$, are consistent with the scaling law predicted by turbulence models, which may suggest that large-scale velocity structures originate from turbulence.

Unified Astronomy Thesaurus concepts: Star formation (1569); Low mass stars (2050); Interstellar medium (847); Protostars (1302); Young stellar objects (1834); Radio astronomy (1338); Millimeter astronomy (1061)

## 1. Introduction

Proplanetary disks are sites of planet formation and are ubiquitously found around T Tauri stars (e.g., Guilloteau & Dutrey 1998). However, the detail of their formation process is still poorly understood. Circumstellar disks are expected to be formed around protostars as the angular momentum of dense cores is transferred to the center during the gravitational collapse of these cores (Terebey et al. 1984). Revealing the gas kinematics around protostars on scales ranging from dense cores to disks is thus essential to understand the physical processes of disk formation.

The gas kinematics around protostars have been observationally investigated on various spatial scales. Earlier observational works on Class 0 and I protostars, especially with radio interferometers, have revealed that infalling and rotational motions are dominant at scales of hundreds to thousands of astronomical units (Hayashi et al. 1993; Ohashi et al. 1997a; Momose et al. 1998; Di Francesco et al. 2001; Hogerheijde 2001; Belloche et al. 2002; Arce & Sargent 2004). Recent observations at high angular resolutions suggest that radial distributions of rotational velocity at radii of ~100–1000 au typically follow $v \propto r^{-1}$, i.e., the specific angular momentum $j$ is constant (e.g., Yen et al. 2013; Harsono et al. 2014; Ohashi et al. 2014), which is interpreted as a consequence of the conservation of angular momentum during core collapse (Takahashi et al. 2016). On the other hand, dense cores exhibit

different velocity structures at scales of a few to tens of thousands of astronomical units. Goodman et al. (1993) have measured the mean specific angular momentum $j=J/M$ of each dense core, where $J$ and $M$ are the total angular momentum and mass of the dense cores, respectively, from velocity gradients across subparsec scales assuming rigid-body rotation. They found that $J/M$ is proportional to $R^{1.6}$, where $R$ is the core size. This is known as the J/M–R relation (see also Caselli et al. 2002; Tatematsu et al. 2016; Punanova et al. 2018). Recently, Pineda et al. (2019) revealed that the radial distributions of the specific angular momentum at radii of ~800–10,000 au in protostellar and prestellar sources follow $j \propto r^{1.8}$, which is similar to the J/M–R relation found in the earlier works.

Such observational works on different spatial scales have suggested that velocity structures at inner radii of ~100–1000 au and outer radii of ~1000–10,000 au of protostellar dense cores are distinct: The specific angular momentum $j$ is almost constant on the smaller scale but increases with increasing radius on the larger scale (Ohashi et al. 1997b; Belloche 2013). Gaudel et al. (2020) detected a transition between the two distinct regimes, i.e., the $j$-constant and $j$-increase regimes, at a radius of ~1600 au for the first time by taking the average of the radial distributions of the specific angular momentum of 12 Class 0 protostars. Sai et al. (2022) reported a similar transition between the two regimes at a radius of ~2900 au around the Class I protostar L1489 IRS, which is about two times larger than the averaged transitional radius reported by Gaudel et al. (2020).

Such a difference found in the two works raises a question of how different the transitional radius between the two regimes is

<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, Sec. 4, Roosevelt Road, Taipei 10617, Taiwan; ohashi@asiaa.sinica.edu.tw</footnote>

<footnote n="3">³AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Citè, F-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, F-38000 Grenoble, France</footnote>
<table><caption>**Table 1**
Source Properties</caption><thead><tr><th>Source</th><th>R.A.<br>(J2000)</th><th>Decl.<br>(J2000)</th><th>T<sub>bol</sub><br>(K)</th><th>L<sub>bol</sub><br>(L<sub>☉</sub>)</th><th>d<br>(pc)</th><th>M<sub>*</sub><br>(M<sub>☉</sub>)</th><th>R<sub>disk</sub><br>(au)</th><th>i<br>(°)</th><th>References</th></tr></thead><tbody><tr><td>IRAS 15398-3359</td><td>15:43:02.24</td><td>-34:09:06.81</td><td>61</td><td>0.92</td><td>150</td><td>0.007</td><td>≥40</td><td>70</td><td>1, 2, 3, 4, 5</td></tr><tr><td>L1527 IRS</td><td>4:39:53.88</td><td>+26:03:09.55</td><td>44</td><td>2.0</td><td>140</td><td>0.45</td><td>74</td><td>85-90</td><td>6, 7, 8, 9</td></tr><tr><td>TMC-1A</td><td>4:39:35.20</td><td>+25:41:44.35</td><td>118</td><td>2.7</td><td>140</td><td>0.68</td><td>100</td><td>55-65</td><td>6, 7, 10, 11</td></tr></tbody></table>

**References.** (1) Froebrich (2005); (2) Comerón (2008); (3) Okoda et al. (2018); (4) Oya et al. (2014); (5) Yen et al. (2017); (6) Kristensen et al. (2012); (7) Zucker et al. (2019); (8) Aso et al. (2017); (9) Tobin et al. (2008); (10) Aso et al. (2015); (11) Harsono et al. (2014).

<table><caption>**Table 2**
Summary of ALMA Observations</caption><thead><tr><th>Source</th><th>R.A.<sup>a</sup><br>(J2000)</th><th>Decl.<sup>a</sup><br>(J2000)</th><th>Baseline Length<sup>b</sup><br>(kλ)</th><th>Date<br>(UT)</th><th>Calibrators<br>(Gain, Bandpass, Flux)</th><th>Δν<br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>IRAS 15398-3359</td><td>15:43:02.16</td><td>-34:09:09.0</td><td>13-478</td><td>2014 Apr 30, May 19, Jun 6</td><td>J1534-3526, J1427-4206, Titan</td><td>0.084</td></tr><tr><td rowspan="2">L1527 IRS</td><td rowspan="2">04:39:53.91</td><td rowspan="2">+25:41:44.4</td><td>12-363</td><td>2015 May 24</td><td>J0510+1800, J0423-0120, J0510+180</td><td>0.042</td></tr><tr><td>21-1499</td><td>2015 Sep 20, 2015 May 24</td><td>J0429+2724, J0510+1800, J0423-013</td><td>0.042</td></tr><tr><td rowspan="2">TMC-1A</td><td rowspan="2">04:39:35.20</td><td rowspan="2">+25:41:44.4</td><td>12-363</td><td>2015 May 24</td><td>J0510+1800, J0423-0120, J0510+180</td><td>0.042</td></tr><tr><td>21-1499</td><td>2015 Sep 20</td><td>J0429+2724, J0510+1800, J0423-013</td><td>0.042</td></tr></tbody></table>

**Notes.**

<sup>a</sup> R.A. and decl. of the phase center of the observations.

<sup>b</sup> Projected baseline lengths at 220 GHz.

from one source to another. Furthermore, the origin of the velocity structure at the *j*-increase regime is not well understood. Observations have shown complex velocity structures of dense cores on ~1000–10,000 au scales: Local velocity gradients within the dense cores show a significant variation in direction. Such complex velocity structures cannot be interpreted as simple rigid-body rotation (Caselli et al. 2002; Tobin et al. 2011; Chen et al. 2019b). On the one hand, simulation works suggest that the *J/M-R* relation found in dense cores is inherited from turbulence in the parental filaments or clouds (Chen & Ostriker 2018; Misugi et al. 2019). The similarity between the *J/M-R* relation (v~j/r~ ∝ r<sup>0.6</sup>) and the scaling law of the cloud-scale turbulence (δv ∝ τ<sup>0.5</sup>, where δv is the velocity deviation and τ is the spatial scale; Larson 1981; McKee & Ostriker 2007) has also been pointed out (Tatematsu et al. 2016; Chen et al. 2019a; Pineda et al. 2019; Gaudel et al. 2020). On the other hand, it is proposed that the *J/M-R* relation arises as a consequence of gravitational contraction with angular momentum loss via turbulent viscosity (Arroyo-Chávez & Vázquez-Semadeni 2022).

Turbulence possesses scale-dependent correlations of fluid variables on the spatial scale. This property of the turbulence is observed as a correlation between the velocity deviation δv and the spatial scale τ. A number of works have studied the correlation in molecular clouds with various methods and reported the scaling law of δv ∝ τ<sup>γ</sup>, where γ ~ 0.3–0.5, on ~0.03–30 pc scales (e.g., Larson 1981; Heyer & Brunt 2004). These indices are often interpreted as indicative of Kolmogorov turbulence, which yields γ = 1/3 for incompressible fluids (Kolmogorov 1941), or Burgers turbulence, which yields γ = 1/2 for highly compressible fluids (Burgers 1974). However, such measurements have not been performed at smaller scales of a few to several thousand astronomical units inside dense cores because of the limit of angular and velocity resolutions.

In this paper, we present observations toward three well-studied protostars, i.e., two Class 0 protostars IRAS 15398–3359 and L1527 IRS (IRAS 04365+2557) and one Class I protostar TMC-1A (IRAS 04365+2535), with the ALMA 12 m array, the ACA 7 m array, and the IRAM-30 m and Atacama Pathfinder Experiment (APEX) telescopes in the C<sup>18</sup>O *J* = 2–1 line. We investigate velocity structures around the protostars over a radius of ~100–10,000 au with the C<sup>18</sup>O maps. All three protostars are located in the nearby star-forming regions: IRAS 15398–3359 is in the Lupus I molecular cloud (d ~ 150 pc; Comerón 2008), and L1527 IRS and TMC-1A are in the Taurus molecular cloud (d ~ 140 pc; Zucker et al. 2019). Previous works have reported that the rotational velocity of their envelopes is approximately proportional to r<sup>-1</sup> (Yen et al. 2013, 2017; Ohashi et al. 2014; Aso et al. 2015, 2017). The Keplerian rotation curves of their disks were also identified from line observations, and dynamical masses were estimated (Tobin et al. 2012b; Ohashi et al. 2014; Aso et al. 2017; Okoda et al. 2018; Maret et al. 2020). The source properties are summarized in Table 1.

The outline of this paper is as follows. Observations and data reduction are summarized in Section 2. Observational results are presented in Section 3. Analyses of the velocity structures are provided in Section 4. The possible interpretation and implication of analysis on the velocity structures are discussed in Section 5. Finally, all results and discussions are summarized in Section 6.

## 2. Observations

### 2.1. ALMA 12 m Array Observations

The three protostars were observed with the Atacama Large Millimeter/submillimeter Array (ALMA) 12 m array in the C<sup>18</sup>O *J* = 2–1 (219.560358 GHz; *E*<sub>2</sub> = 15.8 K; *A*<sub>21</sub> = 6.01 × 10<sup>-7</sup> s<sup>-1</sup>) line emission in ALMA Cycle 2. The ALMA observations are summarized in Table 2. The observations consist of a single field for each source, and the field of view is
<table><caption>**Table 3**
Summary of ACA 7 m Array Observations</caption><thead><tr><th>Source</th><th>Number of Fields</th><th>Date (UT)</th><th>Gain Calibrators</th><th>Bandpass and Flux Calibrators<sup>a</sup></th><th>Δν (km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>IRAS 15398-3359</td><td>28</td><td>2019 Dec 19</td><td>J1534-3526</td><td>J1337-1257</td><td>0.084</td></tr><tr><td>L1527 IRS</td><td>27</td><td>2019 Dec 1-17</td><td>J0426+2327</td><td>J0423-0120/J0725-0054</td><td>0.084</td></tr><tr><td>TMC-1A</td><td>25</td><td>2019 Dec 1-17</td><td>J0426+2327</td><td>J0423-0120/J0725-0054</td><td>0.084</td></tr></tbody></table>

**Note.**

<sup>a</sup> Bandpass and flux were calibrated using the same source, either J1337-1257, J0423-0120, or J0725-0054 depending on the target source and observing date.

26" at 220 GHz. The observations of L1527 IRS and TMC-1A were conducted with a compact configuration with 34 antennas and an extended configuration with 35 antennas. The projected baseline lengths of compact and extended configurations were 12–363 kλ (16–494 m) and 21–1499 kλ (29–2043 m) at 220 GHz, respectively. The spectral window for the C<sup>18</sup>O J = 2–1 line had a bandwidth of 58.6 MHz and a spectral resolution of 30.5 kHz, corresponding to a velocity resolution of 0.042 km s<sup>-1</sup>. The total on-source time was 48 minutes for L1527 IRS and 46 minutes for TMC-1A. The data obtained with the compact and extended antenna configurations were calibrated with the Common Astronomy Software Applications package (CASA; McMullin et al. 2007) version 4.5.0 and 4.5.2, respectively. IRAS 15398–3359 was observed with a single configuration with 34 antennas. The shortest and longest projected baseline lengths were 13 kλ and 478 kλ (18 and 651 m) at 220 GHz, respectively. The spectral window for the C<sup>18</sup>O J = 2–1 line had a bandwidth of 58.6 MHz and a spectral resolution of 61.0 kHz (~0.084 km s<sup>-1</sup>). The on-source time was 88 minutes. The details of the observations are summarized in Yen et al. (2017). The data were calibrated with CASA 4.3.1.

## 2.2. ACA 7 m Array Mosaic Observations

We have conducted mosaic observations of the three protostars in the C<sup>18</sup>O J = 2–1 emission using the 7 m array of the Atacama Compact Array (ACA) during 2019 December 1–19 in the ALMA Cycle 7 with 10 antennas. The summary of the observations is presented in Table 3. The mosaic observations covered ~2' × 2' regions centered at protostellar positions. The spectral window for the C<sup>18</sup>O J = 2–1 line had a bandwidth of 120.0 MHz and a spectral resolution of 61.0 kHz, providing a velocity resolution of 0.084 km s<sup>-1</sup>. The mosaic maps of IRAS 15398, L1527 IRS, and TMC-1A consist of 28, 27, and 25 fields, and the on-source time on each single field was 1.7, 1.1, and 3.2 minutes, respectively. The projected baseline lengths ranged from 5 kλ to 28 kλ (6.8–38 m) for L1527 IRS and TMC-1A, and from 6 kλ to 35 kλ (8.2–48 m) for IRAS 15398–3359 at 220 GHz. The data were calibrated in the pipeline with CASA 5.6.1 and its pipeline version was 42866M (Pipeline-CASA56-P1-B).

## 2.3. Single-dish Observations

Observations of the three protostars with the IRAM-30 m and APEX telescopes have been carried out in the C<sup>18</sup>O J = 2–1 emission in the on-the-fly (OTF) mapping mode to cover 2' square regions. The position-switching method was used for the OTF mapping in all observations. The summary of the observations is presented in Table 4.

L1527 IRS was observed with the IRAM-30 m telescope using the Heterodyne Receiver Array (HERA) receiver and the Versatile Spectrometer Arrays (VESPA) backend with a spectral resolution of 20 kHz (0.027 km s<sup>-1</sup>). The pointing was corrected every two hours with quasars near the target: 0316+413 and 0415+379. The focus was corrected every 4 hr with the quasar 0316+413. The reference position for the position switching was Δα = 6', Δδ = -12' with respect to the source position given in Table 1.

TMC-1A was observed with the IRAM-30 m telescope using the heterodyne receiver Eight MIxer Receiver (EMIR) E230 with the VESPA backend with a spectral resolution of 20 kHz (0.027 km s<sup>-1</sup>). The telescope pointing and focus were corrected every ~2–3 and ~4–5 hr with the quasars 0430 +052 and 0316+413, respectively. The reference position for position switching was Δα = 32', Δδ = 48' with respect to the source position.

Observations of IRAS 15398–3359 were conducted using the APEX telescope on 2017 July 31 to August 2. The APEX-1 heterodyne receiver was used together with the RPG extended bandwidth fast Fourier transform spectrometer (XFFT5) backend. The spectral resolution was 76 kHz (0.1 km s<sup>-1</sup>). The pointing was corrected every ~1–2 hr with RAFGL 4211 and IRAS 15194–5115. The reference position for the position switching was Δα = 6', Δδ = 19' with respect to the protostellar position.

All single-dish data were reduced with the GILDAS software package. The antenna temperature $T_A^*$ was converted to the main-beam temperature $T_{mb}$ according to the relation $T_{mb} = T_A^* F_{eff}/B_{eff}$, where $B_{eff}$ is the main-beam efficiency and $F_{eff}$ is the forward efficiency. We adopted $B_{eff} = 0.60$ and $F_{eff} = 0.92$ for the C<sup>18</sup>O J = 2–1 data obtained with IRAM-30 m, and $B_{eff} = 0.63$ and $F_{eff} = 0.95$ for the APEX data. The half-power beamwidth (HPBW) at the rest frequency of C<sup>18</sup>O J = 2–1 is ~12" for the IRAM-30 m maps and ~30" for the APEX map.

## 2.4. Combined Maps

The interferometric and single-dish data were combined for further analysis with the task feather in CASA 5.6.1, which combines interferometric and single-dish images in Fourier space, weighting them by the spatial frequency response of each image. Two maps with different map sizes and angular resolutions were generated for each source to trace different spatial scales: (i) combined maps of the 12 m array, 7 m array, and single-dish data at high angular resolutions to resolve structures near protostars (within a radius of ~12") and (ii) combined maps of the 7 m array and single-dish data covering wider areas around protostars (~2' × 2'). We compared fluxes of the data in both the image and visibility domains and confirmed that all data are consistent in flux. To compare the
<table><caption>**Table 4**
Summary of Single-dish Observations</caption><thead><tr><th>Source</th><th>Telescope</th><th>Receiver/Backend</th><th>Date (UT)</th><th>Δv (km s<sup>-1</sup>)</th><th>Mapping Area</th></tr></thead><tbody><tr><td>IRAS 15398–3359</td><td>APEX</td><td>APEX-1/XFFTS</td><td>2017 Jul 31–Aug 2</td><td>0.01</td><td>4′.3 × 4′.3</td></tr><tr><td>L1527 IRS</td><td>IRAM-30 m</td><td>HERA/VESPA</td><td>2014 Sep 3–4</td><td>0.027</td><td>2′.2 × 2′.2</td></tr><tr><td>TMC-1A</td><td>IRAM-30 m</td><td>EMIR/VESPA</td><td>2019 Dec 27–31</td><td>0.027</td><td>2′.7 × 2′.5</td></tr></tbody></table>

<table><caption>**Table 5**
Summary of C<sup>18</sup>O Maps</caption><thead><tr><th>Map</th><th></th><th>IRAS 15398–3359</th><th>L1527 IRS</th><th>TMC-1A</th></tr></thead><tbody><tr><td rowspan="3">Small scale</td><td>Angular resolution</td><td>0″77 × 0″72 (75°)</td><td>0″94 × 0″82 (8°1)</td><td>0″93 × 0″81 (7°6)</td></tr><tr><td>rms (mJy beam<sup>-1</sup>)</td><td>5.2</td><td>12</td><td>12</td></tr><tr><td>Velocity resolution (km s<sup>-1</sup>)</td><td>0.11</td><td>0.084</td><td>0.084</td></tr><tr><td rowspan="3">Large scale</td><td>Angular resolution</td><td>8″1 × 4″4 (85°)</td><td>7″5 × 6″4 (−1°5)</td><td>7″7 × 6″4 (−85°)</td></tr><tr><td>rms (mJy beam<sup>-1</sup>)</td><td>160</td><td>290</td><td>300</td></tr><tr><td>Velocity resolution (km s<sup>-1</sup>)</td><td>0.11</td><td>0.084</td><td>0.084</td></tr></tbody></table>

single-dish data in the visibility domain, pseudo-visibilities were generated using the MIRIAD software following the methods described in Takakuwa et al. (2003). The peak positions of the continuum emission observed with ALMA in previous works were adopted as protostellar positions and map centers (Aso et al. 2015, 2017; Yen et al. 2017). The protostellar positions are summarized in Table 1. All imaging processes were performed using CASA 5.6.1. The angular resolutions, rms, and velocity resolutions of the produced maps are summarized in Table 5.

### 2.4.1. Small-scale Maps

Maps covering regions near protostars at high angular resolutions were produced by combining the 12 m array, 7 m array, and single-dish data. We have referred to these maps as small-scale maps throughout this paper. First, the combined maps of 12 m and 7 m array data were produced through joint deconvolutions with the task `tclean`. To enhance the sensitivity to extended structures, the natural weighting and taper with the FWHM of 200 kλ were adopted. Multiscale CLEAN was applied with scale sizes of zero, 1", and 3". Then, the 12 m+7 m maps were combined with the single-dish maps using the following steps: (1) The single-dish and 12 m+7 m images were trimmed to exclude masked regions at the map edges, (2) the single-dish images were multiplied by the 12 m +7 m primary beam response, (3) the modified single-dish maps were combined with the 12 m+7 m maps with the task `feather`, and (4) the combined maps were corrected by the 12 m+7 m primary beam response. As such, the angular resolution, rms noise level, and velocity resolution are ~0"8, ~11 mJy beam<sup>-1</sup>, and 0.09 km s<sup>-1</sup>, respectively, on average. The final map size centered at the protostellar positions is ~24 arcsec<sup>2</sup>.

### 2.4.2. Large-scale Maps

Combined maps of the 7 m array and single-dish data were produced to cover wider regions. We have referred to these maps as large-scale maps throughout this paper. First, mosaic maps of the 7 m array data were produced using the `tclean` task with a Briggs robust parameter of 0.5 and the Hogbom deconvolver. The cleaned 7 m maps were combined with the single-dish maps using the `feather` task and following the

same steps used to produce the small-scale maps. The angular resolution, rms noise level, and velocity resolution of the final maps are ~7", ~25 mJy beam<sup>-1</sup>, and 0.09 km s<sup>-1</sup>, respectively, on average. The maps cover ~2' square regions centered at the protostars.

## 3. Results

The C<sup>18</sup>O *J* = 2–1 integrated intensity and centroid velocity maps of the three protostars are presented in Figure 1, where small-scale and large-scale maps are shown in the left and right columns, respectively. For the integrated intensity maps, the C<sup>18</sup>O emission was integrated over velocity ranges where the emission was detected at least above 3σ. For the centroid velocity maps, only emission detected above 5σ was considered. In both small- and large-scale maps, the emission peaks are located at the protostellar positions. Thus, the C<sup>18</sup>O emission traces the envelopes associated with the protostars. Figures A1–A6 in Appendix A show velocity channel maps for each protostar, and more details of moment 0 maps are presented below.

### 3.1. IRAS 15398–3359

In the small-scale map in Figure 1(a), the distribution of the C<sup>18</sup>O integrated intensity exhibits an X-shape morphology: a part of the emission is elongated from northeast to southwest, while another part is along a direction from southeast to northwest. These elongated structures are ~15" in length across the protostellar position. The elongation from the northeast to the southwest is the same as the direction of the primary outflow associated with the protostar (Oya et al. 2014; Yen et al. 2017), and it exhibits a velocity gradient along the elongation. The velocity gradient showing redshifted velocity on its northeastern side and blueshifted velocity on its southwestern side is also the same as that of the primary outflow. These suggest that the emission elongated from the northeast to the southwest is likely affected by the primary outflow. Weak emission extends to the southeast from the protostar, which is almost the same direction as that of the secondary outflow suggested in a previous work (Okoda et al. 2021). Outside the X-shape morphology, the C<sup>18</sup>O emission extends over the entire map without any clear velocity gradient. In the large-scale map shown in Figure 1(b), a weak velocity
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Integrated intensity (contour) and centroid velocity (color) maps of the C¹⁸O J = 2–1 emission of IRAS 15398–3359, L1527 IRS, and TMC-1A. The left and right columns show the small- and large-scale maps, respectively. In the maps of the left column, contour levels are from 5σ to 20σ in steps of 5σ, and then from 20σ to 80σ in steps of 20σ, where 1σ is the rms noise of a map for a given source, which is summarized in Table 5. In the maps of the right column, the contour steps are the same as those in the maps of the left column but begin from 40σ in the map of IRAS 15398–3359 and from 10σ in the maps of L1527 IRS and TMC-1A. Red and blue arrows with solid lines show directions of (primary) outflows. Blue arrows with dashed lines indicate the direction of the secondary outflow of IRAS 15398–3359. Black crosses at the center and ellipses in the bottom-left corners denote the protostellar positions and the synthesized beam size, respectively.</figcaption></figure>

gradient is seen from northwest to southeast, although the velocity structure appears more perturbed than in the small-scale map.

While the systemic velocity of this system has been estimated based on the rotational velocity of the envelope and disk at radii less than 100 au (Yen et al. 2017; Okoda et al. 2018), the values reported in previous works may not be appropriate at the spatial scales of the current maps because the center of mass could be different due to the very small

protostellar mass of this source (~0.007 $M_\odot$; Okoda et al. 2018). Indeed, Okoda et al. (2018) reported a discrepancy between systemic velocities estimated from the rotational velocity of the disk and the envelope. Thus, we estimate the systemic velocity to be 5.183 ± 0.005 km s⁻¹ using our large-scale map by fitting a Gaussian function to the C¹⁸O spectrum at the protostellar position. We have adopted this centroid velocity of 5.18 km s⁻¹ for the systemic velocity of this system in the following analyses and discussions.
<table><caption>**Table 6**
Results of the Two-dimensional Linear Fitting to the Centroid Velocity Maps</caption><thead><tr><th rowspan="2">r<sub>fit</sub><br>(")</th><th colspan="3">IRAS 15398–3359</th><th colspan="3">L1527 IRS</th><th colspan="3">TMC-1A</th></tr><tr><th>v<sub>0</sub><sup>a</sup><br>(km s<sup>-1</sup>)</th><th>G<br>(km s<sup>-1</sup> pc<sup>-1</sup>)</th><th>θ<br>(°)</th><th>v<sub>0</sub><br>(km s<sup>-1</sup>)</th><th>G<br>(km s<sup>-1</sup> pc<sup>-1</sup>)</th><th>θ<br>(°)</th><th>v<sub>0</sub><br>(km s<sup>-1</sup>)</th><th>G<br>(km s<sup>-1</sup> pc<sup>-1</sup>)</th><th>θ<br>(°)</th></tr></thead><tbody><tr><td>5</td><td>5.18</td><td>14.1 ± 0.2</td><td>69.3 ± 0.8</td><td>5.83</td><td>41.7 ± 0.2</td><td>14.0 ± 0.3</td><td>6.34</td><td>107.8 ± 0.2</td><td>-146.18 ± 0.09</td></tr><tr><td>10</td><td>5.21</td><td>1.05 ± 0.05</td><td>67 ± 3</td><td>5.87</td><td>8.39 ± 0.04</td><td>5.8 ± 0.3</td><td>6.39</td><td>18.50 ± 0.03</td><td>-125.9 ± 0.1</td></tr><tr><td>20</td><td>5.23</td><td>1.6 ± 0.1</td><td>83 ± 4</td><td>5.91</td><td>3.7 ± 0.1</td><td>96 ± 2</td><td>6.40</td><td>8.78 ± 0.08</td><td>-149.6 ± 0.5</td></tr><tr><td>30</td><td>5.22</td><td>0.97 ± 0.05</td><td>96 ± 3</td><td>5.93</td><td>2.26 ± 0.04</td><td>108 ± 1</td><td>6.38</td><td>4.17 ± 0.03</td><td>-144.3 ± 0.4</td></tr><tr><td>40</td><td>5.22</td><td>0.89 ± 0.03</td><td>107 ± 2</td><td>5.94</td><td>2.66 ± 0.02</td><td>112.9 ± 0.5</td><td>6.37</td><td>2.49 ± 0.02</td><td>-135.5 ± 0.4</td></tr><tr><td>50</td><td>5.22</td><td>0.71 ± 0.02</td><td>137 ± 1</td><td>5.94</td><td>2.56 ± 0.01</td><td>113.7 ± 0.3</td><td>6.36</td><td>1.72 ± 0.01</td><td>-129.2 ± 0.4</td></tr><tr><td>60</td><td>5.22</td><td>0.73 ± 0.01</td><td>166 ± 1</td><td>5.95</td><td>2.540 ± 0.009</td><td>115.7 ± 0.2</td><td>6.36</td><td>1.356 ± 0.008</td><td>-125.1 ± 0.3</td></tr></tbody></table>

**Note.**

## 3.2. L1527 IRS

Figures 1(c) and (d) show the integrated intensity and centroid velocity maps of L1527 IRS. The small-scale map in Figure 1(c) shows a flattened structure along the north–south direction within a radius of ~2" around the protostar. A clear velocity gradient is seen along the flattened structure. The velocity gradient consists of velocity components blueshifted and redshifted with respect to the systemic velocity of the protostar (5.8 km s⁻¹; Aso et al. 2017) to the south and north sides of the protostar, respectively. This feature suggests that the emission likely traces the rotational motion of the disk and envelope, as found in previous works (Tobin et al. 2012b; Ohashi et al. 2014; Sakai et al. 2014; Aso et al. 2017). Weak emission is elongated from northeast to southwest over ~20" across the protostellar position outside the flattened emission. Another component is also extended from the protostar toward the southeast. These structures seem to be associated with its outflow cavity, as the outflow has been detected along the east–west direction (Tamura et al. 1996; Tobin et al. 2008). In the large-scale map in Figure 1(d), the emission is elongated and exhibits a weak velocity gradient in the east–west direction, which likely traces the outflow (see also velocity channel maps in Figure A5). On the other hand, the velocity gradient in the direction perpendicular to the outflow becomes less clear in the large-scale map.

## 3.3. TMC-1A

In the small-scale map of TMC-1A in Figure 1(e), the C¹⁸O emission near the protostar has a flattened shape and a clear velocity gradient in the direction from northeast to southwest, which is almost perpendicular to the outflow. The velocity gradient consists of blueshifted and redshifted velocities with respect to the systemic velocity (6.4 km s⁻¹; Aso et al. 2015) on the northeast and southwest sides of the protostar, respectively. These features are consistent with previous ALMA observations, which suggested that these velocity structures originate from the rotating disk and envelope (Aso et al. 2015). Outside the flattened emission, weak emission is extended from northeast to southwest and from southeast to northwest. These components are considered to trace its outflow cavity, as the outflow was detected along the southeast to the northwest (Aso et al. 2015). A weak velocity gradient from the northeast to the southwest is observed in the large-scale map in Figure 1(f), although it appears more perturbed.

## 4. Analysis

### 4.1. Two-dimensional Velocity Gradient

To characterize the observed velocity gradients quantitatively, the following two-dimensional linear function, introduced by Goodman et al. (1993), was fitted to the centroid velocity maps with different aperture areas representing different spatial scales:

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

where Δα and Δδ are the offsets from the protostellar position in R.A. α and decl. δ, respectively, and v₀, a, and b are the fitting parameters. The magnitude and position angle of the velocity gradient are calculated as follows:

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

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

where G is the magnitude, θ is the position angle measured from north to east, and d is the distance of the source. The small-scale maps were used for the fitting with the aperture radius r<sub>fit</sub> of 5"–10", while the fitting with r<sub>fit</sub> of 20"–60" was performed using the large-scale maps.

The fitting results are summarized in Table 6. IRAS 15398–3359 exhibits the weakest velocity gradient among the three protostars on all scales. TMC-1A shows the strongest velocity gradient near the protostar within a radius of 30", while L1527 IRS does on scales larger than 30". Overall, the velocity gradients are stronger near the protostars and weaker at larger radii in all sources. The same fitting has been performed in L1527 IRS by Maret et al. (2020) and Gaudel et al. (2020) using different data sets. Gaudel et al. (2020) report G = 66 km s⁻¹ pc⁻¹ and θ = 22° within a radius of 5" using a C¹⁸O J = 2–1 map of PdBI+IRAM-30 m data, which is slightly different from the values measured in the current work. These differences arise because our map has better sensitivity to extended emission, and more pixels with low velocities close to the systemic velocity are included even in the fitting within a radius of 5". On the other hand, the centroid velocity map presented by Gaudel et al. (2020) has blank pixels at the map edge because of the sensitivity limit, excluding extended, low-velocity components from the fitting. They measure G = 2 km s⁻¹ pc⁻¹ and θ = 113° within a radius of 40" using an IRAM-30 m C¹⁸O J = 2–1 map, which is in agreement with our results. Maret et al. (2020) present G = 671 km s⁻¹ pc⁻¹ within a radius of 2", which also agrees

<footnote><sup>a</sup> Fitting errors of v₀ are less than 0.1%.</footnote>
with the trend that the velocity gradient is stronger at smaller radii.

Scale dependence of the directions of the velocity gradients differs from source to source. The position angles of the velocity gradients of IRAS 15398–3359 within a radius of 10" is ~70°, which is a direction of the associated outflow (Oya et al. 2014; Yen et al. 2017). The directions of the velocity gradients of IRAS 15398–3359 vary significantly from 70° to 170° as the scale becomes larger. In L1527 IRS, the position angle of the velocity gradient within a radius of 5" is ~14°, which is close to the direction of the disk major axis (P. A. = 0°; Tobin et al. 2012a; Ohashi et al. 2014; Aso et al. 2017). L1527 IRS also shows a significant change in the directions of the velocity gradients by ~90° as the scale becomes larger. The velocity gradient over a radius of 60" is along P.A. ~110°, which is close to the outflow direction (P. A. = 90°; Hogerheijde et al. 1998; Tobin et al. 2008; Oya et al. 2015). On the contrary, TMC-1A shows a small dispersion of the directions of the velocity gradients ~±10° over spatial scales of 5"–60". The directions of the velocity gradients are roughly normal to the outflow with P. A. = -25° to -17° (Hogerheijde et al. 1998; Aso et al. 2015).

## 4.2. Peak Velocity Measurements

The small-scale maps of L1527 IRS and TMC-1A exhibit velocity gradients, which likely trace the rotational motion of the disks and envelopes, although such a velocity gradient due to the rotational motion of the disk and envelope is not very clear in IRAS 15398–3359 because of its small rotational velocity compared to those in other sources (Yen et al. 2017). On the other hand, velocity gradients are less clear and more perturbed in the large-scale maps of all three sources. In this subsection, we measure the peak velocity as a function of radius along the disk major axis and examine its radial dependence to investigate how the nature of the velocity structure changes with radius.

For this purpose, position–velocity (PV) diagrams cut along the disk major axis were produced from both small- and large-scale maps. The velocity component along a disk major axis mainly reflects rotational velocity and is less affected by outflows. Therefore, a PV cut along the disk major axis is reasonable to trace rotation motion at inner radii and investigate how it changes with changing radii. Position angles of 140°, 0°, and 70° are adopted for the disk major axes of IRAS 15398–3359, L1527 IRS, and TMC-1A, respectively, based on previous observations of their disks (Tobin et al. 2012b; Ohashi et al. 2014; Oya et al. 2014; Aso et al. 2015, 2017; Yen et al. 2017). The width of the PV cuts is one pixel (~1/10 beam size), which means no averaging pixels perpendicular to the slices.

The generated PV diagrams are presented in Figure 2. Small-scale PV diagrams of all three protostars, shown in Figures 2(a), (c), and (e), exhibit high-velocity components at small offsets within ~±2". The emission peak at each position appears at a higher velocity with respect to the systemic velocity as the offset decreases at offsets of ~3"–8" in L1527 IRS and TMC-1A, although this feature is not very clear in IRAS 15398–3359. These velocity structures can be interpreted as the differential rotation of the envelope, where the rotational velocity increases with decreasing radius. In addition to these velocity structures, the small-scale PV diagram of IRAS 15398–3359 shows an emission whose velocity appears

to linearly increase with the radius at offsets of ~-6 to -3. This velocity structure was also reported in previous ALMA observations by Okoda et al. (2021), who suggest the velocity component originates from the secondary outflow launched in a direction perpendicular to the primary outflow. On the other hand, the velocity gradient is less clear, and the differential rotation feature is not seen at offsets of ≥12" and ≲12" in the large-scale PV diagrams presented in Figures 2(b), (d), and (f). Moreover, the emission peak at each position at those offsets appears at a redshifted velocity on both the north and south sides in L1527 IRS, which would not be explained by rotational motion.

In order to measure the peak velocity as a function of radius, the method used in Sai et al. (2020, 2022) is adopted, i.e., a Gaussian function is fitted to the spectrum at each offset in the PV diagrams. The offset coordinates are sampled in steps of half of the beam size, and the spectrum measured on a single pixel at a given offset is used for the fitting. The fitting is performed only on the spectrum whose detected peak intensity is above 6σ. The ±2 velocity channels around the intensity maximum of the spectrum are used for the fitting to better trace the peak velocity of the non-Gaussian shape spectrum (Sai et al. 2020). The peak velocities are measured at offset ranges of 3" to 12" and -12" to -3" using the small-scale PV diagrams. The offset range of -3" to 3" is excluded from the fitting since the emission peaks appear at low velocities, and the peak velocities no longer trace the feature of differential rotation within the offset range. The peak velocities outside 12" are measured using large-scale PV diagrams. The error of the measured peak velocity is derived as the fitting error.

The measured data points are overlaid on the PV diagram in Figure 2. Figure 2(a) shows that the peak velocities of IRAS 15398–3359 measured on small scales are overall blueshifted on the northwest side and redshifted on the southeast side of the protostar with respect to $v_{sys}$, which may suggest rotational motion despite the very small measured velocity. Peak velocities measured with the large-scale PV diagram in Figure 2(b) show a similar feature, while the peak velocity at offsets larger than 40" on the southeast side is blueshifted. The peak velocity is further blueshifted as the radius increases on the northwest side. Figure 2(c) shows that the peak velocities measured at small scales for L1527 IRS roughly trace the feature of differential rotation within offset ranges of 3"–5" on the south side and -7" to -3" on the north side, while the peak velocities show abrupt changes across the systemic velocity outside of these offsets. The peak velocities measured with the large-scale PV diagram shown in Figure 2(d), on the other hand, are redshifted on both the north and south sides. Moreover, they show more redshifted velocities with respect to the systemic velocity at larger radii. The peak velocities measured for TMC-1A also show abrupt changes at an offset of ~-4" on the northeast side and ~5" on the southwest side across the systemic velocity in the small-scale PV diagram, as presented in Figure 2(e). In the large-scale PV diagram shown in Figure 2(f), most of the peak velocities on the northeast and southwest sides are blueshifted and redshifted, respectively, with respect to the systemic velocity.

Figure 3 presents the measured peak velocities in log r–log v diagrams in the form of the relative velocity, $v_{peak} = |v_{LSR} - v_{sys}|$, to investigate their radial dependence. In addition to the measured data points represented by colored circles, averaged peak velocities, which are calculated by averaging the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** PV diagrams of the C$^{18}$O 2–1 emission cut along the disk major axes. Contour levels are from 3σ to 12σ in steps of 3σ and from 12σ to 60σ in steps of 6σ, where 1σ is the rms noise summarized in Table 5. Green and purple markers show the peak velocities measured on opposite sides of the protostars through the Gaussian fitting. Vertical dashed lines denote the systemic velocities. Vertical and horizontal bars in the bottom-left corners indicate the spatial and velocity resolutions, respectively. Scales of the velocity axis are different in the right and left columns.</figcaption></figure>

velocities obtained at a given offset assuming an azimuthally symmetric velocity structure, are also plotted in black circles in the diagrams.

The radial profile of the averaged peak velocity for IRAS 15398–3359 shown in Figure 3(a), overall, exhibits a systematic trend as a function of radius: The averaged peak velocities decrease and then increase with increasing radius, with a transition around a radius of ~1200 au. In addition, the radial profile shows small bumps at radii of ~900 and ~1500 au. The averaged peak velocities of L1527 IRS presented in Figure 3(b) also decrease and then increase with increasing radius with a transition radius of ~1700 au. Similar radial velocity profiles with two such regimes have been observed for

envelopes around 12 Class 0 protostars by calculating an averaged profile among the sources (Gaudel et al. 2020) and also in the envelope around the Class I protostar L1489 IRS (Sai et al. 2022). TMC-1A also exhibits a break at a radius of ~800 au in the radial profile of the averaged peak velocity, as shown in Figure 3(c). However, the measured peak velocities are more scattered than those of the other two sources. Moreover, the peak velocities measured on the southwest and northeast sides of the protostar are mostly inconsistent with each other at small radii of <~800 au.

To characterize the radial profiles showing two regimes, the following double-power-law function was fitted to the averaged peak velocities with the Markov Chain Monte Carlo (MC) code
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Radial profiles of the peak velocities measured with the PV diagrams. The green and purple circles show the data points measured on opposite sides of the protostars, and the black circles represent the averaged value at each radius. The dashed lines are the extrapolated best-fit power-law functions measured at ~100 au in previous works (Aso et al. 2015, 2017; Yen et al. 2017). The black solid lines show the double-power-law functions with the best-fit parameters for the averaged peak velocities in this work. The green solid line shows the double-power-law function with the best-fit parameter for the peak velocity measured on the southwest side of TMC-1A.</figcaption></figure>

<table><caption>**Table 7**

Best-fit Parameters of the Double-power-law Fitting to the Measured Peak Velocities</caption><thead><tr><th>Source</th><th>v<sub>break</sub> (km s<sup>-1</sup>)</th><th>r<sub>break</sub> (au)</th><th>p<sub>in</sub></th><th>p<sub>out</sub></th></tr></thead><tbody><tr><td>IRAS 15398-3359</td><td>0.013 ± 0.001</td><td>1180 ± 90</td><td>-1.7 ± 0.3</td><td>0.68 ± 0.05</td></tr><tr><td>L1527 IRS</td><td>0.069 ± 0.003</td><td>1700 ± 50</td><td>-1.38 ± 0.03</td><td>0.46 ± 0.05</td></tr><tr><td>TMC-1A</td><td>0.068 ± 0.005</td><td>760 ± 10</td><td>-4.1 ± 0.1</td><td>-0.30 ± 0.06</td></tr><tr><td>TMC-1A (SW side only)</td><td>0.099 ± 0.007</td><td>900 ± 30</td><td>-2.33 ± 0.09</td><td>-0.60 ± 0.07</td></tr></tbody></table>

**Note.** The first three rows show the results of the fitting to the averaged peak velocities, and the last row shows the result of the fitting to the peak velocities measured on the southwest side of TMC-1A.

emcee (Foreman-Mackey et al. 2013):

$$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 $r_{\text{break}}$ is the break radius where the power-law function changes. In the double-power-law fitting, ($v_{\text{break}}$, $r_{\text{break}}$, $p_{\text{in}}$, $p_{\text{out}}$) were set as free parameters and their parameter spaces

were searched. The best-fit parameters and their errors are obtained as the mean and standard deviation of the posterior probability distribution, respectively. The results of the fitting are summarized in Table 7.

For TMC-1A, it is appropriate to assume an azimuthally symmetric envelope and use the averaged peak velocity, as previous observations have suggested an asymmetric molecular gas distribution on ~100–500 au scales (Sakai et al. 2016; Harsono et al. 2021). Hence, we performed fittings using the
double-power-law function to the velocities measured on each side of the protostar separately without averaging as well as to the averaged peak velocities. While the fitting for the data points on the southwest side is converged, we could not obtain any solution for the data points on the northeast side with the double-power-law function. These fitting results are also summarized in Table 7.

# 5. Discussion

## 5.1. Nature of the Radial Profiles of the Peak Velocity

### 5.1.1. IRAS 15398−3359

The power-law index within the break radius is $-1.7 \pm 0.3$ for IRAS 15398−3359. A similar measurement was made by Yen et al. (2017), who estimated a power-law index to be $\sim -1.0$ within a radius of $\sim 100$ au in $C^{18}O J=2-1$ and interpreted the power-law index as the rotational motion of the infalling envelope conserving specific angular momentum.

The measured power-law index in the current work is smaller than that derived in the previous work, although the fitting uncertainty is also large. This difference could be due to the contamination by gas motion on large spatial scales, as all velocity components on different spatial scales along the line of sight are included in the current maps. The peak velocity of the spectrum may be more strongly affected by the large-scale gas motion around the centroid velocity of the foreground gas. This stronger contamination by the large-scale gas motion at a particular velocity would change the slope of the radial profile of the peak velocity. Thus, the measured power-law index of $-1.7 \pm 0.3$ could be interpreted as indicating the rotational motion of the envelope with a degree of contamination by the large-scale gas motion.

It is suggested that the large-scale gas is turbulent in IRAS 15398−3359, as discussed in Section 5.2. The contamination by the turbulent, large-scale gas motion could be also the cause of the two small bumps seen in the radial profile of the peak velocity at radii of $\sim 900$ and $\sim 1500$ au in Figure 3(a). The bumps have a width of $\sim 400$ au, and the deviation of the peak velocities around the bumps is $\sim 0.01 \text{ km s}^{-1}$. Our analysis suggests that the turbulent motion in IRAS 15398−3359 is characterized by the relation $\delta v = 10^{-1.69}(\tau/1000 \text{ au})^{0.6}$, where $\delta v$ is the velocity deviation and $\tau$ is the spatial scale (see Section 5.2 for more details). This relation predicts a velocity difference of $\sim 0.01 \text{ km s}^{-1}$ on a scale of 400 au, which matches the deviation of the peak velocities around the bumps.

One might wonder whether the secondary outflow component suggested in a previous work (Okoda et al. 2021) affects our analysis, as it is present on the southeast side in the cut direction of the PV diagram. The secondary outflow component appears at velocities of $4.5–5.0 \text{ km s}^{-1}$ within an offset range of $-6''$ to $-3''$, while the peak velocities measured within the offset range are around $5.18 \text{ km s}^{-1}$. Thus, our measurements mainly trace different velocity components from the secondary outflow. On the other hand, the secondary outflow could skew the overall shape of the spectrum and cause a shift in the peak velocity. A few data points around $\pm 3''$ offsets ($\sim 450$ au) in the PV diagram in Figure 2(a) have blueshifted velocity on both sides of the protostar, while most of the data points have blueshifted and redshifted velocities on the northwest and southeast sides, respectively, as expected for its envelope rotation (Oya et al. 2014; Yen et al. 2017). Those velocity components that are not expected for pure rotational motion on

the southeast side could be influenced by the secondary outflow.

The power-law index of the radial profile outside the break radius, on the other hand, is $0.68 \pm 0.05$ and significantly different from that inside the break radius. This increase of the peak velocity with radius is consistent with previous measurements of velocity structures within dense cores on scales of $\sim 1000–10,000$ au, i.e., $v \sim j/r \propto r^{0.6}-r^{0.8}$ (Goodman et al. 1993; Caselli et al. 2002; Pineda et al. 2019; Gaudel et al. 2020). The origin of the velocity structures outside the break radius is further discussed in Section 5.2.

### 5.1.2. L1527 IRS

The power-law index within the break radius is $-1.38 \pm 0.03$ for L1527 IRS. Aso et al. (2017) have also measured the radial profile of the gas velocity for L1527 IRS within a radius of $\sim 100$ au, giving a power-law index of $-1.22 \pm 0.04$, which roughly agrees with our measurement. Although the indices of the two measurements are consistent with each other, the measured velocities are not consistent between two measurements, as shown in Figure 3(b). The peak velocities measured in the current work are a few times larger than those expected from the extrapolation of the radial profile measured by Aso et al. (2017) while assuming the same systemic velocity. Such a discrepancy could be caused by the contamination from infalling motions due to the relatively lower angular resolution of our measurements, as discussed by Sai et al. (2022). Gaudel et al. (2020) have also measured the line-of-sight velocity as a function of radius at radii of 50–8000 au in L1527 IRS and fitted a double-power-law function to its apparent specific angular momentum $|j_{\text{app}}| = r \times v$. They report a break radius of $1320 \pm 260$ and power-law indices of the radial dependence of the apparent specific angular momentum of $-0.14 \pm 0.03$ and $1.3 \pm 0.2$ inside and outside the break radius (see their Table F.2), corresponding to $p \sim -1.14$ and $0.3$ of $v \propto r^p$, respectively. The break radius and power-law indices they reported are roughly consistent with the values derived in the current work despite the difference in data sets and fitting methods.

The power-law index we measured within the break radius is smaller than $-1$, which is expected for the rotational motion of the envelope conserving specific angular momentum. We also see abrupt changes in the peak velocity seen in the PV diagram (Figure 2), as discussed in Section 4.2. Importantly, the radii where the abrupt changes of the peak velocities occur are different between the north and south sides and are smaller than the break radius, suggesting that entire motions within the break radius cannot be explained by simple rotational motions. A possible explanation for these velocity structures is the contamination from gas motion on larger scales, as discussed for IRAS 15398−3359 in the previous subsection. The actual spatial scale where the velocity structure follows the rotational motion could be, thus, smaller than the estimated break radius, as Gaudel et al. (2020) reported a smaller break radius of 1320 au.

Similarly to IRAS 15398−3359, the power-law index outside the break radius is very different from that inside the break radius for L1527 IRS. The interpretation of the velocity structure outside the break radius is discussed in more detail in Section 5.2.
### 5.1.3. TMC-1A

The radial profile of the peak velocity for TMC-1A is very different from those of the two other protostars described above. The power inside of the break radius, $-4.1 \pm -0.1$ for the averaged peak velocity and $-2.33 \pm -0.09$ on the southwest side of the protostar, is much smaller than that in IRAS 15398–3559 and L1527 IRS. These values also vary greatly from previous measurements of the rotational velocity made by Aso et al. (2015), reporting the power-law index of $-0.86$ within a radius of $\sim 100$ au. Moreover, the peak velocities measured on the southwest and northeast sides of the protostar are mostly inconsistent with each other within the break radius. These facts suggest that the measured peak velocity in TMC-1A in the current work does not seem to trace the rotational velocity of the infalling envelope even within the break radius. One possible reason for this is the asymmetric structure of its infalling envelope. An asymmetric molecular distribution has been reported on $\sim 100–500$ au scales in TMC-1A in previous works (Sakai et al. 2016; Harsono et al. 2021), and it is explained with an infalling streamer-like structure caused by an infalling cloudlet (Hanawa et al. 2022). In such a case where infalling gas is asymmetric, the line-of-sight velocity along the rotational major axis can represent infalling velocity rather than rotational velocity. The projected infalling velocity at a radius of 400 au on the sky plane is a maximum of $\sim 1.9$ km s$^{-1}$ according to the protostellar mass of $0.68 M_\odot$ (Aso et al. 2015). Thus, it would be possible that inflows produce peak velocities of $\sim 0.4–1$ km s$^{-1}$ around a radius of 400 au depending on the projection angle.

The power outside of the break radius, $-0.30 \pm 0.06$ or $-0.60 \pm 0.07$, is also significantly different from those for the other two sources. There seems to be a less clear velocity gradient compared to the other two sources. One possibility is that the foreground gas of TMC-1A has little motion, and its velocity component close to the systemic velocity is dominant in the current maps due to integration along the line of sight. We also note that the peak velocities outside the break radius are small, particularly at radii larger than 2000 au on the southwest side; most of them are less than $0.03$ km s$^{-1}$ and have relatively large uncertainties of $\sim 50\%$ of the measured peak velocities. Thus, the peak velocities measured at large radii could be close to the accuracy limit.

## 5.2. Origin of the Velocity Structure outside the Break Radius

The velocity structures outside the break radii in IRAS 15398–3359 and L1527 IRS are not likely due to simple rotational motion. A possible origin of the incoherent velocity structures on large scales is the nonaxisymmetric collapse of the dense core. A numerical simulation performed by Verliat et al. (2020) demonstrates that such nonaxisymmetric collapse generates local rotational motion, and thus velocity gradients appear at scales of hundreds to thousands of astronomical units. In their simulation where the initial dense core has density perturbations but no motion, the specific angular momentum calculated using projected velocity as $j = r^2\Omega$, where $\Omega$ is the two-dimensional velocity gradient measured with simulated observational maps, is a few times $10^{-4}$ km s$^{-1}$ pc at a radius of 1000 au. This value agrees with the values of $\sim 10^{-4}–10^{-3}$ km s$^{-1}$ pc calculated from the current observations as $j = r \times v_{peak}$. On the other hand, the specific angular momentum calculated in their simulation at a larger radius from

2000 to 8000 au is $(2–6) \times 10^{-4}$ km s$^{-1}$ pc, which is smaller than the $(5–80) \times 10^{-4}$ km s$^{-1}$ pc, calculated from our observations within the same radius range, by 1 order of magnitude at most. This difference may arise because their calculation simulates an isolated self-standing core with no motion that is not embedded in large-scale structures. The angular momentum at larger scales could also depend on the degree of the initial density perturbation in the simulations. As their study focused on the inner envelope scales, further models treating both the large-scale structures and density perturbations are needed to investigate the velocity structures caused by nonaxisymmetric collapse on scales of $\ge 2000$ au in more detail.

Another possibility to explain the kinematic structure outside the break radius is turbulence motion. It is suggested that scale-dependent turbulence in dense cores can produce a velocity gradient on subparsec scales in projected maps (Burkert & Bodenheimer 2000). The turbulence in molecular clouds is known to follow Larson's law, $\Delta v \propto L^{-0.5}$, where L is the cloud size and $\Delta v$ is the velocity deviation over the cloud (Larson 1969; Heyer & Brunt 2004; McKee & Ostriker 2007). The relation of $v_{peak} \propto r^{0.5}-r^{0.6}$ measured in IRAS 15398–3359 and L1527 IRS is similar to the scaling law, as some previous works pointed out the similarity of Larson's law and the relation of $J/M \propto R^{1.6}$ (corresponding to $v \propto r^{0.6}$) measured in dense cores (Goodman et al. 1993; Ohashi et al. 1997b; Tatematsu et al. 2016; Gaudel et al. 2020).

To examine whether the velocity structures on a large scale originate from turbulent motion, spatial correlations of the velocity deviation are examined with the large-scale moment 1 maps (i.e., the centroid velocity) through the second-order structure function (SF) and the autocorrelation function (ACF). The slope of the spatial correlation of the velocity deviation, $\gamma_{2D}$, where $\delta v \propto \tau^{\gamma_{2D}}$, is compared with that observed in molecular clouds and predictions of turbulence models. Note that the slope measured from the centroid velocity is denoted as $\gamma_{2D}$ since it does not necessarily equal the slope $\gamma_{3D}$, which is measured in three dimensions, as discussed later. The large-scale maps include velocity components of infalling envelopes and outflows, as was reviewed in Section 3. To exclude these components, one-dimensional slices of the large-scale moment 1 maps along the disk major axis are used with masks excluding the velocity components inside the break radii. The one-dimensional slices of the moment 1 maps are sampled with half of the beam size along the cut direction. The effect of masking on the slope $\gamma_{2D}$ is examined using model turbulent velocity fields in Appendix B, demonstrating that masking does not affect the value of the slope for the turbulent velocity field significantly.

The discrete, one-dimensional SF and ACF can be defined, respectively, as follows (Miesch & Bally 1994):

$$S(\tau) = N(\tau)^{-1} \sum_x [v_c(x) - v_c(x + \tau)]^2 \quad (5)$$

and

$$C(\tau) = \sigma_{v_c}^{-2} N(\tau)^{-1} \sum_x [v_c(x) - \mu][v_c(x + \tau) - \mu], \quad (6)$$

where $\mu$ is the mean centroid velocity and defined as

$$\mu \equiv \frac{\sum v_c(x)}{N_{\text{tot}}}, \quad (7)$$
<figure><img src="image_6.png" /><figcaption>Figure 4. ACFs calculated from the one-dimensional slices of the large-scale moment 1 maps.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5. Spatial correlations of the velocity deviation around the three protostars. Red solid lines indicate the best-fit power-law relations. Vertical dashed lines show $\tau_0$, within which the fitting was performed.</figcaption></figure>

and $\sigma_{v_c}^2$ is the variance of the centroid velocity and defined as

$$ \sigma_{v_c}^2 = \frac{\sum[v_c(x) - \mu]^2}{N_{tot}}. \qquad (8) $$

Here, $S(\tau)$ is the SF, $C(\tau)$ is the ACF, $x$ is the coordinate in the one-dimensional slice of the moment 1 map, $v_c$ is the centroid velocity, $\tau$ is the spatial separation between two pixels and $N$ is the number of pixels that used to calculate correlations. The spatial correlation of the velocity deviation is obtained as $\delta v(\tau) = \sqrt{S(\tau)}$. ACFs provide information on the effective, largest scale of turbulence $\tau_0$, where $C(\tau_0) = 0$ (Miesch & Bally 1994; Miville-Deschenes et al. 1995). Because the velocity deviation is decorrelated with the spatial scale at scales larger than $\tau_0$, the spatial correlation of the velocity deviation is examined within $\tau_0$, following previous studies (e.g., Miesch & Bally 1994). Uncertainties of ACFs and SFs are estimated through error propagation.

Figure 4 presents ACFs calculated from one-dimensional slices. The ACFs of IRAS 15398–3359 and L1527 IRS almost monotonically decrease until they reach zero, providing $\tau_0$ of ~4000 au and 4400 au, respectively. Calculated $\delta v(\tau)$ for the two sources are presented in Figures 5(a) and (b). The spatial scale and velocity deviation exhibit a strong positive correlation within $\tau_0$. A linear

function is fitted to $\delta v(\tau)$ in the log-log space through $\chi^2$ fitting. The best-fit power-law relation is

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

for IRAS 15398–3359, and

$$ \delta v = 10^{-1.59 \pm 0.04} \left( \frac{\tau}{1000 \text{ au}} \right)^{0.65 \pm 0.08} \quad (\text{km s}^{-1}) \qquad (10) $$

for L1527 IRS. We also performed the same analysis using one-dimensional cuts with position angles that deviate from the position angle of the major axis within ±10° and found that the uncertainty of the slope associated with the cut direction is ~0.1 within the range of the position angle.

We also calculated a spatial correlation of the velocity deviation for TMC-1A, although the break radius of TMC-1A seems not to represent the transition between the infalling envelope with a constant specific angular momentum and the outside. The ACF for TMC-1A also monotonically decreases with the spatial scale, providing $\tau_0$ of 3600 au, as shown in Figure 4(c). The calculated $\delta v(\tau)$ for TMC-1A has a steeper profile compared to the other two sources, as seen in Figure 5(c). As such, the best-fit power-law function is $\delta v = 10^{-1.41 \pm 0.03}(\tau/1000 \text{ au})^{0.93 \pm 0.08}$. All fitting results are summarized in Table 8.
<table><caption>**Table 8**
Results of the Fitting to the Spatial Correlation of the Velocity Deviation</caption><thead><tr><th>Source</th><th>r<sub>mask</sub> (au)</th><th>q<sup>a</sup></th><th>γ<sub>2D</sub><sup>a</sup></th></tr></thead><tbody><tr><td rowspan="3">IRAS 15398–3359</td><td>1200</td><td>-1.69 ± 0.08</td><td>0.6 ± 0.2</td></tr><tr><td>1700</td><td>-1.69 ± 0.08</td><td>0.6 ± 0.2</td></tr><tr><td>2400</td><td>-1.68 ± 0.08</td><td>0.5 ± 0.3</td></tr><tr><td rowspan="2">L1527 IRS</td><td>1700</td><td>-1.59 ± 0.04</td><td>0.65 ± 0.08</td></tr><tr><td>2400</td><td>-1.58 ± 0.04</td><td>0.6 ± 0.1</td></tr><tr><td rowspan="3">TMC-1A</td><td>800</td><td>-1.41 ± 0.03</td><td>0.93 ± 0.08</td></tr><tr><td>1700</td><td>-1.52 ± 0.04</td><td>0.83 ± 0.09</td></tr><tr><td>2400</td><td>-1.50 ± 0.04</td><td>0.74 ± 0.09</td></tr></tbody></table>

**Note.**

$$a \delta v(\tau) = 10^q(\tau/1000 \text{ au})^{\gamma_{2D}}.$$

Spatial correlations of the velocity deviation (or SFs) have been studied in molecular clouds using centroid velocity maps. An earlier work by Miesch & Bally (1994) reported a wide range of $\gamma_{2D}$ from 0.23 to 0.79, with a mean and standard deviation of 0.50 and 0.16, respectively, at spatial scales of ~0.1–1 pc in 14 molecular clouds. Ossenkopf & Mac Low (2002) applied the $\Delta$-variance method, which is similar to the SF but calculates the average variance on a certain spatial scale with the filtered map, to the centroid maps of the Polaris Flare cloud, suggesting that $\gamma_{2D} \sim 0.32–0.68$ within a scale of ~1500 au to 0.75 pc. The slope $\gamma_{2D}$ derived for IRAS 15398–3359 and L1527 IRS in the current work is similar to the larger values derived in molecular clouds mentioned above, although that of TMC-1A is larger than the values reported in the previous works.

In order to compare the slopes derived using the centroid velocity with theoretical predictions, the projection effect has to be considered. The $\gamma_{2D}$ measured with centroid velocity maps is expected to be larger than the $\gamma_{3D}$ measured in three dimensions, representing the true nature of turbulence, since the velocity fluctuation along the line of sight averages out in centroid velocity maps (O'dell & Castaneda 1987; Brunt et al. 2003; Miville-Deschênes et al. 2003). Brunt & Mac Low (2004) studied this effect using HD and MHD simulations in detail, taking into account density inhomogeneity, and derived an empirical equation for the conversion from $\gamma_{2D}$ to $\gamma_{3D}$:

$$\gamma_{2D} = \gamma_{3D} + 0.5 + \delta\kappa/2, \qquad (11)$$

where $\delta\kappa$ is an empirically derived constant. They suggested that, when turbulent energy is being injected (driven turbulence), $\delta\kappa$ ranges from -1.5 to -0.5 and statistically approaches ~-1 (i.e., $\gamma_{2D} \sim \gamma_{3D}$). This is because density fluctuation causes additional small-scale fluctuation in the centroid velocity through intensity weighting, and the projection-smoothing effect is canceled out. The coefficient $\delta\kappa$, on the other hand, decreases to zero when the injection of turbulent energy is turned off and turbulence is decaying (decaying turbulence). We, therefore, estimate the possible $\gamma_{2D}$ values for both cases of driven turbulence ($\delta\kappa \sim -1.5$ to -0.5) and decaying turbulence ($\delta\kappa = 0$) with Equation (11).

Widely accepted theoretical models of isotropic turbulence —Kolmogorov turbulence and Burgers turbulence—expect $\gamma_{3D} = 1/3$ for incompressible fluids (Kolmogorov 1941) and $\gamma_{3D} = 1/2$ for compressible supersonic turbulence (Burgers 1974). In these models, $\gamma_{2D}$ is expected to be ~0.831.0 for the case of decaying turbulence, which is much larger than

those measured for IRAS 15398–3359 and L1527 IRS. On the other hand, the possible ranges of $\gamma_{2D}$ for the case of driven turbulence are ~0.08–0.75. The slopes derived in IRAS 15398–3359 and L1527 IRS are consistent with this range. A simulation work suggests that outflows associated with low-mass protostars can drive and maintain turbulence in an isolated core over an entire lifetime of Class 0 protostars of ~10⁵ yr (Offner & Arce 2014; Offner & Chaban 2017). Hence, it could be possible that turbulence is being driven around these protostars.

TMC-1A shows a spatial correlation of the velocity deviation with a steeper slope than the other two protostars, i.e., the velocity field around TMC-1A has a relatively larger gradient across a large spatial length compared to perturbations on small spatial scales. However, this could be due to systemic motions around the protostar, such as rotational and infalling motions, as the measured break radius would not correspond to the transition between the two regimes seen in IRAS 15398–3359 and L1527 IRS. Indeed, as seen in Figure 3(c), the radial profile of the peak velocity shows a large velocity difference between the inner 400–700 au radius and the outer thousands of astronomical units radius. Gas motions with a high velocity around the protostar would cause a relatively large velocity gradient and thus make the slope $\gamma_{2D}$ steeper. On the other hand, the slope $\gamma_{2D}$ is not expected to change significantly depending on the mask size for the turbulent velocity field. Thus, we calculated SFs with larger mask sizes, $r_{\text{mask}} = 1700$ au and 2400 au, to test whether the steeper slope found in TMC-1A is due to systemic motions around the protostar. The fitting results are summarized in Table 8. The obtained $\gamma_{2D}$ with the largest mask size of 2400 au is much smaller than that derived with a mask of 800 au. We also measured $\gamma_{2D}$ for IRAS 15398–3359 and L1527 IRS with the largest mask size of 2400 au and obtained 0.5 and 0.6, respectively, which are comparable to the values derived with the smaller masks within fitting errors. These results, thus, suggest that the steep slope of 0.93 obtained with a mask size of 800 au for TMC-1A is due to systemic motions around the protostar that have relatively high velocity. The obtained $\gamma_{2D}$ for TMC-1A with the mask size of 2400 au is consistent with the values measured for the other two sources and values reported for molecular clouds in previous works. However, the velocity structure could still originate from envelope components, as mentioned above, and it is difficult to conclude what motion is dominant for the velocity structure on ~1000–10,000 au scales around TMC-1A.

While the derived slopes for IRAS 15398–3359 and L1527 IRS are consistent with a turbulent velocity field, we note that these results do not necessarily rule out the presence of any rotational motion on these scales. It would be possible that the velocity field of a hybrid of rigid-body rotation and turbulence exhibits a spatial correlation of the velocity deviation with a slope similar to the observed slopes, as presented in Figure B1 in Appendix B. Hence, in order to assess the degree of possible rotational motions, a linear function is fitted to the one-dimensional velocity profile and SFs are calculated after subtracting the linear velocity gradients. Linear velocity gradients provide reasonable estimates of the projected rotational velocity, even in cases where the velocity field consists of rotational and turbulent motions (Stewart & Federrath 2022). The results of the fitting to the one-dimensional velocity profiles and spatial correlations of the
<table><caption>**Table 9**
Results of the Fitting to the One-dimensional Cuts of the Moment 1 Maps and the Spatial Correlation of the Velocity Deviation</caption><thead><tr><th>Source</th><th>r<sub>mask</sub><br>(au)</th><th>G<sup>a</sup><br>(km s<sup>-1</sup> pc<sup>-1</sup>)</th><th>v<sub>0</sub><sup>a</sup><br>(km s<sup>-1</sup>)</th><th>q<sup>b</sup></th><th>γ<sub>2D</sub><sup>b</sup></th></tr></thead><tbody><tr><td>IRAS 15398–3359</td><td>1200</td><td>0.5 ± 0.2</td><td>5.203 ± 0.005</td><td>-1.69 ± 0.08</td><td>0.5 ± 0.3</td></tr><tr><td>L1527 IRS</td><td>1700</td><td>1.2 ± 0.1</td><td>5.951 ± 0.004</td><td>-1.59 ± 0.04</td><td>0.5 ± 0.2</td></tr></tbody></table>

**Notes.**

$$^{a}v(x) = Gx + v_0.$$

$$^{b}\delta v(\tau) = 10^q(\tau/1000 \text{ au})^{\gamma_{2D}}.$$

velocity deviation are summarized in Table 9. The slope becomes slightly shallower compared to the slope without subtraction of the linear velocity gradients, as the overall velocity gradient responsible for the velocity deviation on a large scale is suppressed. The slopes are still consistent with a turbulent velocity field for the driven turbulence case. We calculated the square of the velocity deviation caused by the rotation ($\delta v_{\text{rot}}^2$) and turbulence ($\delta v_{\text{turb}}^2$) at a radius of 4000 au from the fitting results. For IRAS 15398–3359 and L1527 IRS, $\delta v_{\text{turb}}^2$ is calculated to be $2 \times 10^{-3} \text{ km}^2 \text{ s}^{-2}$ and $3 \times 10^{-3} \text{ km}^2 \text{ s}^{-2}$ and $\sim 5–17$ times larger than the $\delta v_{\text{rot}}^2$ of $10^{-4} \text{ km}^2 \text{ s}^{-2}$ and $5 \times 10^{-4} \text{ km}^2 \text{ s}^{-2}$, respectively. Although the rotational velocity is not corrected by the inclination angle, it affects $\delta v_{\text{rot}}^2$ by only 0.8%–10%, assuming the core rotational axis is the same as those of the outflow and disk, as the inclination angles of the outflows and disks in IRAS 15398–3359 and L1527 IRS are estimated to be $\sim 70^\circ$ and $\sim 85^\circ$, respectively ($0^\circ$ for a pole-on configuration Tobin et al. 2008; Oya et al. 2014; Aso et al. 2017). Thus, turbulent motions would be more dominant than rotational motions, even when the dense cores also possess a rotational motion.

In two sources, IRAS 15398–3359 and L1527 IRS, the direction of the velocity gradient changes greatly depending on the spatial scale, as seen in Section 4.1. This trend is also consistent with turbulent velocity fields, where the directions of the specific angular momentum axis are expected to vary as a function of radius (Joos et al. 2013; Matsumoto et al. 2017), as is discussed in the next subsection in more detail. A recent simulation work by Misugi et al. (2019) suggests that velocity fluctuation inside the filament, described with one-dimensional Kolmogorov power spectra, explains the $J/M-R$ relation measured in dense cores (e.g., Goodman et al. 1993; Tatematsu et al. 2016). Our results agree with this scenario. Hacar et al. (2013) measured the turbulent velocity of filaments to be typically sonic on a 0.5 pc scale. Assuming the Kolmogorov scaling law, it is expected that turbulent velocity is about $0.06 \text{ km s}^{-1}$ at 4000 au, which is comparable to the values observed in IRAS 15398–3359 and L1527 IRS.

We note that, in another possible explanation of the velocity structures outside the break radius, Arroyo-Chávez & Vázquez-Semadeni (2022) recently demonstrated with numerical simulations that the $J/M-R$ relation observed in dense cores can be explained as a consequence of gravitational contraction with angular momentum loss via turbulent viscous. This scenario may be also applied to the relation of $j \propto r^{-1.6}$ found in the radial profile of the specific angular momentum measured in individual protostellar cores. However, as they calculated the angular momentum and the dense core radius in three dimensions, and placed focus on the spatial scale of $\sim 0.1–1$ pc, it is difficult to make a direct comparison between their simulations and our observed maps on $\sim 100–10,000$ au

scales. Simulated projected maps on $\sim 1000–10,000$ au scales will allow us to investigate how the velocity structures in individual dense cores look like in this scenario as well as whether they are consistent with our observational results.

## 5.3. Comparison with Other Sources

### 5.3.1. Break Radius

The two Class 0 sources in our sample, L1527 IRS and IRAS 15398–3359, show a break at a radius of $\sim 1200$ and $1700$ au in the radial profile of the peak velocity, respectively. These breaks may correspond to the transition between the $j$-constant and $j$-increase regimes discussed in previous works (Gaudel et al. 2020; Sai et al. 2022), although the velocity structures outside the break radius around the two sources are not likely explained as rotational motion. Gaudel et al. (2020) reported a break radius of 1600 au, which is similar to that of IRAS 15398–3359 and L1527 IRS, based on an averaged radial profile of the specific angular momentum of 12 Class 0 sources. Sai et al. (2022) reported that the Class I protostar L1489 IRS exhibits a break at a radius of $\sim 2900$ au, which is about two times larger than that in L1527 IRS and IRAS 15398–3359. Pineda et al. (2019) measured radial profiles of the specific angular momentum distributions within a radius of 800–10,000 au around two Class 0 sources (HH 211 and IRAS 03282+3035) and one first hydrostatic core candidate (L1451-mm), obtaining $j \propto r^{1.8}$ without a break. Lee (2010) found that the Class I protostar HH 111 exhibits rotation motion with a constant specific angular momentum within a radius of 2000–7000 au, suggesting that the transition between the two regimes occurs at a radius $\gtrsim 7000$ au.

These comparisons of the break radius, including lower and upper limits and summarized in Table 10, suggest that the break radius between the two regimes is typically less than 2000 au for Class 0 sources. It is also found that the break radii of Class I sources are larger than those of Class 0 sources, although there are only two samples of Class I sources. This possible trend suggests that the extent of the infalling envelope with a constant specific angular momentum increases as the protostar evolves. This finding is consistent with the results of an analytic core-collapse model calculated by Takahashi et al. (2016). They reported that the infalling gas element is extended in the radial direction during core collapse when the angular momentum is conserved, as the infalling velocity is faster at inner radii, which results in the flat radial distribution of the specific angular momentum (see also Yen et al. 2011, 2017). These calculations predict that the extent of the flat radial distribution of the specific angular momentum increases with evolution, as larger specific angular momentum is brought in from larger radii at later evolutionary stages. It should be noted, however, that these models assume a very simple initial
<table><caption>**Table 10**
Break Radii Estimated in Protostellar and Prestellar Sources</caption><thead><tr><th>Evolutionary Stage</th><th>Source/Number of Sources</th><th>r<sub>break</sub> (au)</th><th>References</th></tr></thead><tbody><tr><td>FHSC<sup>a</sup> candidate/Class 0</td><td>3 sources</td><td>&lt;1000</td><td>Pineda et al. (2019)</td></tr><tr><td>Class 0</td><td>12 sources</td><td>1600 ± 300</td><td>Gaudel et al. (2020)</td></tr><tr><td>Class 0</td><td>IRAS 15398-3359</td><td>1180 ± 90</td><td>This work</td></tr><tr><td>Class 0</td><td>L1527 IRS</td><td>1700 ± 50</td><td>This work</td></tr><tr><td>Class I</td><td>L1489 IRS</td><td>2900 ± 200</td><td>Sai et al. (2022)</td></tr><tr><td>Class I</td><td>HH 111</td><td>&gt;7000</td><td>Lee (2010)</td></tr></tbody></table>

**Note.**

<sup>a</sup> First hydrostatic core.

condition, such as a dense core rotating like a rigid body without any turbulence or magnetic field. Numerical simulations with more realistic conditions are required to further investigate whether the above picture is valid in turbulent dense cores.

### 5.3.2. Velocity Structure outside the Break Radius

The velocity structures outside the break radius around IRAS 15398–3359 and L1527 IRS are consistent with turbulent motion, as discussed in Section 5.2. Gaudel et al. (2020) also revealed that 12 Class 0 protostars exhibit large dispersions of the directions of velocity gradients over radii of 100–5000 au and proposed the hypothesis of turbulent motions as the best interpretation of the velocity structures at scales of >1600 au. On the other hand, the velocity structure outside the break radius around L1489 IRS is interpreted as rotational motion by Sai et al. (2022), where the magnitude and direction of the velocity gradient of L1489 IRS were measured with the same method used in Section 4.1. For data points of L1489 IRS at $r_{fit} \le 10''$, we performed the same fitting as used in Section 4.1 to the small-scale map of L1489 IRS, which is presented in Appendix C.

To compare the nature of velocity structures among different sources, the directions of the velocity gradients are measured as a function of the spatial scale, as shown in Figure 6, based on the measurements presented in Section 4.1. We exclude TMC-1A from comparison here, as we could not measure the break radius, which likely corresponds to a transition between the *j*-constant and *j*-increase regimes in the source. For L1489 IRS, data points at $r_{fit}$ of 20"–60" are compiled from measurements in Sai et al. (2022), where the magnitude and direction of the velocity gradient of L1489 IRS were measured with the same method used in Section 4.1. For data points of L1489 IRS at $r_{fit} \le 10''$, we performed the same fitting as used in Section 4.1 to the small-scale map of L1489 IRS, which is presented in Appendix C.

As seen in Figure 6, the directions of the velocity gradients of IRAS 15398–3359 and L1527 IRS significantly vary by about 90°–100° from the smallest to largest scales. The dispersion of the directions is 34° for IRAS 15398–3359 and 45° for L1527 IRS. On the contrary, the dispersion of the directions of the velocity gradients of L1489 IRS is only $\sim \pm 5°$ across its break radius of 2900 au ($\sim 21''$).

One might wonder whether the large dispersion of the directions of the velocity gradients in IRAS 15398–3359 and L1527 IRS is because outflow components are dominant at particular spatial scales. Indeed, the position angles of the velocity gradient measured in IRAS 15398–3359 with the small-scale map are $\sim 70°$ and similar to that of the outflow (Yen et al. 2017). The position angles of the velocity gradient in the large-scale map of L1527 IRS are $\sim 110°$ and are also close to the position angle of the outflow ($\sim 90°$; Tobin et al. 2008).

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Directions of velocity gradients measured within aperture areas with a radius of $r_{fit}$. In the top panel, P.A$_{cent}$ is the P.A. of the velocity gradient measured within $r_{fit} = 5''$. Vertical dashed lines denote the break radius measured in each source.</figcaption></figure>

To measure the orientations of the velocity gradients in the two sources minimizing the contribution of the outflows, we performed the same fitting to maps where the area within $\pm 60°$ from the outflow axis is masked. The position angle of the outflow is assumed to be 50° for IRAS 15398–3359 and 90° for L1527 IRS, based on previous observations (Tobin et al. 2008; Yen et al. 2017). The fitting results with the masks are presented in Figure 7. The two sources still show a larger dispersion of the directions of the velocity gradients than L1489 IRS, even when the emission associated with the outflows is not considered in the analysis. The dispersion measured in IRAS 15398–3359 with the outflow masking is 19°, which is smaller than when measured without the outflow masking but still larger than in L1489 IRS. The direction of the velocity gradient varies by 50° at maximum. An inversion of the direction by 180° is seen at a 20" scale in L1527 IRS, as discussed in Section 5.1. This is because an overall velocity gradient caused by the outflow is less significant. This leads to a larger dispersion of the orientation of 110° compared to when it is measured without the outflow masking.

Gaudel et al. (2020) presented similar measurements in their Figure 19, showing that directions of velocity gradients vary more than 50° from 100 to 5000 au in all 12 Class 0 sources.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Same as Figure 6 but with measurements where outflows are masked for IRAS 15398–3359 and L1527 IRS.</figcaption></figure>

Simulations of turbulent dense cores suggest that the directions of the specific angular momentum axis tend to vary as a function of radius (Joos et al. 2013; Matsumoto et al. 2017), as would be the case for IRAS 15398–3359 and L1527 IRS. These imply that the small dispersion of the direction of the velocity gradient outside the break radius around L1489 IRS, which is interpreted as rotational motion, is different from those around the other sources.

A possible explanation for the difference in the velocity structures outside the break radius between L1489 IRS and the other sources is the difference in the initial condition. Two out of three protostellar or prestellar sources, where Pineda et al. (2019) have measured $j \propto r^{1.8}$, exhibit clear velocity gradients in directions perpendicular to their outflows. The velocity structure on several thousands of astronomical unit scales around the isolated Class 0 protostar B335 is well explained by the rigid-body rotation of the dense core (Saito et al. 1999; Kurono et al. 2013). These observations suggest that some Class 0 sources have an ordered velocity structure with less variation in the direction of local velocity gradients on several thousands of astronomical unit scales, which would be typically outside the infalling envelope with a constant specific angular momentum in the Class 0 phase. On the other hand, all 12 sources studied by Gaudel et al. (2020) show a dispersion of the directions of velocity gradients at scales >1600 au. Tobin et al. (2011) investigated the gas kinematics around 18 protostars on 1000 au to 0.1 pc scales, reporting that velocity gradients over ~0.1 pc scales are often not normal to their outflows. These works with a large sample number suggest that the ordered velocity structure described by rotational motion is rare on ~1000–10,000 au scales in the Class 0 phase. Furthermore, a warped disk structure is suggested around L1489 IRS, which is expected to form with a turbulent initial dense core (Sai et al. 2020). The other possibility explaining the difference in velocity structures is the evolutionary effect. Rotational motion could be more dominant even outside the break radius at later evolutionary stages since energy injection by outflows is expected to decrease with time at the Class I

phase by orders of magnitude (Machida & Hosokawa 2013) and turbulence would decay. If the difference in velocity structures outside the break radius is due to the evolutionary effect, Class I sources will tend to show more coherent velocity structures than Class 0 sources. Currently, L1489 IRS is the only example of a Class I protostar where the break radius is measured. A comparison of the velocity structures between Class 0 and I sources with more samples is needed to investigate the evolution of the velocity structures in more detail.

## 6. Summary and Conclusion

We have observed the three protostars, IRAS 15398–3359, L1527 IRS, and TMC-1A, with the ALMA 12 m array, ACA 7 m array, APEX, and the IRAM-30m telescope in the C$^{18}$O $J=2-1$ emission. The kinematic structures on ~100–10,000 au scales have been investigated with two different types of combined maps: small-scale maps covering regions within a radius of ~800 au around the protostars with spatial resolutions of ~110 au, and large-scale maps covering wider areas within a radius of ~9000 au with spatial resolutions of ~1000 au. Our main results are summarized as follows:

1. In the small-scale maps, two out of the three sources, L1527 IRS and TMC-1A, exhibit clear velocity gradients that likely originate from the rotational motions of the disks and envelopes. IRAS 15398–3359 does not show a clear velocity gradient in the direction perpendicular to the primary outflow but shows a velocity gradient along the primary outflow direction. In large-scale maps, on the other hand, velocity gradients are less clear, and velocity structures are more perturbed in all the sources.

2. To investigate the radius at which the dominant motion changes from the rotational motion of the envelope to others, the peak velocity $v_{peak} = |v_{LSR} - v_{sys}|$ was measured as a function of radius from ~400 to 9000 au using PV diagrams cut along the disk major axis. IRAS 15398–3359 and L1527 IRS exhibit a break at a radius of ~1200 and 1700 au, respectively, in the radial profiles of the peak velocity. Fitting a double-power-law function results in the power-law index of $-1.7 \pm 0.3$ and $-1.38 \pm 0.03$ inside the break radius, which is smaller than the power-law index of -1 expected for the infalling envelope under conservation of specific angular momentum. These slopes would be interpreted as indicating rotational motions of the envelope with contamination by gas motion on larger spatial scales along the line of sight. The power-law indices outside the break radius are, on the other hand, ~0.46–0.68, which are similar to the slope of the $J/M-R$ relation of dense cores ($v \sim J/M/R \propto R^{0.6}$). TMC-1A shows a break in the radial profile of its averaged peak velocity at a radius of ~800 au, and the power-law indices of ~-4.1 and ~-0.30 inside and outside the break radius, respectively. These slopes are consistent with neither the slopes expected for the infalling envelope conserving specific angular momentum nor the $J/M-R$ relation of dense cores.

3. In order to probe the origin of the velocity structures outside the break radius, spatial correlations of the velocity deviation were investigated through the second-order SF of the centroid velocity. The spatial correlations of the velocity deviation, $\delta v \propto \tau^{-0.6}$, were
found in IRAS 15398-3359 and L1527 IRS. The values of the derived slopes are consistent with predictions of turbulence models within an uncertainty associated with the projection effect, assuming that turbulence is being driven around the protostars. These results may suggest that the velocity structures outside the break radius around the two sources originate from turbulence. TMC-1A exhibits a steeper slope, 0.93, which is likely due to the contamination by systemic motions around the protostar.

4. The break radius estimated for IRAS 15398-3359 and L1527 IRS is similar to that derived in previous works. A comparison of the break radii between Class 0 and I sources suggests a possible evolutionary trend that the extent of the infalling envelope with a constant specific angular momentum increases with time, although the number of Class I samples is limited. The velocity structure around the Class I source L1489 IRS appears more coherent across the break radius compared to those around other Class 0 sources. This difference would be due to the difference in the initial condition or turbulence decay at later evolutionary stages. Further observations of more samples at different evolutionary stages are needed to distinguish these two possibilities.

This paper used the following ALMA data: ADS/JAO. ALMA \#2013.1.00879.S, 2013.1.01086.S, and 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 Nos. 044-14 and 129-19 with the IRAM 30-m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). This publication is based on data acquired with the Atacama Pathfinder Experiment (APEX). APEX is a collaboration between the Max-Planck-Institut fur Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory. N.O. is supported by the National Science and Technology Council (NSTC) in Taiwan through the grant NSTC 109-2112-M-001-051 and 110-2112-M-001-031. H.-W.Y. acknowledges support from the Ministry of Science and Technology (MOST) in Taiwan through the grant MOST 110-2628-M-001-003-MY3 and from the Academia Sinica Career Development Award (AS-CDA-111-M03). A.M. is supported by the European Research Council (ERC Starting Grant MagneticYSOs with grant agreement No. 679937).

Facilities: ALMA, IRAM-30 m, APEX.

Software: CASA (McMullin et al. 2007), GILDAS (http://www.iram.fr/IRAMFR/GILDAS), MIRIAD (Sault et al. 1995), Numpy (Oliphant 2006; van der Walt et al. 2011), Scipy (Virtanen et al. 2020), Astropy (Astropy Collaboration et al. 2013, 2018), Matplotlib (Hunter 2007), emcee (Foreman-Mackey et al. 2013).

## Appendix A
### Velocity Channel Maps

Velocity channel maps of the C$^{18}$O J = 2-1 emission of IRAS 15398-3359, L1527 IRS, and TMC-1A are presented in Figures A1-A6.

<figure><img src="image_10.png" /><figcaption>**Figure A1.** Small-scale channel maps of the C$^{18}$O J = 2-1 emission of IRAS 15398-3359. Contour levels are 3, 6, 12, 24, ... × σ, where 1σ = 5.2 mJy beam$^{-1}$. The labels in the top-left corner indicate the local standard of rest (LSR) velocity of each channel in km s$^{-1}$. Crosses at the center and a filled ellipse in the bottom-left corner denote the protostellar position and the beam size, respectively. The systemic velocity of IRAS 15398-3359 is 5.18 km s$^{-1}$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure A2.** Same as Figure A1 but for L1527 IRS with $1\sigma$ = 12 mJy beam$^{-1}$. Maps are shown in steps of two channels. The systemic velocity of L1527 IRS is 5.8 km s$^{-1}$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure A3.** Same as Figure A1 but for TMC-1A with $1\sigma$ = 12 mJy beam$^{-1}$. Maps are shown in steps of two channels. The systemic velocity of TMC-1A is 6.4 km s$^{-1}$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure A4.** Large-scale channel maps of the C$^{18}$O 2–1 emission of IRAS 15398–3359. Contour levels are 3, 6, 12, 24, ... × σ, where 1σ = 0.16 Jy beam$^{-1}$. Dashed boxes show the size of the small-scale map (12" in radius). The labels in the top-left corner indicate the LSR velocity of each channel in km s$^{-1}$. Crosses at the center and a filled ellipse in the bottom-left corner denote the protostellar position and the beam size, respectively. The systemic velocity of IRAS 15398–3359 is 5.18 km s$^{-1}$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure A5.** Same as Figure A4 but for L1527 IRS with $1\sigma = 0.29$ Jy beam$^{-1}$. Maps are shown in steps of two channels. The systemic velocity of L1527 IRS is 5.8 km s$^{-1}$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure A6.** Same as Figure A4 but for TMC-1A with $1\sigma = 0.30$ Jy beam$^{-1}$. Maps are shown in steps of two channels. The systemic velocity of TMC-1A is 6.4 km s$^{-1}$.</figcaption></figure>

## Appendix B
One-dimensional Model of the Turbulent Velocity Field

The slope of the spatial correlation of the velocity deviation obtained from one-dimensional velocity fields is examined using models of the one-dimensional turbulent velocity field in order to assess the robustness of the discussion in Section 5.2.

The model turbulent velocity field is computed based on the method described by Dubinski et al. (1995). The Kolmogorov-

like power spectrum of the turbulent velocity field is given by

$$P_v(k) = \langle |v(k)|^2 \rangle = C(k^2 + k_{\min}^2)^{p_k/2}, \quad (B1)$$

where k is the spatial frequency, $k_{\min}$ is the spatial frequency of the largest scale of turbulence and $p_k$ is the slope of the power spectrum. Here, $\langle \cdot \rangle$ denotes the ensemble average. The given power spectrum results in the spatial correlation of the velocity deviation, $\delta v \propto \tau^{(p_k-n)/2}$, where n is the dimension. The
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**Figure B1.** (Left) Examples of one-dimensional velocity fields of pure rotational motion, pure turbulent motion, and hybrids of rotational and turbulent motions. (Right) Spatial correlations of the velocity deviation for the one-dimensional velocity fields shown in the left panel. Solid lines show the structure functions calculated from the velocity fields, the vertical dashed-dotted lines show $\tau_0$ estimated from their ACFs, and the dashed lines show the best-fit power-law relations. Labels at the top-left corner indicate the best-fit slopes.</figcaption></figure>

velocity field in $k$ space is calculated as follows (Myers & Gammie 1999):

$$v(k) = v_k \exp(ikx)\exp(i\phi_k), \quad (B2)$$

where $v_k$ is the amplitude determined by the power spectrum and $\phi_k$ is a random phase uniformly distributed from 0 to $2\pi$. The velocity field in real space is obtained as the inverse Fourier transform of $v(k)$.

In our calculations, we adopted $p_k = 2$, corresponding to Larson's law $\delta v \propto \tau^{0.5}$. To simulate the current observations, we have set the domain size and the grid number to $L = 17,000$ au and 51, respectively, which are comparable to those of the sampled one-dimensional slice of the moment 1 maps. The largest scale of turbulence, $\lambda_{max} = 1/k_{min}$, is taken to be $2L$. The constant $C$ in the power spectrum is chosen so that the velocity deviation over the entire length is ~0.03 km s$^{-1}$, which is comparable to that of the observations. An example of the turbulent velocity field and the correlation between the spatial scale and velocity deviation is shown in Figure B1.

We also present the one-dimensional velocity field and the spatial correlation of the velocity deviation for the cases of rigid-body rotation and a hybrid of turbulence and rigid-body rotation in Figure B1. The model calculations demonstrate that the slope of the spatial correlation of the velocity deviation for the rigid-body rotation is 1 and much steeper than that for the turbulent velocity field. A hybrid velocity field of rigid-body rotation and turbulence, on the other hand, can show a slope similar to that of the turbulent velocity field.

To assess the effect of masking, we applied the same analysis as that performed in Section 5.2 to the model velocity field with and without a mask. The uncertainty of the measured slope is evaluated with the MC method by performing the analysis with 1000 model velocity fields, which have random phases and obtained as the 68% highest density interval of the posterior probability distribution. Figure B2 shows the relation between the mask size and obtained slopes of the spatial correlation of the velocity deviation. A larger mask size increases the uncertainty of the slope, but the slope is always around the expected value of 0.5, which indicates that the

<figure><img src="image_18.png" /><figcaption>**Figure B2.** Comparison between the mask size and obtained slopes of the spatial correlations of the velocity deviation for the one-dimensional turbulent velocity field.</figcaption></figure>

masking does not significantly change the obtained slope. According to these calculations, the uncertainty of the slope associated with the mask size for IRAS 1539–3359 and L1527 IRS is about 0.1, assuming that their true slope is around 0.5.

## Appendix C
### The Small-scale Map of L1489 IRS

For comparison of velocity gradients, we produced the small-scale map of the Class I protostar L1489 IRS combining the observational data of ALMA 12 m array, ACA 7 m array, and IRAM-30 m, whose observational details are presented in Sai et al. (2020, 2022). The imaging and data combining were performed with the same parameters and method described in
<figure><img src="image_19.png" /><figcaption>Figure C1. Same as Figure 1 but for the small-scale map of L1489 IRS.</figcaption></figure>

<table><caption>Table C1
Result of the Two-dimensional Linear Fitting to the Small-scale Centroid Velocity Map of L1489 IRS</caption><thead><tr><th>r<sub>fit</sub><br>(")</th><th>v<sub>0</sub><sup>a</sup><br>(km s<sup>-1</sup>)</th><th>G<br>(km s<sup>-1</sup> pc<sup>-1</sup>)</th><th>θ<br>(°)</th></tr></thead><tbody><tr><td>5</td><td>7.20</td><td>283.8 ± 0.2</td><td>-120.61 ± 0.04</td></tr><tr><td>10</td><td>7.17</td><td>77.13 ± 0.04</td><td>-122.94 ± 0.03</td></tr></tbody></table>

**Note.**

a Fitting errors of $v_0$ are less than 0.1%.

Section 2.4. The velocity resolution and rms of the resultant map are 0.084 km s⁻¹ and 8.1 mJy beam⁻¹, respectively. The synthesized beam size is 0"94 × 0"81 (−16°). The integrated intensity and centroid velocity maps are presented in Figure C1.

We measured the magnitude and direction of the velocity gradient in the small-scale map of L1489 IRS with the method described in Section 4.1. The fitting result is summarized in Table C1.

## ORCID iDs

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

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

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

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

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

## References

* Arce, H. G., & Sargent, A. I. 2004, ApJ, 612, 342

* Arroyo-Chávez, G., & Vázquez-Semadeni, E. 2022, ApJ, 925, 78

* 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

* Belloche, A., André, P., Despois, D., & Blinder, S. 2002, A&A, 393, 927

* Brunt, C. M., Heyer, M. H., Vázquez-Semadeni, E., & Pichardo, B. 2003, ApJ, 595, 824

* Brunt, C. M., & Mac Low, M.-M. 2004, ApJ, 604, 196

* Burgers, J. M. 1974, The Nonlinear Diffusion Equation (1st ed.; Dordrecht Springer)

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

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

* Chen, C.-Y., & Ostriker, E. C. 2018, ApJ, 865, 34

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

* Chen, H. H.-H., Pineda, J. E., Goodman, A. A., et al. 2019b, ApJ, 877, 93

* Comerón, F. 2008, in The Lupus Clouds, ed. B. Reipurth, Vol. 5 (San Francisco, CA: ASP), 295

* Di Francesco, J., Myers, P. C., Wilner, D. J., Ohashi, N., & Mardones, D. 2001, ApJ, 562, 770

* Dubinski, J., Narayan, R., & Phillips, T. G. 1995, ApJ, 448, 226

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

* Froebrich, D. 2005, ApJS, 156, 169

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

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

* Guilloteau, S., & Dutrey, A. 1998, A&A, 339, 467

* Hacar, A., Tafalla, M., Kauffmann, J., & Kovács, A. 2013, A&A, 554, A55

* Hanawa, T., Sakai, N., & Yamamoto, S. 2022, ApJ, 932, 122

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

* Harsono, D., van der Wiel, M. H. D., Bjerkeli, P., et al. 2021, A&A, 646, A72

* Hayashi, M., Ohashi, N., & Miyama, S. M. 1993, ApJL, 418, L71

* Heyer, M. H., & Brunt, C. M. 2004, ApJL, 615, L45

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

* 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

* Kolmogorov, A. 1941, DoSSR, 30, 301

* Kristensen, L. E., van Dishoeck, E. F., Bergin, E. A., et al. 2012, A&A, 542, A8

* Kurono, Y., Saito, M., Kamazaki, T., Morita, K.-I., & Kawabe, R. 2013, ApJ, 765, 85

* Larson, R. B. 1969, MNRAS, 145, 271

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

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

* Machida, M. N., & Hosokawa, T. 2013, MNRAS, 431, 1719

* 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

* Miesch, M. S., & Bally, J. 1994, ApJ, 429, 645

* Misugi, Y., Inutsuka, S.-i., & Arzoumanian, D. 2019, ApJ, 881, 11

* Miville-Deschenes, M.-A., Joncas, G., & Durand, D. 1995, ApJ, 454, 316

* Miville-Deschênes, M. A., Levrier, F., & Falgarone, E. 2003, ApJ, 593, 831

* Momose, M., Ohashi, N., Kawabe, R., Nakano, T., & Hayashi, M. 1998, ApJ, 504, 314

* Myers, P. C., & Gammie, C. F. 1999, ApJL, 522, L141

* O'dell, C. R., & Castaneda, H. O. 1987, ApJ, 317, 686

* Offner, S. S. R., & Arce, H. G. 2014, ApJ, 784, 61

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

* Ohashi, N., Hayashi, M., Ho, P. T. P., & Momose, M. 1997a, ApJ, 475, 211

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

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

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

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

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

* Ossenkopf, V., & Mac Low, M. M. 2002, A&A, 390, 307

* Oya, Y., Sakai, N., Lefloch, B., et al. 2015, ApJ, 812, 59

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

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

* Punanova, A., Caselli, P., Feng, S., et al. 2018, ApJ, 855, 112

* Sai, J., Ohashi, N., Maury, A. J., et al. 2022, ApJ, 925, 12

* Sai, J., Ohashi, N., Saigo, K., et al. 2020, ApJ, 893, 51
Saito, M., Sunada, K., Kawabe, R., Kitamura, Y., & Hirano, N. 1999, ApJ, 518, 334

Sakai, N., Oya, Y., López-Sepulcre, A., et al. 2016, ApJL, 820, L34

Sakai, N., Oya, Y., Sakai, T., et al. 2014, ApJL, 791, L38

Sault, R. J., Teuben, P. J., & Wright, M. C. H. 1995, in ASP Conf. Ser. 77, Astronomical Data Analysis Software and Systems IV, ed. R. A. Shaw, H. E. Payne, & J. J. E. Hayes (San Francisco, CA: ASP), 433

Stewart, M., & Federrath, C. 2022, MNRAS, 509, 5237

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

Takakuwa, S., Kamazaki, T., Saito, M., & Hirano, N. 2003, ApJ, 584, 818

Tamura, M., Ohashi, N., Hirano, N., Itoh, Y., & Moriarty-Schieven, G. H. 1996, AJ, 112, 2076

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

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

Tobin, J. J., Hartmann, L., Bergin, E., et al. 2012a, ApJ, 748, 16

Tobin, J. J., Hartmann, L., Calvet, N., & D'Alessio, P. 2008, ApJ, 679, 1364

Tobin, J. J., Hartmann, L., Chiang, H.-F., et al. 2011, ApJ, 740, 45

Tobin, J. J., Hartmann, L., Chiang, H.-F., et al. 2012b, Natur, 492, 83

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

Verliat, A., Hennebelle, P., Maury, A. J., & Gaudel, M. 2020, A&A, 635, A130

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

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

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