# Early Planet Formation in Embedded Disks (eDisk) VI:
Kinematic Structures around the Very Low Mass Protostar IRAS 16253-2429

Yusuke Aso Korea Astronomy and Space Science Institute, 776 Daedeok-daero,
Yuseong-gu, Daejeon 34055, Republic of Korea  Woojin Kwon Department of Earth
Science Education, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul
08826, Republic of Korea  SNU Astronomy Research Center, Seoul National
University, 1 Gwanak-ro, Gwanak-gu, Seoul 08826, Republic of Korea  Nagayoshi
Ohashi Academia Sinica Institute of Astronomy and Astrophysics, 11F of
Astronomy-Mathematics Building, AS/NTU, No.1, Sec. 4, Roosevelt Rd, Taipei
10617, Taiwan  Jes K. Jørgensen Niels Bohr Institute, University of
Copenhagen, Øster Voldgade 5–7, DK 1350 Copenhagen K., Denmark  John J. Tobin
National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA
22903 USA  Yuri Aikawa Department of Astronomy, Graduate School of Science,
The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan  Itziar
de Gregorio-Monsalvo European Southern Observatory, Alonso de Cordova 3107,
Casilla 19, Vitacura, Santiago, Chile  Ilseung Han Division of Astronomy and
Space Science, University of Science and Technology, 217 Gajeong-ro, Yuseong-
gu, Daejeon 34113, Republic of Korea  Korea Astronomy and Space Science
Institute, 776 Daedeok-daero, Yuseong-gu, Daejeon 34055, Republic of Korea
Miyu Kido Department of Physics and Astronomy, Graduate School of Science and
Engineering, Kagoshima University, 1-21-35 Korimoto, Kagoshima, Kagoshima
890-0065, Japan  Patrick M. Koch Academia Sinica Institute of Astronomy and
Astrophysics, 11F of Astronomy-Mathematics Building, AS/NTU, No.1, Sec. 4,
Roosevelt Rd, Taipei 10617, Taiwan  Shih-Ping Lai Institute of Astronomy,
National Tsing Hua University, No. 101, Section 2, Kuang-Fu Road, Hsinchu
30013, Taiwan  Center for Informatics and Computation in Astronomy, National
Tsing Hua University, No. 101, Section 2, Kuang-Fu Road, Hsinchu 30013, Taiwan
Department of Physics, National Tsing Hua University, No. 101, Section 2,
Kuang-Fu Road, Hsinchu 30013, Taiwan  Academia Sinica Institute of Astronomy
and Astrophysics, 11F of Astronomy-Mathematics Building, AS/NTU, No.1, Sec. 4,
Roosevelt Rd, Taipei 10617, Taiwan  Chang Won Lee Division of Astronomy and
Space Science, University of Science and Technology, 217 Gajeong-ro, Yuseong-
gu, Daejeon 34113, Republic of Korea  Korea Astronomy and Space Science
Institute, 776 Daedeok-daero, Yuseong-gu, Daejeon 34055, Republic of Korea
Jeong-Eun Lee Department of Physics and Astronomy, Seoul National University,
1 Gwanak-ro, Gwanak-gu, Seoul 08826, Republic of Korea  Zhi-Yun Li University
of Virginia, 530 McCormick Rd., Charlottesville, Virginia 22904, USA  Zhe-Yu
Daniel Lin University of Virginia, 530 McCormick Rd., Charlottesville,
Virginia 22904, USA  Leslie W. Looney Department of Astronomy, University of
Illinois, 1002 West Green St, Urbana, IL 61801, USA  Suchitra Narayanan
Institute for Astronomy, University of Hawai‘i at Mānoa, 2680 Woodlawn Dr.,
Honolulu, HI 96822, USA  Nguyen Thi Phuong Korea Astronomy and Space Science
Institute, 776 Daedeok-daero, Yuseong-gu, Daejeon 34055, Republic of Korea
Department of Astrophysics, Vietnam National Space Center, Vietnam Academy of
Science and Technology, 18 Hoang Quoc Viet, Cau Giay, Hanoi, Vietnam  Jinshi
Sai (Insa Choi) Academia Sinica Institute of Astronomy and Astrophysics, 11F
of Astronomy-Mathematics Building, AS/NTU, No.1, Sec. 4, Roosevelt Rd, Taipei
10617, Taiwan  Kazuya Saigo Department of Physics and Astronomy, Graduate
School of Science and Engineering, Kagoshima University, 1-21-35 Korimoto,
Kagoshima, Kagoshima 890-0065, Japan  Alejandro Santamaría-Miranda European
Southern Observatory, Alonso de Cordova 3107, Casilla 19, Vitacura, Santiago,
Chile  Rajeeb Sharma Niels Bohr Institute, University of Copenhagen, Øster
Voldgade 5–7, DK 1350 Copenhagen K., Denmark  Shigehisa Takakuwa Department of
Physics and Astronomy, Graduate School of Science and Engineering, Kagoshima
University, 1-21-35 Korimoto, Kagoshima, Kagoshima 890-0065, Japan  Academia
Sinica Institute of Astronomy and Astrophysics, 11F of Astronomy-Mathematics
Building, AS/NTU, No.1, Sec. 4, Roosevelt Rd, Taipei 10617, Taiwan  Travis J.
Thieme Institute of Astronomy, National Tsing Hua University, No. 101, Section
2, Kuang-Fu Road, Hsinchu 30013, Taiwan  Center for Informatics and
Computation in Astronomy, National Tsing Hua University, No. 101, Section 2,
Kuang-Fu Road, Hsinchu 30013, Taiwan  Department of Physics, National Tsing
Hua University, No. 101, Section 2, Kuang-Fu Road, Hsinchu 30013, Taiwan
Kengo Tomida Astronomical Institute, Graduate School of Science, Tohoku
University, Sendai 980-8578, Japan  Jonathan P. Williams Institute for
Astronomy, University of Hawai‘i at Mānoa, 2680 Woodlawn Dr., Honolulu, HI
96822, USA  Hsi-Wei Yen Academia Sinica Institute of Astronomy and
Astrophysics, 11F of Astronomy-Mathematics Building, AS/NTU, No.1, Sec. 4,
Roosevelt Rd, Taipei 10617, Taiwan

###### Abstract

Precise estimates of protostellar masses are crucial to characterize the
formation of stars of low masses down to brown-dwarfs (BDs;
$M_{*}<0.08~{}M_{\sun}$). The most accurate estimation of protostellar mass
uses the Keplerian rotation in the circumstellar disk around the protostar. To
apply the Keplerian rotation method to a protostar at the low-mass end, we
have observed the Class 0 protostar IRAS 16253-2429 using the Atacama Large
Millimeter/submillimeter Array (ALMA) in the 1.3 mm continuum at an angular
resolution of $0\farcs 07$ (10 au), and in the 12CO, C18O, 13CO ($J=2-1$), and
SO ($J_{N}=6_{5}-5_{4}$) molecular lines, as part of the ALMA Large Program
Early Planet Formation in Embedded Disks (eDisk). The continuum emission
traces a non-axisymmetric, disk-like structure perpendicular to the associated
12CO outflow. The position-velocity (PV) diagrams in the C18O and 13CO lines
can be interpreted as infalling and rotating motions. In contrast, the PV
diagram along the major axis of the disk-like structure in the 12CO line
allows us to identify Keplerian rotation. The central stellar mass and the
disk radius are estimated to be $\sim$0.12-0.17 $M_{\sun}$ and $\sim$13-19 au,
respectively. The SO line suggests the existence of an accretion shock at a
ring ($r\sim 28$ au) surrounding the disk and a streamer from the eastern side
of the envelope. IRAS 16253-2429 is not a proto-BD but has a central stellar
mass close to the BD mass regime, and our results provide a typical picture of
such very low-mass protostars.

Circumstellar disks (235) — Protostars (1302) — Low mass stars (2050)

††facilities: ALMA††software: astropy (Astropy Collaboration et al., 2013,
2018), CASA (McMullin et al., 2007), eDisk script (Tobin et al., 2023),
ptemcee (Foreman-Mackey et al., 2013; Vousden et al., 2016), SLAM (Aso & Sai,
2023)

##  1 Introduction

The low-mass end of star formation is connected to brown dwarf (BD) formation.
Investigating such low-mass regimes is crucial to comprehensively
understanding star formation. BDs are characterized by their masses that are
not enough to fuse hydrogen ($M_{*}<0.08~{}M_{\sun}$) and are as numerous as
hydrogen-burning stars (Chabrier, 2002). The formation process of BDs can be
similar to or different from those of low-mass stars (Padoan & Nordlund,
2004). Theoretical studies have suggested various mechanisms for BD formation,
such as turbulent fragmentation of a molecular cloud (e.g., Bate, 2012), disk
fragmentation (e.g., Stamatellos & Whitworth, 2009), ejection from multiple
young stellar systems (e.g., Basu et al., 2012), photo-erosion of a prestellar
core by OB stars (e.g., Whitworth & Zinnecker, 2004), and eroding outflows
(e.g., Machida et al., 2009). A point of view for testing these scenarios is
whether or not physical quantities in star formation and BD formation follow
the same scaling laws. Kim et al. (2019) reported that their sample of 15
proto-BD candidates have different scaling laws, than the laws among $\sim 60$
low-mass protostars, between the outflow force versus the luminosity and
between the outflow force and the envelope mass. Protostars around the BD
threshold i.e., proto-BDs or very low-mass protostars, are then important to
determine down to which mass such a scaling law of protostars holds. Recent
observational studies have aimed to identify and characterize proto-BDs, as
well as pre-BDs (de Gregorio-Monsalvo et al., 2016; Huélamo et al., 2017;
Santamaría-Miranda et al., 2021). To study star formation in this very low
mass regime, it is, therefore, necessary to establish a method for precisely
estimating the central mass of each protostar down to the mass regime around
$M_{*}\sim 0.1~{}M_{\sun}$.

The most direct method for estimating a central protostellar mass is to
identify the Keplerian rotation in a circumstellar disk, as demonstrated in
previous observations toward young stellar objects in the typical low mass
regime ($M_{*}\gtrsim 0.2~{}M_{\sun}$; Yen et al., 2017). In contrast,
previous studies of two representative proto-BD candidates, IC348-SMM2E and
L328-IRS, estimated masses in indirect methods, which do not verify that the
rotational velocity has the radial profile of the Keplerian rotation but
merely assume that the observed gas motion is the Keplerian rotation or use
other kinematics. The central mass of IC348-SMM2E is estimated to be $\sim
0.002-0.024~{}M_{\sun}$ from its luminosity, mass accretion rate, and
efficiencies of mass accretion (Palau et al., 2014). This range includes a
dynamical mass of $\sim 16~{}M_{\rm Jup}$ estimated on the assumption that the
velocity gradient in the C18O $J=2-1$ line is due to the Keplerian rotation,
based on Submillimeter Array (SMA) observations (Palau et al., 2014). The
central mass of L328-IRS is estimated to be $\sim 0.012-0.023~{}M_{\sun}$ from
its mass accretion rate converted from its outflow force (Lee et al., 2018).
The outflow force is estimated using ALMA observations in the 12CO $J=2-1$
line and converted to the mass accretion rate, assuming an entrainment
efficiency of 0.25, a ratio of mass loss and accretion rates of 0.1, and a
wind velocity of $15~{}{\rm km~{}s^{-1}}$. Meanwhile, Lee et al. (2018)
explained the C18O and 13CO $J=2-1$ emission in L328-IRS by the Keplerian
rotation and suggested a central mass, $\sim 0.27$-$0.30~{}M_{\sun}$,
significantly larger than the BD mass regime. In addition to the disregard of
Keplerian-disk identification, numerical simulations predict that the disk
mass can be comparable to the central stellar mass during the early phase with
$M_{*}\lesssim 0.1~{}M_{\sun}$ (Machida et al., 2010, 2014). Such a massive
disk may cause the rotational velocity to deviate from the Keplerian rotation
determined only by $M_{*}$, preventing us from directly estimating $M_{*}$ in
observations. It is thus crucial to verify whether the direct method with the
Keplerian disk identification can be applied for protostars down to the
$M_{*}\sim 0.1$-$M_{\sun}$ regime as is for the typical low mass protostars.

###  1.1 Target IRAS 16253-2429

The Class 0 protostar IRAS 16253-2429 (hereafter I16253) in the Ophiuchus
star-forming region is a good target for the verification in that mass regime.
I16253 has a bolometric luminosity of $L_{\rm bol}=0.16~{}L_{\sun}$ and a
bolometric temperature of $T_{\rm bol}=42$ K (Ohashi et al., 2023). Its
internal luminosity is estimated to be $L_{\rm int}\sim 0.08~{}L_{\sun}$ from
the infrared luminosity of $L_{\rm IR}=0.046~{}L_{\sun}$ measured in
$1.25-70~{}\micron$ and an empirical ratio of $L_{\rm int}/L_{\rm IR}\sim 1.7$
(Dunham et al., 2008). This protostar is thus classified as a very low
luminosity object (VeLLO), which suggests that this protostar may eventually
evolve into a very-low mass star or a BD. We observed this protostar as a part
of the ALMA large program “Early Planet Formation in Embedded Disks” (eDisk)
(see the overview paper by Ohashi et al., 2023). The main goal of the program
is to reveal signs of planet formation in disks in the course of protostellar
evolution. In addition to this main goal, I16253 was observed to investigate
the protostellar evolution and the disk formation down to the BD mass regime.

The central stellar mass of this target in previous works has been estimated
to be a wide range of $\sim 0.02$ to $\sim 0.12~{}M_{\sun}$, depending on
methods. Tobin et al. (2012), who observed I16253 in the N2H+
$J,F_{1},F=1,2,2-0,1,1$ line using the Combined Array for Research in
Millimeter-wave Astronomy (CARMA) at a resolution of $\sim 9\arcsec$,
estimated the stellar mass to be $0.1~{}M_{\sun}$ based on comparisons of a
position-velocity (PV) diagram across the associated outflow on a 6000-au
scale between the observations and a model of the rotating, collapsing
envelope (UCM model; Ulrich, 1976; Cassen & Moosman, 1981). Yen et al. (2017)
estimated the stellar mass to be $0.02^{+0.02}_{-0.01}~{}M_{\sun}$ based on
C18O $J=2-1$ observations of the protostar using ALMA at an angular resolution
of $\sim 1\arcsec$; they reproduced the observed PV diagram along the major
and minor axes of the C18O envelope with a model envelope having a rotation
conserving its specific angular momentum and a radial free-fall motion. Hsieh
et al. (2019) used a similar model of the envelope and estimated the stellar
mass to be $\sim 0.028~{}M_{\sun}$ by reproducing PV diagrams of 12CO and C18O
$J=2-1$ obtained with ALMA at an angular resolution of $0\farcs 1-0\farcs 4$.
Additionally, Hsieh et al. (2019) also suggested the stellar mass of
$0.12~{}M_{\sun}$ with the assumption that the 12CO emission arises from a
Keplerian disk; they found offsets of the 12CO emission, $\pm 7.8$ au away
from the center, at two velocity channels of $V_{\rm LSR}-V_{\rm sys}=\pm
3.4~{}{\rm km~{}s^{-1}}$ by 2D Gaussian fitting. This mass derived on the
assumption of Keplerian rotation can be different from the one derived on the
assumption of free-fall motion, if the infall velocity is slower than the
free-fall velocity, as reported in other protostars (Ohashi et al., 2014; Aso
et al., 2015, 2017; Sai et al., 2022). These previous works show controversy
about the central stellar mass of I16253 around the BD mass regime (i.e.,
$M_{*}<0.08~{}M_{\sun}$).

Previous observations using the Walraven Photometer estimated the distance of
the Ophiuchus star-forming region to be $125\pm 25$ pc (de Geus et al., 1989),
while the same observation estimated the distance at a location of
$(l,b)=(353.161\arcdeg,15.936\arcdeg)$ closest to I16253
$(353.2\arcdeg,16.5\arcdeg)$ to be 140 pc. This distance is consistent with an
estimate using Hipparcos ($120-160$ pc; Knude & Hog, 1998) and a new estimate
by Gaia DR2, $139^{+9}_{-10}$ pc, at a location of
$(353.2\arcdeg,16.6\arcdeg)$ closest to I16253 (Zucker et al., 2020). In
addition to the I16253 distance, the average distance of the whole Ophiuchus
region was updated by Gaia DR2 to $144\pm 9$ pc (Zucker et al., 2019). We
adopt the Gaia DR2 distance of 139 pc as the distance of I16253 in this paper.

This work aims to estimate the central mass of the Class 0 protostar I16253
more accurately and precisely than the previous works and reveal the kinematic
structures around the protostar. Our observations have $>2.5$ times better
sensitivity in the CO isotopologue lines than that of previous observations in
the same lines at similar angular and velocity resolutions to those of our
observations. We aim to reveal kinematic structures, such as a rotating disk,
an infalling envelope, and an outflow, using the molecular line emission
detected at the higher signal to noise ratio. The rest of the paper has the
following structure. Section 2 describes key points of the settings and data
processing of our I16253 observations. Section 3 shows the observed results in
the 1.3-mm dust continuum and the 12CO, 13CO, C18O, and SO lines. We analyze
the asymmetry in the continuum image and radial profiles of the rotational
velocity in Section 4. The rotational velocity is found to follow the
Keplerian rotation, which allows us to directly estimate the central mass of
I16253. We discuss kinetic structures around I16253, showing a schematic
picture (Figure 18) of those structures, in Section 5. A summary of our
results and interpretation is in Section 6.

##  2 Observations

The details of the observing strategy, spectral setups, and data reduction
process are presented in Ohashi et al. (2023)111The scripts used for
reduction, including the self-calibration, can be found at
https://github.com/jjtobin/edisk (Tobin et al., 2023).. Here we briefly
describe the key points specific to I16253 and summarize the observational
parameters in Table 1. We observed I16253 using ALMA in Cycle 8 with the
antenna configuration of C-8 on 2021 October 5, 26, 27, and 28
(2019.1.00261.L). The total observing time with C-8 is $\sim 310$ min (5.2
hr), while the on-source observing time is $\sim 60$ min. The number of
antennas with C-8 was 46, and the projected baseline length ranges from 52 to
10540 m. Any Gaussian component with FWHM$\gtrsim 2.8\arcsec$ ($\sim 400$ au)
was resolved out by $\gtrsim 63\%$ with this shortest baseline length (Wilner
& Welch, 1994). The phase center is $(\alpha_{\rm ICRS},\delta_{\rm
ICRS})=(16^{\rm h}28^{\rm m}$21$\fs$6, $-24^{\circ}36\arcmin 23\farcs 4)$.

To recover the more extended structure, we also observed I16253 using ALMA in
Cycle 8 with the antenna configuration of C-5 on 2022 June 14 and 15
(2019.A.00034.S). The total observing time with C-5 is $\sim 170$ min (2.8
hr), while the on-source observing time is $\sim 30$ min. The number of
antennas with C-5 was 43, and the projected baseline length ranges from 14 to
1290 m. The resolving out scale is $10\arcsec$ ($\sim 1400$ au) with this
shortest baseline length. The phase center is the same as with C-8.

The ALMA observations were set up to cover spectral windows including CO
isotopologues ($J=2-1$), SO ($J_{N}=6_{5}-5_{4}$), and other molecular lines
at Band 6. The spectral window for the 12CO ($J=2-1$) line has 3840 channels
covering a 940 MHz bandwidth at an original frequency resolution of 224 kHz.
Spectral windows for the 13CO, C18O ($J=2-1$), SO ($J_{N}=6_{5}-5_{4}$), and
H2CO ($3_{21}-2_{20}$) lines have 960 channels covering a 59 MHz bandwidth at
an original frequency resolution of 61 kHz. Spectral windows for the other
lines have 3840 channels covering a 1.9 GHz bandwidth at an original frequency
resolution of 488 kHz. Channels were binned to produce the velocity
resolutions of 0.32, 0.17, 0.17, $0.17~{}{\rm km~{}s^{-1}}$, respectively, to
make maps in the 13CO, C18O, SO, and H2CO lines. The velocity resolution for
other lines is $1.34~{}{\rm km~{}s^{-1}}$ (Appendix A). Continuum maps were
made using line-free channels of the spectral windows including two wide
spectral windows with a $\sim 1.8$ GHz bandwidth centering at $\sim 220$ to
$\sim 230$ GHz. The absolute flux density accuracy of ALMA is $\sim 10\%$ in
this frequency band.

All the imaging procedure was carried out with Common Astronomical Software
Applications (CASA) (McMullin et al., 2007) version 6.2.1. The visibilities
were Fourier transformed and cleaned with Briggs weighting and a robust
parameter of $-2.0$, 0.0, 0.5, or 2.0, using the CASA task tclean at a pixel
size of $0\farcs 02$. Continuum images adopt robust$=-2.0$, 0.0, and 2.0 to
show the most compact, intermediate, and most extended components,
respectively. The line images adopt robust$=0.5$ and 2.0 to show compact
components and extended components, respectively. The line images do not have
signal-to-noise (S/N) ratios high enough for robust$=-2$. The continuum images
with robust$=-2$ and 2 are produced with tapering parameters of 3 and 1
M$\lambda$ ($\sim 0\farcs 06$ and $\sim 0\farcs 18$), respectively. The
tapering parameter of 3 M$\lambda$ was selected to focus more on extended
components with robust$=2$, while the tapering parameter of 1 M$\lambda$ was
selected to increase the S/N ratio of the image with robust$=-2$. All the line
images are produced with a tapering parameter of 2 M$\lambda$ ($\sim 0\farcs
09$) to increase the S/N ratio. The resultant angular resolution is $0\farcs
04-0\farcs 29$ for the continuum emission and $0\farcs 17-0\farcs 38$ for the
line emission. We also performed self-calibration for the continuum data using
tasks in CASA (tclean, gaincal, and applycal). While all the maps are primary-
beam corrected in this paper with the primary beam size of $26\arcsec$, the
root-mean-square noise levels of the line maps were measured in emission-free
channels/areas before the primary-beam correction, which indicates the
sensitivity achieved toward the phase center.

Table 1: Summary of the parameters of our ALMA observations toward the Class 0
protostar I16253.

Date | 2021 October 5, 26, 27, & 28 | 2022 June 14 & 15
---|---|---
Projected baseline range | 52–10540 m | 14–1290 m
Maximum Recoverable Scale | ${0\farcs 62}$ | ${2\farcs 9}$
Bandpass/flux calibrator | J1517-2422 | J1517-2422
Check source | J1650-2943 | J1650-2943
Phase calibrator | J1633-2557 | J1700-2610
| Continuum |  12CO $J=2-1$ |  13CO $J=2-1$ | C18O $J=2-1$ | SO $J_{N}=6_{5}-5_{4}$
Frequency (GHz) | 225 | 230.5380000 | 220.3986842 | 219.5603541 | 219.9494420
Freq./vel. width |  $\sim 2$ GHz | $0.32~{}{\rm km~{}s^{-1}}$ | $0.17~{}{\rm km~{}s^{-1}}$ | $0.17~{}{\rm km~{}s^{-1}}$ | $0.17~{}{\rm km~{}s^{-1}}$
Beam (P.A.) | $0\farcs 29\times 0\farcs 24\ (86\arcdeg)^{a}$ | $0\farcs 35\times 0\farcs 25\ (76\arcdeg)^{a}$ | $0\farcs 37\times 0\farcs 26\ (77\arcdeg)^{a}$ | $0\farcs 38\times 0\farcs 27\ (78\arcdeg)^{a}$ | $0\farcs 37\times 0\farcs 26\ (78\arcdeg)^{a}$
noise rms | $23~{}{\rm\mu Jy~{}beam^{-1}}^{a}$ | $1.4~{}{\rm mJy~{}beam^{-1}}^{a}$ | $2.5~{}{\rm mJy~{}beam^{-1}}^{a}$ | $1.9~{}{\rm mJy~{}beam^{-1}}^{a}$ | $2.3~{}{\rm mJy~{}beam^{-1}}^{a}$
Beam (P.A.) |  | $0\farcs 073\times 0\farcs 048\ (77\arcdeg)^{c}$
---
$0\farcs 043\times 0\farcs 034\ (73\arcdeg)^{d}$
$0\farcs 17\times 0\farcs 13\ (88\arcdeg)^{b}$ | $0\farcs 17\times 0\farcs 13\ (87\arcdeg)^{b}$ | $0\farcs 18\times 0\farcs 14\ (87\arcdeg)^{b}$ | $0\farcs 18\times 0\farcs 14\ (87\arcdeg)^{b}$
noise rms |  | $22~{}{\rm\mu Jy~{}beam^{-1}}^{c}$
---
$100~{}{\rm\mu Jy~{}beam^{-1}}^{d}$
$1.8~{}{\rm mJy~{}beam^{-1}}^{b}$ | $3.0~{}{\rm mJy~{}beam^{-1}}^{b}$ | $2.1~{}{\rm mJy~{}beam^{-1}}^{b}$ | $2.5~{}{\rm mJy~{}beam^{-1}}^{b}$

Note. — The beam and noise rms values are shown for different robust
parameters with superscripts of (a) Robust$=2$, (b) 0.5, (c) 0, and (d) $-2$.

##  3 Results

###  3.1 1.3 mm continuum

Figure 1 shows continuum images with different robust and tapering parameters.
Figure 1(a) shows that the continuum emission with the robust parameter of 0
traces a disk-like structure on a $\sim 40$ au scale. This emission appears
more extended to the southeast than to the northwest from the emission peak.
This asymmetry is investigated in detail in Section 4.1. Two dimensional (2D)
Gaussian fitting to this robust$=0$ continuum image provides the central
position of $(\alpha_{\rm ICRS},\ \delta_{\rm ICRS})=(16^{\rm h}28^{\rm
m}21\fs 615,\ -24\arcdeg 36\arcmin 24\farcs 33$) with an uncertainty of (0.6
mas, 0.3 mas). The fitting, including the error bar calculation, used the CASA
task imfit. The fitting result also provides the deconvolved size of $107\pm
2$ mas $\times$ $40\pm 2$ mas with a major axis in the direction of P.A. of
$113\arcdeg\pm 1\arcdeg$. We adopt this direction as the major axis of the
protostellar system, whici is perpendicular to the outflow direction,
P.A.$=200\arcdeg-205\arcdeg$ (Yen et al., 2017). The aspect ratio corresponds
to an inclination angle of ${\rm arccos}(40/107)\sim 68\arcdeg$, assuming a
geometrically thin disk.

The peak intensity and flux density derived from the Gaussian fitting are
$5.10\pm 0.02~{}{\rm mJy~{}beam^{-1}}$ and $11.5\pm 0.2$ mJy, respectively.
The flux densities integrated over $0\farcs 4\times 0\farcs 4$ (the compact
emission) and $10\arcsec\times 10\arcsec$ (the maximum recoverable scale of
our ALMA observations) regions are $12.4\pm 0.3$ and $25\pm 2$ mJy,
respectively. The flux density $F_{\nu}$ can be converted to a lower limit of
gas mass, assuming that the emission is optically thin: $M_{\rm
gas}=d^{2}F_{\nu}/\kappa_{\nu}B_{\nu}(T)$, where $d$, $\kappa_{\nu}$, and
$B_{\nu}$ are the target distance, dust mass opacity, and the Planck function
with an average temperature $T$, respectively. For this conversion, two
different temperatures are adopted: $T=20$ and 27 K. $T=20$ K is derived from
the relation between 850-$\micron$ continuum fluxes and model gas masses for
44 young stellar objects (Andrews & Williams, 2005). $T=27$ K is calculated
from the empirical relation, $T=43~{}{\rm K}~{}(L_{\rm
bol}/1~{}L_{\sun})^{0.25}$, from Tobin et al. (2020). With these temperatures,
$d=139$ pc, $\kappa_{\nu}=2.3~{}{\rm cm}^{2}~{}{\rm g}^{-1}$ (Beckwith et al.,
1990), and a gas-to-dust mass ratio of 100, $F_{\nu}\simeq 12$ mJy corresponds
to a gas mass of $M_{\rm gas}\simeq(1.4-2.1)\times 10^{-3}~{}M_{\sun}$.

Figure 1(b) adopts the robust parameter of $-2$ and the 3-M$\lambda$ taper to
focus on the central region with a reasonable S/N ratio. The plotted spatial
range is half of Figure 1(a). This image shows an extension or a secondary
peak to the southeast, separated from the central peak by $\sim 0\farcs 1$
roughly along the major axis. This is consistent with the extension seen in
the image with robust$=0$ (Figure 1a).

Figure 1(c), adopting the robust parameter of 2 and the 1-M$\lambda$ taper,
shows the extended emission over $\sim 600$ au from the continuum emission
peak to the northwest, as well as the compact ($\sim 100$ au) emission. This
large-scale structure is also reported in Yen et al. (2017).

Figure 1: Images of the 1.3 mm continuum emission produced with (a) the robust
parameter of 0.0 and no taper, (b) the robust parameter of $-2.0$ and a 3
M$\lambda$ ($\sim 0\farcs 06$) taper, (c) the robust parameter of 2.0 and a 1
M$\lambda$ ($\sim 0\farcs 18$) taper. The brightness temperature $T_{b}$ is
calculated from the Rayleigh-Jeans approximation. The contour levels are
$3,6,12,24,48,96,192\sigma$, where $1\sigma$ is (a) 0.022, (b) 0.100, (c)
$0.023~{}{\rm mJy~{}beam^{-1}}$ (Table 1). The filled ellipses at the bottom-
left corners denote the synthesized beams: (a) $0\farcs 073\times 0\farcs 048$
($77\arcdeg$), (b) $0\farcs 043\times 0\farcs 034$ ($73\arcdeg$), and (c)
$0\farcs 29\times 0\farcs 24$ ($86\arcdeg$). The diagonal lines denote the
major and minor axes of the emission in the image in panel (a):
P.A.$=113\arcdeg$ and $23\arcdeg$. The square in panel (a) shows the plotted
range of panel (b).

###  3.2 12CO $J=2-1$

Figure 2(a) shows the integrated intensity (moment 0) and the mean velocity
(moment 1) maps in the 12CO emission produced with the robust parameter of 2.
Positive and negative intensities are integrated for the moment 0 map, while
intensities above the $3\sigma$ level are integrated for the moment 1 map. The
noise level of each moment 0 map is calculated by the noise propagation from
the noise level of the data cube (Table 1) used to make the moment 0 map:
$\sigma_{\rm mom0}=\sigma_{\rm cube}dv\sqrt{N_{\rm ch}}$, where $\sigma_{\rm
mom0}$, $\sigma_{\rm cube}$, $dv$, and $N_{\rm ch}$ are the noise level of the
moment 0 map, the noise level of the cube data, the velocity width per channel
(Table 1), and the number of integrated channels, respectively. The moment 0
and 1 maps of the other lines in this paper are also made in the same method.
The 12CO emission traces a clear bipolar outflow whose axis is perpendicular
to the major axis of the disk-like structure. An additional eastern component
results from a large-scale structure around the systemic velocity, which is
velocity resolved out by the interferometric observation. While the emission
in the northern outflow lobe is overall blueshifted, redshifted emission can
also be found on the northern side near the outflow axis. Similarly, the
southern lobe is mostly redshifted but partly blueshifted. This is expected if
a given outflow lobe crosses the plane of the sky that contains the central
source, with the outflow axis close to the plane of sky (Cabrit, 1989). A part
of the lobe in the front of the plane of the sky appears blue-shifted, while
the other part of the lobe behind the plane of the sky appears red-shifted.
Using a radially-expanding parabolic outflow model, Yen et al. (2017)
estimated the outflow shape and inclination by fitting a moment 0 map and a PV
diagram along the outflow axis of the 12CO emission observed at an angular
resolution of $\sim 1\arcsec$. Their best model provides an inclination angle
of $60\arcdeg-65\arcdeg$ ($0\arcdeg$ means pole-on) and an opening angle of
$60\arcdeg-70\arcdeg$, calculated from the parabolic shape at the distance
from the protostar along the outflow axis of $z\sim 5\arcsec-10\arcsec$. This
inclination is consistent with that of the disk-like structure detected in the
dust continuum emission (Section 4.1). Thus, their model is consistent with
our observations which reveal blue- and red-shifted outflow components to the
northeast and southwest, respectively.

Figure 2: Moment 0 (integrated intensity) and 1 (mean velocity) maps in the
12CO $J=2-1$ line. The contour maps show the moment 0 map, while the color
images show the moment 1 map. (a) Moment maps with the robust parameter of
2.0, integrated from $V_{\rm LSR}=-6.0$ to $14.0~{}{\rm km~{}s^{-1}}$. (b)
Zoom-in view of the moment maps with the robust parameter of 0.5, integrated
from $V_{\rm LSR}=-6.0$ to $14.0~{}{\rm km~{}s^{-1}}$. The contour levels are
in (a) $12\sigma$ steps from $12\sigma$ with $1\sigma=3.2~{}{\rm
mJy~{}beam^{-1}}~{}{\rm km~{}s^{-1}}$ and (b) $6\sigma$ steps from $6\sigma$
with $1\sigma=3.5~{}{\rm mJy~{}beam^{-1}}~{}{\rm km~{}s^{-1}}$. The filled
ellipse at the bottom-left corner shows the synthesized beam: (a) $0\farcs
35\times 0\farcs 25$ ($76\arcdeg$) and (b) $0\farcs 17\times 0\farcs 13$
($88\arcdeg$). The diagonal lines are the major and minor axes of the
continuum emission, P.A.$=113\arcdeg$ and $23\arcdeg$. The green contour in
panel (b) shows the 5$\sigma$ level of the continuum emission (Figure 1a).

Figure 2(b) shows a zoom-in view of the moment 0 and 1 maps with a different
robust parameter of 0.5. The green contour shows the 5$\sigma$ level of the
continuum emission (Figure 1a) for comparison. The moment 0 map shows local
peaks within $\sim 0\farcs 4$ from the protostellar position, which is
reported in Hsieh et al. (2019) as a quadruple peak. In addition to the
velocity gradient along the outflow axis on a $1\arcsec$ scale, this figure
also shows a velocity gradient perpendicular to the outflow axis near the
midplane on a $0\farcs 3$ scale, i.e., transition from a blueshifted part at
P.A.$=113\arcdeg$ to a redshifted part at P.A.$=-67\arcdeg$. This velocity
gradient is expected for rotating gas and analyzied in more detail in Section
4.2.

Figure 4 shows a position-velocity (PV) diagram of the 12CO emission along the
outflow axis from the robust parameter 2 images and a width of $0\farcs 6$
($\sim 2\times$ beam size). The width allows us to increase the S/N ratio but
still focus on the velocity structure on the outflow axis. Because of this
width, the noise level of this PV diagram is better than the one in Table 1.
This figure also shows that both northern and southern lobes have both blue-
and red-shifted parts. The outflow velocity is $|V-V_{\rm sys}|\sim
1.5-2.5~{}{\rm km~{}s^{-1}}$ at any position except for the central
$1\arcsec-2\arcsec$ region. The central component appears not to be a part of
the outflow because of its $\sim 3$ times larger line widths. The outer
($>2\arcsec$) emission shows several ($\sim 5-9$) local peaks or extensions
within $10\arcsec$ ($\sim 1400$ au) in each lobe, based on visual inspection
(short lines in Figure 4). This may be suggestive of episodic mass ejection
(Lee, 2020). The interval angular scale, $10\arcsec/7=1.4\arcsec$, is
significantly larger than the beam size of $0\farcs 35$. With the inclination
angle of $65\arcdeg$ (Yen et al., 2017), the deprojected length $l$ and the
deprojected velocity $v$ of the outflow can be calculated as
$l=l^{\prime}/\cos i$ and $v=v^{\prime}/\sin i$, respectively, where
$l^{\prime}$ and $v^{\prime}$ are the projected length and velocity. The
interval $t$ of the mass ejection is a time required for material to move from
a peak to the next peak. This can be calculated with the number of peaks $n$
within the length: $t=l^{\prime}/n/v^{\prime}$. Then, the values measured
above provide a typical interval of $1400~{}{\rm au}/\cos 65\arcdeg/(5\ {\rm
to}\ 9)/(1.5\ {\rm to}\ 2.5~{}{\rm km~{}s^{-1}}/\sin 65\arcdeg)=600\ {\rm to}\
2000$ year. The episodic mass ejection is thought to be caused by episodic
accretion bursts. The interval of accretion bursts is investigated in
statistical studies using FUor objects. Park et al. (2021) reported accretion-
burst intervals of $500-2300$ yr in a statistical study of FUor outbursts
based on the fraction of outbursting sample in a set of NEOWISE observations.
The fraction of variable sources is similar between VeLLOs like I16253 and
more luminous objects in their sample. The uncertainty comes from the number
of outbursts ranging from 2 to 9 and the total number of sample protostars
ranging from 735 to 1059. This accretion-burst interval is consistent with the
mass-ejection interval of I16253, implying that accretion bursts likely
occurred in I16253.

Figure 3:  Figure 4: Position-velocity diagram in the 12CO $J=2-1$ line along
the outflow axis (P.A.$=23\arcdeg$) with a width of $0\farcs 6$. The positive
offset corresponds to the northern side. The contour levels are in $6\sigma$
steps from $6\sigma$, where $1\sigma$ is $0.85~{}{\rm mJy~{}beam^{-1}}$. The
dashed lines denote the systemic velocity ($V_{\rm LSR}=4~{}{\rm
km~{}s^{-1}}$) and the protostellar position. The orange short lines denote
the position of local peaks or extensions visually identified.

This PV diagram shows a strong emission component at offsets of $<-10\arcsec$.
This component is not confined around the outflow axis but extended over the
redshifted lobe. This may imply that there could be more material on the
southern side than on the northern side, although our observation did not
target such a large spatial scale.

###  3.3 C18O and 13CO $J=2-1$

Figure 5 shows the moment 0 and 1 maps in the C18O and 13CO lines in
robust$=0.5$ and 2. The two lines show an overall structure extended along the
major axis (P.A.$=113\arcdeg$) and a velocity gradient primarily along the
same direction. The integrated intensity maps with robust$=0.5$ show double
peaks with a separation of $\sim 0\farcs 2-0\farcs 3$, and the integrated
intensity is stronger at the eastern, blueshifted peaks than at the western,
redshifted peaks, by $3\sigma$ in both lines. The integrated intensity maps
with robust$=2$ show more clearly that the eastern emission is stronger than
the western emission. The detected emission size is larger in the C18O line
than in the 13CO line with both robust parameters. This is because the 13CO
emission tends to be optically thicker and more extended than the C18O
emission, and thus the 13CO emission is resolved out more severely around the
systemic velocity.

Figure 5: Moment 0 and 1 maps in the C18O $J=2-1$ line emission and the 13CO
$J=2-1$ line emission. The contour maps show the moment 0 map, while the color
images show the moment 1 map. The contour levels are in $3\sigma$ steps in the
robust$=0.5$ images and $6\sigma$ steps in the robust$=2$ images starting from
$3\sigma$. The C18O emission is integrated from 0.8 to $7.1~{}{\rm
km~{}s^{-1}}$. These moment 0 maps have $1\sigma$ of 2.2 (robust$=0.5$) and
$2.0~{}{\rm mJy~{}beam^{-1}}~{}{\rm km~{}s^{-1}}$ (robust$=2$). The 13CO
emission is integrated from 0.7 to $7.3~{}{\rm km~{}s^{-1}}$. These moment 0
maps have $1\sigma$ of 2.8 (robust$=0.5$) and $2.7~{}{\rm
mJy~{}beam^{-1}}~{}{\rm km~{}s^{-1}}$ (robust$=2$). The beam sizes are (a)
$0\farcs 18\times 0\farcs 14$ ($87\arcdeg$), (b) $0\farcs 38\times 0\farcs 27$
($78\arcdeg$) (c) $0\farcs 17\times 0\farcs 13$ ($87\arcdeg$), and (d)
$0\farcs 37\times 0\farcs 26$ ($77\arcdeg$). The diagonal lines are the major
and minor axes of the continuum emission, P.A.$=113\arcdeg$ and $23\arcdeg$.

###  3.4 SO $J_{N}=6_{5}-5_{4}$

Figure 6(a) shows the integrated intensity and mean velocity maps in the SO
line produced with robust$=0.5$ to focus on the central compact emission. The
SO emission shows a double-peaked structure with a separation of $\sim 0\farcs
2-0\farcs 3$ and the velocity gradient along the major axis, like the C18O and
13CO lines. On the other hand, the mean velocity of the SO emission is larger
outside the peaks than inside the peaks, unlike the C18O and 13CO lines.

Figure 6(b) shows the maps produced with robust$=2$ to focus on extended
emission. The central double peaks are not spatially resolved in this map. The
velocity structure in the central $1\arcsec$ region is overall the same as
that in Figure 6(a): eastern blueshifted emission and western redshifted
emission. In addition, this figure shows an extended structure to the east
from the central protostar. This structure is extended ($\sim 3\arcsec$)
beyond the C18O and 13CO emission ($\sim 1\arcsec$) in Figure 5. An inner part
($\lesssim 1\arcsec$ from the protostellar position) of this extended
structure shows velocities blueshifted from the systemic velocity, $V_{\rm
LSR}=4~{}{\rm km~{}s^{-1}}$, which is consistent with the C18O and 13CO
velocities. In contrast, the outer part ($\gtrsim 1\arcsec$) shows slightly
redshifted velocities ($V_{\rm LSR}-V_{\rm sys}\sim 0.3~{}{\rm km~{}s^{-1}}$).
This redshifted velocity on the eastern side is different from the velocity
gradient seen in the C18O and 13CO lines and also different from the velocity
in the eastern cavity wall of the 12CO outflow (Figure 2a). This SO component
is discussed in more detail in Section 5.2.

Figure 6: Moment 0 and 1 maps in the SO $J_{N}=6_{5}-5_{4}$ line emission. The
contour maps show the moment 0 map, while the color images show the moment 1
map. (a) The robust parameter is 0.5. The emission is integrated from 1.7 to
$6.3~{}{\rm km~{}s^{-1}}$. The contour levels are in $3\sigma$ steps from
$3\sigma$ with $1\sigma=2.3~{}{\rm mJy~{}beam^{-1}}~{}{\rm km~{}s^{-1}}$. The
beam size is $0\farcs 18\times 0\farcs 14$ ($87\arcdeg$). (b) The robust
parameter is 2. The emission is integrated from 1.7 to $6.3~{}{\rm
km~{}s^{-1}}$. The contour levels are in $3\sigma$ steps from $3\sigma$ until
$15\sigma$ and then in $10\sigma$ steps with $1\sigma=2.1~{}{\rm
mJy~{}beam^{-1}}~{}{\rm km~{}s^{-1}}$. The beam size is $0\farcs 37\times
0\farcs 26$ ($78\arcdeg$). The diagonal lines are the major and minor axes of
the continuum emission, P.A.$=113\arcdeg$ and $23\arcdeg$.

##  4 Analysis

###  4.1 Non-axisymmetry of the continuum image

The 1.3 mm continuum emission with robust$=0$ (Figure 1a) is more extended in
the southeastern side ($\sim 30$ au) than in the northwestern side ($\sim 20$
au). We verify that the extended components are real by comparing images made
with various robust parameters. To extract the extended component, an
axisymmetric model of the dust continuum emission is constructed in this
subsection. The model is made from a radial profile of intensities in the unit
of Jy pixel-1, before beam convolution, set with free parameters of
$I(r=0),I(dr),I(2dr),...,I(11dr)$. The grid separation $dr$ is half of the
beam minor-axis. The largest radius $11dr=0\farcs 26$ covers the entire
continuum emission with robust$=0$. Then, the circular 2D intensity
distribution is projected by the inclination factor of $\cos i$ in the
direction of the minor axis and rotated by the position angle $pa$; these two
angles are also free parameters. This elliptical 2D intensity in the unit of
Jy pixel-1 is convolved with the observational beam to make the model image in
the unit of ${\rm Jy~{}beam^{-1}}$. Hence, this model does not assume any
shape, such as a Gaussian function but only the axisymmetry. This modeling
method follows the one in Aso et al. (2021). The central position is not a
free parameter but fixed at the observed peak position, $(\alpha_{\rm ICRS},\
\delta_{\rm ICRS})=(16^{\rm h}28^{\rm m}21\fs 6153,\ -24\arcdeg 36\arcmin
23\farcs 325)$, which coincides with the Gaussian center derived in Section
3.1. The best-fit parameters are obtained by minimizing $\chi^{2}$ between the
observed and model images, divided by the number of pixels in the beam,
through the Markov chain Monte Carlo (MCMC) methods with the public Python
package ptemcee. The numbers of free parameters, steps, and walkers per
parameter are 14, 8000, and 16, respectively. The first half (4000) steps were
removed for the burn-in. The derived best-fit angles are $i=64.1\arcdeg\pm
0.5\arcdeg$ and $pa=114.3\arcdeg\pm 0.6\arcdeg$, and Figure 8 shows the
posterior distribution of these two angles produced by the MCMC fitting. The
uncertainties of these angles are derived as the 16 and 84 percentiles.

Figure 7:  Figure 8: Corner plot of the MCMC fitting to the continuum image.
Two angles, $i$ and $pa$, are plotted among the 14 free parameters. The dashed
lines show the 16, 50, and 84 percentiles.

Figure 9(a) shows the comparison between the best-fit axisymmetric model and
the observed image. The overall structure in the observed image is reproduced,
but the model emission is weaker than the observed emission in the southeast,
while it is stronger in the northwest, as expected. The contour map in Figure
9(b) shows the residual image after subtracting the model from the observed
image with robust$=0$, while the color map shows the observed continuum image
with robust$=-2$ (same as Figure 1c). The strong residual is located in the
southeast. This residual overlaps with the extension in the robust$=-2$ image,
supporting that this excess from the symmetric component is not due to the
specific robust parameter. The direction of the residual and extension
coincides with the direction of the SO extended structure and the stronger
peak on the eastern side in the C18O and 13CO moment 0 maps.

Figure 9: (a) Axisymmetric model (red) fitted to the observed 1.3 mm continuum
emission (black) with the robust parameter of 0.0. The contour levels are the
same as Figure 1(b): 3,6,12,24,48,96,192$\sigma$. (b) Residual (contour map)
between the model and observation in panel (a), overlaid on the robust$=-2$
continuum image (color map; Figure 1c). The contour levels are the same:
$3,6,12,24\sigma$.

###  4.2 Keplerian rotation in the 12CO emission

Previous observational studies have not yet identified any part of gas motion
in I16253 having a radial profile of the Keplerian rotation but assumed that
the observed gas follows Keplerian rotation. The Keplerian disk is crucial to
estimate the central stellar mass in the protostellar phase and verify whether
I16253 is a proto-BD or a very low-mass protostar. To tackle this problem, we
analyze the position-velocity (PV) diagrams along the major axis. Figure 10(a)
shows the PV diagrams in the 12CO (blue contours), C18O (red contours), and
13CO (color) lines. The C18O and 13CO lines show a typical shape of an
infalling plus rotating motion: a distorted diamond shape with emission in all
the four quadrants, as reported in Hsieh et al. (2019). In contrast, the 12CO
emission is concentrated on the first (upper right) and third (lower left)
quadrants (the western redshifted emission and the eastern blueshifted
emission), showing a clear velocity gradient, in $\pm(0\farcs 1-0\farcs 2)$ at
$|V-V_{\rm sys}|\gtrsim 2~{}{\rm km~{}s^{-1}}$. Figure 10(b) shows the PV
diagrams in the three lines along the minor axis. The C18O and 13CO emission
can be seen in all four quadrants in the minor-axis PV diagram with a diamond
shape, which is typical for the infall motion (e.g., Ruíz-Rodríguez et al.,
2022). The 12CO emission is mainly concentrated around the center, with
additional components in the second and fourth quadrants (the northern
blueshifted emission and the southern redshifted emission). While the
additional components show a velocity gradient due to the outflow motion, the
main component shows no velocity gradient along the minor axis. These PV
diagrams suggest that the 12CO line traces a purely rotating, i.e., Keplerian
disk, while the C18O and 13CO lines trace the infalling rotating envelope
previously identified. This difference among individual molecular lines can be
understood by the different strengths of emission and the different missing
fluxes in the interferometric observations. The 12CO emission is strong enough
to be detected even in the expected compact disk ($r\sim 0\farcs 1-0\farcs 2$
in the PV diagram as well as in the continuum image). The 13CO and C18O
emission is not detected sufficiently in this compact region at the high
velocities where the 12CO emission is detected. The nondetection in the 13CO
and C18O lines suggests the 12CO emission is not optically thick at the high
velocities. In contrast to the high velocities, the 12CO emission is more
extended and strongly resolved out at the low velocities where the 13CO and
C18O emission traces the envelope.

Figure 10: Position-velocity diagrams in the 12CO (blue; robust$=0.5$), C18O
(red; robust$=2.0$), and 13CO (color; robust$=2.0$) $J=2-1$ lines along the
(a) major and (b) minor axes. The positive offset corresponds to the western
and northern sides in panels (a) and (b), respectively. The contour levels are
in $6\sigma$ steps from $6\sigma$, where $1\sigma$ is $1.4~{}{\rm
mJy~{}beam^{-1}}$ for 12CO, $2.1~{}{\rm mJy~{}beam^{-1}}$ for C18O, and
$3.0~{}{\rm mJy~{}beam^{-1}}$ for 13CO. The beam sizes in these diagrams are
$\sim 0\farcs 15$ for 12CO, $\sim 0\farcs 33$ for C18O, and $\sim 0\farcs 32$
for 13CO (Table 1).

In order to verify whether the 12CO line traces the Keplerian rotation, we
first find the emission ridge (e.g., Yen et al., 2013; Aso et al., 2015; Sai
et al., 2020) and edge (e.g., Seifried et al., 2016; Alves et al., 2017) in
the major-axis PV diagram and fit them with power-law functions. The ridge is
a center of the emission along the positional axis at each velocity, while the
edge is an outer boundary of the emission along the positional axis at each
velocity. The analysis and fitting process was performed through the Python
open package pvanalysis in Spectral Line Analysis/Modeling (SLAM; Aso & Sai,
2023)222https://github.com/jinshisai/SLAM, and the detail is described in the
overview paper of the eDisk project (Ohashi et al., 2023). We thus mention
here the settings specific to the case of I16253. The ridge point is defined
as the intensity-weighted mean position $x_{m}$ at each velocity $v$:
$x_{m}(v)=\int I(x,v)xdx/\int I(x,v)dx$. $x_{m}$ is calculated using the
intensities above the $5\sigma$ level along the positional axis at each
velocity in the range of $|V-V_{\rm sys}|>2~{}{\rm km~{}s^{-1}}$. This
velocity range is selected because the emission is located in the first and
third quadrants in this velocity range. Similarly, using the 1D intensity
profile, the edge point is defined, at each velocity in the same velocity
range, as the outer position where the intensity is at the $5\sigma$ level. In
addition, the tool uses velocities higher than or equal to the velocity at
which the derived edge/ridge radius is largest so that the relation between
radius and velocity can be consistent with spin-up rotation. The ridge and
edge tend to under- and overestimate, respectively, the radius at a given
velocity (Maret et al., 2020, the ridge and edge here correspond to their
“centroid” and “first emission contour”, respectively), and thus we adopt the
former and the latter as the lower and upper limits of radius in this paper.

Figure 11 shows the estimated ridge and edge points overlaid on the PV diagram
in the linear and logarithmic scales. The logarithmic diagram is made by
avaraging the emission in the first and third quadrants of the linear diagram.
The edge and ridge radii are separately fitted with a power-law relation
between either edge or ridge radius $R$ and the velocity $V$:

$\displaystyle V=V_{b}\left(\frac{R}{R_{b}}\right)^{-p}+V_{\rm sys},$ (1)
$\displaystyle p=p_{\rm in}\ {\rm if}\ R<R_{b}\ {\rm else}\ p_{\rm in}+dp.$
(2)

This fitting uses the MCMC method with the tool of emcee, where the numbers of
walker per free parameter, burn-in steps, adopted steps are 16, 2000, and
1000, respectively. The error bar of each free parameter is defined as the 16
and 84 percentiles. First, we adopt a single power-law function, where the
free parameters are the power-law index $p$ and the radius $R_{b}$ at a fixed
middle velocity $V_{b}$. $dp=0$ and $V_{\rm sys}=4.0~{}{\rm km~{}s^{-1}}$ are
fixed. The fixed parameter $V_{b}$ is an arbitrary reference velocity where
$R_{b}$ is derived from the fitting. The central stellar mass $M_{*}$ is
calculated from the radius with the fixed inclination angle $i=65\arcdeg$ as
$M_{*}=R(V)(V/\sin i)^{2}/G$, where $R(V)$ is calculated from Equation (1) as
a function of $V$. The calculated $M_{*}$ can thus depend on $V$ if $p\neq
0.5$. The uncertainty of $M_{*}$ is calculated through the error propagation.
The inclination angle is adopted from an outflow model ($60\arcdeg-65\arcdeg$;
Yen et al., 2017) and our continuum analysis ($64.3\arcdeg$ in Section 4.1).
Even with $i=60\arcdeg$, the calculated stellar mass is only 10% higher. The
best-fit parameters are $p=0.5$ and $M_{*}=0.09~{}M_{\sun}$ with the ridge
points, while they are $p=0.67$ and $M_{*}=0.3-0.4~{}M_{\sun}$ with the edge
points. The range of edge $M_{*}$, $0.3-0.4~{}M_{\sun}$ is because $p_{\rm
in}\neq 0.5$. This range is wider than the statistical uncertainty. The
fitting results are summarized in Table 2 including statistical uncertainties.
The uncertainty of the distance adds a relative uncertainty of $\sim 7\%$. The
best-fit functions are plotted in Figure 11. Figure 11(a) shows that the best-
fit functions trace the ridge and edge of the observed PV diagram. Figure
11(b) clearly shows that the ridge and edge indices are close to each other,
as well as to that of Keplerian rotation $(p=0.5)$. The ridge and edge points
were also fitted with a double power-law function. The best-fit parameters are
summarized in Table 2. The double power-law fitting yields similar power-law
indices $p=0.5-0.6$ for the most part of the fitted velocity range, except for
only the lowest velocity channel. These results indicate that the 12CO line
traces the Keplerian disk around I16253.

Figure 11: Edge and ridge points, estimated for each velocity channel,
overlaid on the 12CO major-axis PV diagram (robust$=0.5$) in the (a) linear
and (b) logarithmic scales. The positive radius corresponds to the western
side. The velocity is the relative velocity to the systemic velocity of
$V_{\rm LSR}=4~{}{\rm km~{}s^{-1}}$. The white curves/lines are the best-fit
power-law functions (Section 4.2). The contour levels are in the $3\sigma$
steps from $3\sigma$, where $1\sigma$ is $1.4~{}{\rm mJy~{}beam^{-1}}$.  Table
2: Power-law fitting to the edge and ridge points in the major-axis 12CO PV
diagram.

| Single power
---|---
|  $R_{b}$ (au) |  $V_{b}$ (${\rm km~{}s^{-1}}$) | $p_{\rm in}$ | $dp$ | $M_{*}(M_{\sun})$
Edge | $24.7\pm 0.7$ | 3.18 (fixed) | $0.67\pm 0.04$ | 0 (fixed) | $0.34\pm 0.01-0.43\pm 0.02$
Ridge | $6.2\pm 0.7$ | 3.20 (fixed) | $0.5\pm 0.1$ | 0 (fixed) | $0.09\pm 0.02$
| Double power
|  $R_{b}$ (au) |  $V_{b}$ (${\rm km~{}s^{-1}}$) | $p_{\rm in}$ | $dp$ | $M_{*}(M_{\sun})$
Edge | $45\pm 10$ | $2.2\pm 0.5$ | $0.63\pm 0.07$ | $2\pm 4$ | $0.3\pm 0.2-0.4\pm 0.2$
Ridge | $13\pm 2$ | $2.3\pm 0.2$ | $0.5\pm 0.1$ | $4\pm 4$ | $0.09\pm 0.04$

Note. — The uncertainties of $R_{b}$ and $M_{b}$ in this table are before the
uncertainty of the distance, $\sim 7\%$, is incorpolated.

By identifying the Keplerian rotation, we have estimated the central stellar
mass to be $M_{*}=0.09-0.34~{}M_{\sun}$ directly (kinematically) for the first
time. Even with the lower limit, I16253 already has sufficient mass to evolve
into a low-mass star, beyond the brown-dwarf mass regime
($M_{*}>0.08~{}M_{\sun}$). Maret et al. (2020) demonstrated that the central
stellar masses derived from the ridge and edge points under- and overestimate
the actual stellar mass by $\sim 30\%$ and $\sim 100\%$ in the case where the
beam size is a few tens percent larger than the observed disk size. The reason
why the ridge underestimates the stellar mass is that a part of the gas inside
the Keplerian radius has the same line-of-sight velocity as the gas at the
Keplerian radius, and this inner emission shifts the ridge inward. This effect
is also discussed quantitatively by Aso et al. (2015). The reason why the edge
overestimates the stellar mass is that the observational beam causes emission
outside the outermost, i.e., Keplerian radius, and this outer emission shifts
the edge outward. The condition in Maret et al. (2020) is applicable to our
case (the disk radius of I16253 is discussed in Section 5.1). Then, the
central stellar mass is likely within $M_{*}=0.09\times 1.3$ to
$0.34/2~{}M_{\sun}=0.12$ to $0.17~{}M_{\sun}$. This mass is consistent with a
value derived from the 12CO $J=2-1$ emission at two velocity channels by
assuming Keplerian rotation ($\sim 0.12~{}M_{\sun}$; Hsieh et al., 2019).

###  4.3 Linear velocity gradient in the SO emission

Figure 13 shows the PV diagram in the SO emission along the major axis. The SO
emission is detected within $|V-V_{\rm sys}|\sim\pm 2~{}{\rm km~{}s^{-1}}$,
which is similar to the velocity range where the 13CO and C18O emission is
present (see Figure 14 for the comparison of the SO and C18O emission). The
shape of the SO PV diagram is represented by a linear velocity gradient from
the eastern blueshifted emission to the western redshifted emission in the
velocity range of $|V-V_{\rm sys}|<1.3~{}{\rm km~{}s^{-1}}$, unlike the CO
isotopologues. The SO emission outside this velocity range appears to overlap
a component traced by the CO isotopologue emissions (see Figure 14 for the
comparison), although it is difficult to discuss it given the limited number
of the channels at the common velocity range. The linear velocity gradient can
also be seen in the moment 1 map (Figure 6a), where the high velocities are
located at the outermost parts along the major axis. To evaluate the linear
velocity gradient quantitatively, the ridge points are found at each velocity,
as plotted in Figure 13, and fitted with a linear function of radius $r$
passing $(r,v)=(0,0)$, using emcee with the same condition as that in Section
4.2. The best-fit gradient is estimated to be $g=0.056\pm 0.005~{}{\rm
km~{}s^{-1}}~{}{\rm au}^{-1}$, after the correction for the inclination angle
of $i=65\arcdeg$. This value is slightly smaller than a previous result of
$0.082\pm 0.004~{}{\rm km~{}s^{-1}}~{}{\rm au}^{-1}$ ($d=139$ pc) by Yen et
al. (2017). This difference is mainly due to the different angular
resolutions, $\sim 0\farcs 2$ in the present observations and $\sim 1\arcsec$
in Yen et al. (2017). Their result may thus include SO emission at larger
radii and slower velocities. The SO emission is located within the velocities
of $|V-V_{\rm sys}|<2~{}{\rm km~{}s^{-1}}$, which are lower than those showing
the Keplerian rotation in the 12CO PV diagram (Figure 11), implying that the
SO emission traces a different part from the 12CO emission. The SO velocity
structure and the velocity gradient are discussed in more detail in Section
5.1 along with the infalling and rotating motions traced in the C18O and 13CO
emission.

Figure 12:  Figure 13: Ridge points, found for each velocity, overlaid on the
SO PV diagram (robust$=2$) along the major axis in the linear scale. The
positive radius corresponds to the western side. The velocity is relative to
the systemic velocity of $V_{\rm LSR}=4~{}{\rm km~{}s^{-1}}$. The white line
is the best linear function (Section 4.3). The contour levels are in the
$3\sigma$ steps from $3\sigma$, where $1\sigma$ is $2.3~{}{\rm
mJy~{}beam^{-1}}$.

##  5 Discussion

###  5.1 Disk, Ring, and Envelope

We suggest a picture of the I16253 system based on our results and the
previous results. Hsieh et al. (2019) estimated the specific angular momentum
of the envelope to be $j\sim 45~{}{\rm km~{}s^{-1}}{\rm au}$ ($d=139~{}pc$),
by reproducing the PV diagrams along the major and minor axes in the C18O and
12CO lines with an infalling and rotating envelope model. When the specific
angular momentum and the central stellar mass are given, the velocity field of
a protostellar envelope can be predicted by the UCM envelope model (Ulrich,
1976; Cassen & Moosman, 1981), where the velocity field consists of ballistic,
parabolic flows from an outer boundary in the rigid-body rotation. Figure 14
compares the maximum line-of-sight velocity of the UCM envelope model, as well
as that of the Keplerian disk model, at each position and the observed PV
diagrams in the C18O, SO, and 12CO lines along the major and minor axes. The
model parameters are the specific angular momentum of $j\sim 45~{}{\rm
km~{}s^{-1}}{\rm au}$ and the central stellar mass of $M_{*}\sim
0.14~{}M_{\sun}$ (middle of $0.12-0.17~{}M_{\sun}$ in Section 4.2). With the
specific angular momentum and the central stellar mass, the centrifugal radius
is calculated to be $R_{c}=j^{2}/GM_{*}\sim 16$ au ($0\farcs 12$);
$M_{*}=0.12-0.17~{}M_{\sun}$ corresponds to $R_{c}=13-19$ au. $R_{c}=16$ au is
used as the disk radius to draw the model curves in Figure 14. The C18O PV
diagrams are consistent with the model maximum velocity both in the major and
minor axes in the velocity range lower than the disk velocities, $|V-V_{\rm
sys}|<2~{}{\rm km~{}s^{-1}}$. The obtained $R_{c}$ is close to the radius of
the continuum emission (the red $6\sigma$ contour in Figure 9a), except for
the southeastern extension, supporting that this is the disk radius. This is
also consistent with a relation suggested by Aso & Machida (2020) that the
Gaussian deconvolved radius of the 1.3 mm continuum emission ($\sim 7.4$ au
for I16253; Section 3.1) is $\sim 0.5$ as large as the disk radius for an
evolutionary phase from $M_{*}\sim 0.1$ to $0.4~{}M_{\sun}$; this relation is
based on their synthetic observations of a magnetohydrodynamics simulation of
protostellar evolution.

Figure 14: Comparison between the predicted maximum line-of-sight velocities
(white dashed and dotted curves) and the observed PV diagrams (C18O in red, SO
in green, and 12CO in blue). (a) PV diagram along the major axis. The dashed
curves represent the UCM envelope model with a central stellar mass of
$M_{*}=0.14~{}M_{\sun}$ (Section 4.2) and the specific angular momentum
$j=45~{}{\rm km~{}s^{-1}}~{}{\rm au}$ (Hsieh et al., 2019). The envelope
maximum velocity is calculated only up to $|V-V_{\rm sys}|=2~{}{\rm
km~{}s^{-1}}$. The dotted curves show the maximum line-of-sight velocities of
the Keplerian disk with a disk radius of 16 au. (b) PV diagram along the minor
axis. The dashed curves represent the same UCM envelope as for the major axis.

The SO emission shows a different shape from the envelope and from the disk in
the major-axis PV diagram, whereas the detected velocity range is similar to
that of the C18O (and 13CO) emission as shown in Figure 14. The linear
velocity gradient in the major-axis PV diagram, along with the double peak
structure in the moment 0 map at a higher angular resolution (Figure 6a),
suggests that the SO emission traces a ring close to edge-on due to accretion
shock between the infalling envelope and the disk, as reported in other
protostellar systems (e.g., Yen et al., 2014; Ohashi et al., 2014; Sakai et
al., 2016). The linear velocity gradient in the SO PV diagram $0.056~{}{\rm
km~{}s^{-1}}~{}(r/1~{}{\rm au})$ intersects with the envelope rotation
$45~{}{\rm km~{}s^{-1}}~{}(r/1~{}{\rm au})^{-1}$ at $r=28$ au. This radius is
close to the outermost radius of the SO ridge points (Figure 13). Recent
theoretical studies predict that the accretion shock around the disk occurs at
the radius of $\sim 1.5R_{c}$ (Shariff et al., 2022). The inferred radius of
the SO ring in I16253 is approximately consistent with this prediction. In
conclusion, we suggest that I16253 has the Keplerian disk with $r\sim 16\pm 3$
au traced in the 12CO line, which is surrounded by the shock at $r\sim 28$ au
traced in the SO line due to the mass accretion from the envelope traced in
the C18O and 13CO lines.

###  5.2 Streamer in the SO emission

In addition to the shocked ring, the SO emission also shows an extended
structure to the east as shown in Figure 6(b). Its slightly redshifted
velocity in the outer ($\gtrsim 1\arcsec$) region cannot be explained by the
rotation, the infall motion on the midplane nor the outflow seen in our
observations toward I16253, implying a different motion. A possible
explanation is a streamer as reported in other protostellar systems on
$\lesssim 1000$ au scales (e.g., Yen et al., 2019; Garufi et al., 2022; Thieme
et al., 2022) (see also Kido et al., 2023) as well as on larger scales (e.g.,
Pineda et al., 2020). Hence, we constructed a streamer model by extracting
ballistic, parabolic flows from the UCM envelope model used in Section 5.1.
The free parameters are the two directional angles $(\theta_{0},\phi_{0})$ to
specify the initial polar and azimuthal angle of the streamer trajectory,
respectively. The polar angle $\theta_{0}$ is $0\arcdeg$ at the disk northern
pole and $90\arcdeg$ at the midplane. The azimuthal angle $\phi_{0}$ is
$0\arcdeg$ in the direction of the observer on the midplane and $90\arcdeg$ on
the right seen from the observer. The inclination and position angles of the
midplane are fixed at $65\arcdeg$ and $113\arcdeg$, respectively. By visual
inspection, we found that
$(\theta_{0},\phi_{0})=(60\arcdeg,240\arcdeg-270\arcdeg)$ can reasonably
reproduce the observed SO distribution and velocity. This suggests that the SO
gas is gravitationally bound if it came from these $(\theta_{0},\phi_{0})$
angles, since the UCM envelope model describes how material infalls toward the
central gravitational source. Figure 16 shows the distribution of the model
streamer projected on the plane-of-sky and the line-of-sight velocity of the
streamer model and the rotating ring at a radius of 28 au. The ring rotates in
the same velocity as the UCM envelope model. The model ring appears to
reproduce the size and line-of-sight velocity of the observed compact
component (Figures 6a), while the model streamer appears to reproduce the
extended structure to the east and its line-of-sight velocity (Figure 6b).
Meanwhile, the observed extended structure appears less confined than the
model. This is probably because our model did not consider any beam blurring
effect or radiative transfer effect. More sophisticated modeling including
such effects will help to verify whether the anisotropy in the SO line is
caused by the streamer.

Figure 15:  Figure 16: Model map of the line-of-sight velocity of a streamer
in the UCM envelope ($M_{*}=0.14~{}M_{\sun}$ and $j=45~{}{\rm
km~{}s^{-1}}~{}{\rm au}$ are the same as for Figure 14) and a Keplerian ring
at 28 au. The streamer comes from $30\arcdeg$ above the midplane (i.e.,
$\theta_{0}=60\arcdeg$) and $0\arcdeg-20\arcdeg$ away from the left to the
front seen from the observer around the rotational axis (i.e.,
$\phi_{0}=240\arcdeg-270\arcdeg$).

Ballistic streamers are simulated in a theoretical study of the cloudlet
capture in a protostellar system (Hanawa et al., 2022). The simulation shows
that the streamers can be recognized as enhanced molecular line emission. Such
an enhancement is seen in the eastern peak ($3\sigma$ higher than the western
peak) of the C18O and 13CO moment 0 maps in I16253 (Figure 5), as well as in
the SO extended emission. The eastern extension in the continuum emission
(Figures 1b and 9b) could also be explained by the enhancement due to the
streamer since it is located on the eastern side.

Tobin et al. (2010) show the envelope of I16253 in the $8~{}\micron$
extinction observations at a $\sim 2\arcsec$ resolution using InfraRed Array
Camera (IRAC) on the Spitzer telescope. The extinction is mainly concentrated
in the $\sim 40\arcsec$-long outflow cavity wall, and the northeastern wall
shows higher extinction than the other walls. This may suggest the existence
of inhomogeneities in the ambient material that could produce the anisotropic
streamer observed in our SO result. On the other hand, if the anisotropy in
the SO line is due to an inhomogeneous density structure, similar anisotropy
could also be seen in the C18O and 13CO lines, unlike our results. This may
suggest that the anisotropy in the SO line could be caused by other factors
than density, such as temperature, shock, or chemistry. Figure 18 summarizes
the structures identified around I16253 in our ALMA observations, along with
the large-scale envelope, denoting the radius or length of each structure.

Figure 17:  Figure 18: Schematic picture of the Keplerian disk, ring,
envelope, streamer, outflow identified in the CO isotopologue and SO lines,
and the eastern excess in the 1.3 mm continuum emission in the Class 0
protostar IRAS16253-2429. The gray background represents the large-scale
envelope identified in 8 $\micron$ extinction (Tobin et al., 2010). The radii
or lengths for each of the structures are indicated as well in units of au.

###  5.3 Mass Accretion and Other Quantities in I16253

The central stellar mass of I16253 has been directly estimated to be
$0.12-0.17~{}M_{\sun}$ by identifying the Keplerian rotation in its disk in
this paper, which ranged from $\sim 0.02$ to $\sim 0.12~{}M_{\sun}$ in
previous works (Section 1). In contrast, some previous works estimated the
central stellar mass of a protostar or a proto-BD from the mass accretion rate
from the disk to the central object $\dot{M}_{\rm acc}$. We discuss the
uncertainties for the method using $\dot{M}_{\rm acc}$ along with the
accretion rate derived from the updated stellar mass.

The central stellar mass of I16253 that we have estimated requires a mass
accretion rate of $\dot{M}_{\rm acc}=L_{\rm bol}R_{*}/GM_{*}=(0.9-1.3)\times
10^{-7}~{}M_{\sun}~{}{\rm yr}^{-1}$, where the bolometric luminosity is
$L_{\rm bol}=0.16~{}L_{\sun}$, the stellar radius is assumed to be
$R_{*}=3~{}R_{\sun}$, $G$ is the gravitational constant, and the stellar mass
is $M_{*}=0.12-0.17~{}M_{\sun}$. The stellar radius is predicted in numerical
simulations to be within $\sim 30\%$ of $3~{}R_{\sun}$ (e.g., Masunaga &
Inutsuka, 2000; Vorobyov & Basu, 2015). In comparison, Hsieh et al. (2016)
estimated the mass accretion rate of I16253 from its outflow force. The
outflow force, $F_{\rm out}\sim 7\times 10^{-7}~{}M_{\sun}~{}{\rm
km~{}s^{-1}}~{}{\rm yr}^{-1}$, was measured through IRAM 30-m and APEX
observations in 12CO lines at an angular resolution of $\sim 11\arcsec$. This
force was converted to a mass accretion rate of $\dot{M}_{\rm acc}\sim 5\times
10^{-7}~{}M_{\sun}~{}{\rm yr}^{-1}$ (after the distance correction from 125 to
139 pc) as $\dot{M}_{\rm acc}=F_{\rm out}/\epsilon V_{W}f_{\rm ent}$, where
the ratio between mass loss and accretion rates
$\epsilon=\dot{M}_{W}/\dot{M}_{\rm acc}$ is assumed to be 0.1, the wind (jet)
velocity $V_{W}$ is assumed to be 150 ${\rm km~{}s^{-1}}$, and the entrainment
efficiency $f_{\rm ent}$ is assumed to be 0.1. $F_{\rm out}$ is estimated from
observations in the 12CO $J=2-1$ line using IRAM 30-m and 12CO $J=7-6$ and
$J=6-5$ lines using APEX, with the optical depth correction and assumptions
that the kinetic temperature is 40 K and the H2 number density is
$10^{5}~{}{\rm cm}^{-3}$. This mass accretion rate is $\sim 5$ times higher
than the one derived from the updated stellar mass, whereas being similar to
that of two proto-BD candidates, IC348-SMM2D and L328-IRS, $\sim 2-9\times
10^{-7}~{}M_{\sun}~{}{\rm yr}^{-1}$.

This difference can be explained by large uncertainties in the conversion from
the outflow force to the mass accretion rate. The mass loss and accretion
ratio $\epsilon=\dot{M}_{\rm W}/\dot{M}_{\rm acc}$, the jet velocity $V_{W}$,
and the efficiency $f_{\rm ent}$ can vary from $\sim 0.01$ to $\sim 0.5$
(Ellerbroek et al., 2013; Podio et al., 2021), $\sim 30$ to $\sim 160~{}{\rm
km~{}s^{-1}}$ (jets in SiO $5-4$; Podio et al., 2021), and 0.1 to 0.25 (André
et al., 1999), respectively, in the protostellar phase. A theoretical model
called X-wind suggests that $\epsilon$ and $V_{W}$ are anti-correlated and the
factor of $\epsilon V_{W}$ varies only around $\sim 50$ to $\sim 70~{}{\rm
km~{}s^{-1}}$ (Najita & Shu, 1994) in the T Tauri phase. However, these
velocities are $\gtrsim 4$ times higher than the typical velocity for
protostars, $\sim 15~{}{\rm km~{}s^{-1}}$, as adopted in Hsieh et al. (2016),
and the anti-correlation is not observationally confirmed in the protostellar
phase. Those large uncertainties imply that $\dot{M}_{\rm acc}$ calculated by
Hsieh et al. (2016) could be an order of magnitude lower. The abovementioned
extreme values provide a range of $4\times 10^{-8}$ to $2\times
10^{-5}~{}M_{\sun}~{}{\rm yr}^{-1}$. For this reason, if the central stellar
mass is estimated from the Keplerian rotation, the mass accretion rate could
also be estimated better than from the outflow observations.

The updated mass accretion rate of I16253, $\sim 1\times
10^{-7}~{}M_{\sun}~{}{\rm yr}^{-1}$, cannot provide $M_{*}=0.14~{}M_{\sun}$
within the lifetime of Class 0 $\lesssim$0.2 Myr. In other words, the updated
central stellar mass is large against the luminosity of I16253, as the
required time can be estimated to be $M_{*}/\dot{M}_{\rm
acc}=GM_{*}^{2}/R_{*}L_{\rm bol}\sim 1$ Myr. This suggests that I16253 likely
experienced an accretion burst (or stronger accretion) in the past. The
presence of an accretion burst is supported by the episodic mass ejection as
seen in the 12CO PV diagram along the outflow axis (Section 3.2). The Class 0
age is also long enough to experience bursts with a typical interval of 2400
years in the Class 0 phase (Hsieh et al., 2019). Once an accretion burst
occurs, the luminosity of a protostar increases resulting in higher
temperatures in the surrounding envelope at a given radius. This causes CO
molecules to sublimate from the icy grain mantles in an extended region of the
envelope. After the protostar has returned to its quiescent state and the
luminosity decreased again, the molecules remain in the gas-phase for a period
of time before freezing out again. This freeze-out time-scale depends on the
density at a given radius (Rodgers & Charnley, 2003) but is typically of order
$10^{4}$ years where CO sublimation due to an accretion burst occurs around
Solar-type protostars (e.g., Jørgensen et al., 2015). From this point of view,
the size of the C18O emission observed in I16253 (Figure 5a), $\gtrsim 150$
au, also supports that an accretion burst happened because this size is much
larger than expected from the current luminosity of this protostar: The C18O
emission radius is predicted to be $\sim 30$ au with the current luminosity of
I16253, $0.16~{}L_{\sun}$ (Jørgensen et al., 2015) (see also Lee, 2007).

Based on the measured $M_{*}$ and the above discussion about the mass
accretion and the luminosity, we conclude that I16253 is not a proto-BD even
though its luminosity and mass accretion rate are similar to those of the two
proto-BD candidates. The core (or envelope) mass of I16253 is estimated to be
$\sim 1~{}M_{\sun}$ from 1.1 mm observations with Bolocam on the CSO telescope
at an angular resolution of $31\arcsec$ (Young et al., 2006) and $8~{}\micron$
observations with IRAC on the Spitzer telescope at an angular resolution of
$2\arcsec$ (Tobin et al., 2010). If this envelope mass accretes onto I16253
with a typical star formation efficiency ($30\%\pm 10\%$ calculated from the
dense core mass function and initial mass function; Alves et al., 2007), a
mass of 0.2 to 0.4 $M_{\sun}$ will be added to the central stellar mass. This
also supports an idea that I16253 will ultimately obtain mass high enough to
fuse hydrogen (i.e., $M_{*}>0.08~{}M_{\sun}$).

The Keplerian disk around I16253 has been kinematically identified for the
first time. Without any clear identification of the Keplerian disk, previous
works have attempted to estimate the central stellar mass using the infalling
motion or the outflow force (through the mass accretion rate) in I16253 and
suggested this protostar as a proto-BD candidate. Similarly, the central
stellar mass of proto-BD candidates was estimated in those methods in some
previous works (Section 1), without identifying a Keplerian disk. Our study
demonstrates that a proto-BD and a very low mass protostar must be identified
through the dynamical central stellar mass derived from the Keplerian rotation
of its disk, rather than its infall motion or outflow force.

The accurate mass estimation enables us to compare physical quantities in the
I16253 system with scaling relations among young stellar objects. For example,
its disk radius, $16-19$ au, appears consistent with or slightly lower than
the scaling relation between the disk radius and the central stellar mass
found with the protostellar sample of Yen et al. (2017). Its mass accretion
rate, $(0.9-1.3)\times 10^{-7}~{}M_{\sun}~{}{\rm yr}^{-1}$, and disk mass,
$\gtrsim 2\times 10^{-3}~{}M_{\sun}$, are consistent with the scaling relation
found with the Class I sample of Fiorellino et al. (2022). More observations
around the BD mass threshold with accurate mass estimations will be necessary
to determine whether these scaling relations hold down to the BD mass regime
and bridge star formation and BD formation.

##  6 Conclusions

As a part of the ALMA large project eDisk, we have observed the Class 0
protostar IRAS16253-2429, which has been suggested to be a proto-brown dwarf
candidate in previous works, in the 1.3 mm continuum, 12CO $J=2-1$, C18O
$J=2-1$, 13CO $J=2-1$, SO $J_{N}=6_{5}-5_{4}$, and other molecular lines at an
angular resolution of $0\farcs 07$ ($\sim 10$ au). Our results provide a
typical picture of protostars with a very low stellar mass close to the brown
dwarf threshold (Figure 18). The main results are summarized below.

  1. 1.

The continuum emission shows structures from a $\sim 600$-au scale down to a
$\sim 15$-au scale. The main component shows a disk-like structure with a
radius of $\sim 20$ au. The emission is extended to the southeast along the
major axis. These extensions can be interpreted as an enhancement due to a
streamer from the east and the near-side (southwestern side) wall of the disk-
like structure.

  2. 2.

The 12CO emission overall traces a clear bipolar outflow up to a $\sim
3000$-au scale. The outflow suggests episodic mass ejections. Furthermore, the
12CO emission on the midplane shows a velocity gradient along the disk’s major
axis, implying rotation of the disk.

  3. 3.

We analyzed the 12CO major-axis position-velocity diagram by the edge and
ridge methods and identified a Keplerian disk in IRAS16253-2429 for the first
time. From this identification of the Keplerian rotation in both methods, the
central stellar mass is estimated to be $0.12-0.17~{}M_{\sun}$. This mass
leads us to conclude that IRAS16253-2429 is unlikely a proto-brown dwarf but a
very low mass protostar.

  4. 4.

The C18O and 13CO emissions trace the infalling and rotating envelope as
reported in previous observational works. The major- and minor-axis position-
velocity diagrams in these lines are consistent with the UCM envelope model
with the central stellar mass of $0.14~{}M_{\sun}$ and the disk (centrifugal)
radius of $r\sim 16$ au. Their moment 0 maps exhibit stronger emission
intensities on the eastern peak, which could result from the streamer
mentioned below.

  5. 5.

The SO emission shows a different velocity structure from the CO
isotopologues. Its double peaks in the moment 0 map and linear velocity
gradient in the major-axis position-velocity diagram suggest that this
emission traces a ring ($r\sim 28$ au) due to the accretion shock between the
disk and the envelope. This emission also shows a streamer from the eastern
side, which can be explained with a ballistic, parabolic flow at $30\arcdeg$
above the midplane extracted from the UCM envelope model.

## Acknowledgements

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2019.1.00261.L
and ADS/JAO.ALMA#2019.A.00034.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. The National Radio Astronomy Observatory is a facility of the National
Science Foundation operated under cooperative agreement by Associated
Universities, Inc. W.K. was supported by the National Research Foundation of
Korea (NRF) grant funded by the Korea government (MSIT)
(NRF-2021R1F1A1061794). N.O. acknowledges support from National Science and
Technology Council (NSTC) in Taiwan through the grants NSTC 109-2112-M-001-051
and 110-2112-M-001-031. J.K.J. acknowledge support from the Independent
Research Fund Denmark (grant No. 0135-00123B). J.J.T. acknowledges support
from NASA XRP 80NSSC22K1159. N.O. acknowledges support from National Science
and Technology Council (NSTC) in Taiwan through the grants NSTC
109-2112-M-001-051 and 110-2112-M-001-031. Y.A. acknowledges support by NAOJ
ALMA Scientific Research Grant code 2019-13B, Grant-in-Aid for Scientific
Research (S) 18H05222, and Grant-in-Aid for Transformative Research Areas (A)
20H05844 and 20H05847. IdG acknowledges support from grant
PID2020-114461GB-I00, funded by MCIN/AEI/10.13039/501100011033. PMK
acknowledges support from NSTC 108-2112- M-001-012, NSTC 109-2112-M-001-022
and NSTC 110-2112-M-001-057. SPL and TJT acknowledge grants from the National
Science and Technology Council of Taiwan 106-2119-M-007-021-MY3 and
109-2112-M-007-010-MY3. J.K.J. acknowledges support from the Independent
Research Fund Denmark (grant No. 0135-00123B). C.W.L. is supported by the
Basic Science Research Program through the National Research Foundation of
Korea (NRF) funded by the Ministry of Education, Science and Technology (NRF-
2019R1A2C1010851), and by the Korea Astronomy and Space Science Institute
grant funded by the Korea government (MSIT; Project No. 2022-1-840-05). JEL
was supported by the National Research Foundation of Korea (NRF) grant funded
by the Korean government (MSIT) (grant number 2021R1A2C1011718). ZYL is
supported in part by NASA 80NSSC20K0533 and NSF AST-1910106. ZYDL acknowledges
support from NASA 80NSSC18K1095, the Jefferson Scholars Foundation, the NRAO
ALMA Student Observing Support (SOS) SOSPA8-003, the Achievements Rewards for
College Scientists (ARCS) Foundation Washington Chapter, the Virginia Space
Grant Consortium (VSGC), and UVA research computing (RIVANNA). LWL
acknowledges support from NSF AST-2108794. S.N. acknowledges support from the
National Science Foundation through the Graduate Research Fellowship Program
under Grant No. 2236415. Any opinions, findings, and conclusions or
recommendations expressed in this material are those of the authors and do not
necessarily reflect the views of the National Science Foundation. R.S.
acknowledge support from the Independent Research Fund Denmark (grant No.
0135-00123B). S.T. is supported by JSPS KAKENHI grant Nos. 21H00048 and
21H04495, and by NAOJ ALMA Scientific Research grant No. 2022-20A. JPW
acknowledges support from NSF AST-2107841. H.-W.Y. acknowledges support from
the National Science and Technology Council (NSTC) in Taiwan through the grant
NSTC 110-2628-M-001-003-MY3 and from the Academia Sinica Career Development
Award (AS-CDA-111-M03).

##  Appendix A Other molecular lines

Figure 19 shows the moment 0 and 1 maps of the lines observed in the eDisk
project toward I16253, except for the 12CO, C18O, 13CO, and SO lines. These
lines were observed from the wide spectral windows otherwise intended for
continuum measurements. The three H2CO lines have different upper state
energies, depending on the rest frequencies (10.5 K; 218.22 GHz, 57.6 K;
218.47 GHz, 57.6 K; 218.76 GHz). The three cyclic C3H2 lines also have
different upper state energies (28.2 K; 217.82 GHz, 25.0 K; 217.94 GHz, 24.9
K; 218.16 GHz). The CH3OH, DCN, and SiO lines have upper state energies of
35.0, 10.4 and 20.8 K, respectively. In comparison, the 12CO, C18O, 13CO, and
SO lines have upper state energies of 5.5, 5.3, 5.3, and 24.4 K, respectively.
All the panels in Figure 19 are made from three channels centered at the
systemic velocity ($V_{\rm LSR}=2.66-5.34~{}{\rm km~{}s^{-1}}$) and show the
same spatial and velocity ranges.

Figure 19: Moment 0 and 1 maps in 5 other molecular lines (9 transitions)
observed in the ALMA observations of the eDisk project. These images are made
with the robust parameter of 2. The emission is integrated from 2.66 to
$5.34~{}{\rm km~{}s^{-1}}$ (3 channels) for all lines. The contour levels are
in $3\sigma$ steps from $3\sigma$, where $1\sigma$ is $1.5~{}{\rm
mJy~{}beam^{-1}}~{}{\rm km~{}s^{-1}}$. The ellipse at the lower left corner in
each panel denotes the synthesized beam, $\sim 0\farcs 38\times 0\farcs 27\
(78\arcdeg)$. The diagonal lines are the major and minor axes of the continuum
emission, P.A.$=113\arcdeg$ and $23\arcdeg$.

## References

  * Alves et al. (2017) Alves, F. O., Girart, J. M., Caselli, P., et al. 2017, A&A, 603, L3. doi:10.1051/0004-6361/201731077
  * Alves et al. (2007) Alves, J., Lombardi, M., & Lada, C. J. 2007, A&A, 462, L17. doi:10.1051/0004-6361:20066389
  * André et al. (1999) André, P., Motte, F., & Bacmann, A. 1999, ApJ, 513, L57. doi:10.1086/311908
  * Andrews & Williams (2005) Andrews, S. M. & Williams, J. P. 2005, ApJ, 631, 1134. doi:10.1086/432712
  * Aso & Machida (2020) Aso, Y. & Machida, M. N. 2020, ApJ, 905, 174. doi:10.3847/1538-4357/abc6fc
  * Aso & Sai (2023) Aso, Y. & Sai, J. 2023, jinshisai/SLAM: First Release of SLAM, v1.0.0, Zenodo, doi:10.5281/zenodo.7783868
  * Aso et al. (2017) Aso, Y., Ohashi, N., Aikawa, Y., et al. 2017, ApJ, 849, 56. doi:10.3847/1538-4357/aa8264
  * Aso et al. (2015) Aso, Y., Ohashi, N., Saigo, K., et al. 2015, ApJ, 812, 27. doi:10.1088/0004-637X/812/1/27
  * Aso et al. (2021) Aso, Y., Kwon, W., Hirano, N., et al. 2021, ApJ, 920, 71. doi:10.3847/1538-4357/ac15f3
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123. doi:10.3847/1538-3881/aabc4f
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33. doi:10.1051/0004-6361/201322068
  * Basu et al. (2012) Basu, S., Vorobyov, E. I., & DeSouza, A. L. 2012, First Stars IV - from Hayashi to the Future -, 1480, 63. doi:10.1063/1.4754329
  * Bate (2012) Bate, M. R. 2012, MNRAS, 419, 3115. doi:10.1111/j.1365-2966.2011.19955.x
  * Beckwith et al. (1990) Beckwith, S. V. W., Sargent, A. I., Chini, R. S., et al. 1990, AJ, 99, 924. doi:10.1086/115385
  * Cassen & Moosman (1981) Cassen, P. & Moosman, A. 1981, Icarus, 48, 353. doi:10.1016/0019-1035(81)90051-8
  * Chabrier (2002) Chabrier, G. 2002, ApJ, 567, 304. doi:10.1086/324716
  * Cabrit (1989) Cabrit, S. 1989, European Southern Observatory Conference and Workshop Proceedings, 33, 119
  * de Geus et al. (1989) de Geus, E. J., de Zeeuw, P. T., & Lub, J. 1989, A&A, 216, 44
  * de Gregorio-Monsalvo et al. (2016) de Gregorio-Monsalvo, I., Barrado, D., Bouy, H., et al. 2016, A&A, 590, A79. doi:10.1051/0004-6361/201424149
  * Dunham et al. (2008) Dunham, M. M., Crapsi, A., Evans, N. J., et al. 2008, ApJS, 179, 249. doi:10.1086/591085
  * Ellerbroek et al. (2013) Ellerbroek, L. E., Podio, L., Kaper, L., et al. 2013, A&A, 551, A5. doi:10.1051/0004-6361/201220635
  * Fiorellino et al. (2022) Fiorellino, E., Tychoniec, Ł., Manara, C. F., et al. 2022, ApJ, 937, L9. doi:10.3847/2041-8213/ac8fee
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., et al. 2013, PASP, 125, 306. doi:10.1086/670067
  * Garufi et al. (2022) Garufi, A., Podio, L., Codella, C., et al. 2022, A&A, 658, A104. doi:10.1051/0004-6361/202141264
  * Hanawa et al. (2022) Hanawa, T., Sakai, N., & Yamamoto, S. 2022, ApJ, 932, 122. doi:10.3847/1538-4357/ac6e6a
  * Hsieh et al. (2019) Hsieh, T.-H., Hirano, N., Belloche, A., et al. 2019, ApJ, 871, 100. doi:10.3847/1538-4357/aaf4fe
  * Hsieh et al. (2019) Hsieh, T.-H., Murillo, N. M., Belloche, A., et al. 2019, ApJ, 884, 149. doi:10.3847/1538-4357/ab425a
  * Hsieh et al. (2016) Hsieh, T.-H., Lai, S.-P., Belloche, A., et al. 2016, ApJ, 826, 68. doi:10.3847/0004-637X/826/1/68
  * Huélamo et al. (2017) Huélamo, N., de Gregorio-Monsalvo, I., Palau, A., et al. 2017, A&A, 597, A17. doi:10.1051/0004-6361/201628510
  * Jørgensen et al. (2015) Jørgensen, J. K., Visser, R., Williams, J. P., et al. 2015, A&A, 579, A23. doi:10.1051/0004-6361/201425317
  * Kido et al. (2023) Kido, M., Takakuwa, S., Saigo, K., et al. 2023, ApJ, 953, 190. doi:10.3847/1538-4357/acdd7a
  * Kim et al. (2019) Kim, G., Lee, C. W., Maheswar, G., et al. 2019, ApJS, 240, 18. doi:10.3847/1538-4365/aaf889
  * Knude & Hog (1998) Knude, J. & Hog, E. 1998, A&A, 338, 897
  * Lee (2020) Lee, C.-F. 2020, A&A Rev., 28, 1. doi:10.1007/s00159-020-0123-7
  * Lee et al. (2018) Lee, C. W., Kim, G., Myers, P. C., et al. 2018, ApJ, 865, 131. doi:10.3847/1538-4357/aadcf6
  * Lee (2007) Lee, J.-E. 2007, Journal of Korean Astronomical Society, 40, 83. doi:10.5303/JKAS.2007.40.4.083
  * Machida et al. (2009) Machida, M. N., Inutsuka, S.-. ichiro ., & Matsumoto, T. 2009, ApJ, 699, L157. doi:10.1088/0004-637X/699/2/L157
  * Machida et al. (2010) Machida, M. N., Inutsuka, S.-. ichiro ., & Matsumoto, T. 2010, ApJ, 724, 1006. doi:10.1088/0004-637X/724/2/1006
  * Machida et al. (2014) Machida, M. N., Inutsuka, S.-. ichiro ., & Matsumoto, T. 2014, MNRAS, 438, 2278. doi:10.1093/mnras/stt2343
  * Maret et al. (2020) Maret, S., Maury, A. J., Belloche, A., et al. 2020, A&A, 635, A15. doi:10.1051/0004-6361/201936798
  * Masunaga & Inutsuka (2000) Masunaga, H. & Inutsuka, S.-. ichiro . 2000, ApJ, 531, 350. doi:10.1086/308439
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., et al. 2007, Astronomical Data Analysis Software and Systems XVI, 376, 127
  * Najita & Shu (1994) Najita, J. R. & Shu, F. H. 1994, ApJ, 429, 808. doi:10.1086/174365
  * Ohashi et al. (2014) Ohashi, N., Saigo, K., Aso, Y., et al. 2014, ApJ, 796, 131. doi:10.1088/0004-637X/796/2/131
  * Ohashi et al. (2023) Ohashi, N., Tobin, J. J., Jørgensen, J. K., et al. 2023, ApJ, 951, 8. doi:10.3847/1538-4357/acd384
  * Padoan & Nordlund (2004) Padoan, P. & Nordlund, Å. 2004, ApJ, 617, 559. doi:10.1086/345413
  * Palau et al. (2014) Palau, A., Zapata, L. A., Rodríguez, L. F., et al. 2014, MNRAS, 444, 833. doi:10.1093/mnras/stu1461
  * Park et al. (2021) Park, W., Lee, J.-E., Contreras Peña, C., et al. 2021, ApJ, 920, 132. doi:10.3847/1538-4357/ac1745
  * Pineda et al. (2020) Pineda, J. E., Segura-Cox, D., Caselli, P., et al. 2020, Nature Astronomy, 4, 1158. doi:10.1038/s41550-020-1150-z
  * Podio et al. (2021) Podio, L., Tabone, B., Codella, C., et al. 2021, A&A, 648, A45. doi:10.1051/0004-6361/202038429
  * Rodgers & Charnley (2003) Rodgers, S. D. & Charnley, S. B. 2003, ApJ, 585, 355. doi:10.1086/345497
  * Ruíz-Rodríguez et al. (2022) Ruíz-Rodríguez, D. A., Williams, J. P., Kastner, J. H., et al. 2022, MNRAS, 515, 2646. doi:10.1093/mnras/stac1879
  * Sai et al. (2022) Sai, J., Ohashi, N., Maury, A. J., et al. 2022, ApJ, 925, 12. doi:10.3847/1538-4357/ac341d
  * Sai et al. (2020) Sai, J., Ohashi, N., Saigo, K., et al. 2020, ApJ, 893, 51. doi:10.3847/1538-4357/ab8065
  * Sakai et al. (2016) Sakai, N., Oya, Y., López-Sepulcre, A., et al. 2016, ApJ, 820, L34. doi:10.3847/2041-8205/820/2/L34
  * Santamaría-Miranda et al. (2021) Santamaría-Miranda, A., de Gregorio-Monsalvo, I., Plunkett, A. L., et al. 2021, A&A, 646, A10. doi:10.1051/0004-6361/202039419
  * Seifried et al. (2016) Seifried, D., Sánchez-Monge, Á., Walch, S., et al. 2016, MNRAS, 459, 1892. doi:10.1093/mnras/stw785
  * Shariff et al. (2022) Shariff, K., Gorti, U., & Melon Fuksman, J. D. 2022, MNRAS, 514, 5548. doi:10.1093/mnras/stac1186
  * Stamatellos & Whitworth (2009) Stamatellos, D. & Whitworth, A. P. 2009, 15th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, 1094, 557. doi:10.1063/1.3099172
  * Thieme et al. (2022) Thieme, T. J., Lai, S.-P., Lin, S.-J., et al. 2022, ApJ, 925, 32. doi:10.3847/1538-4357/ac382b
  * Tobin et al. (2023) Tobin, J. 2023, eDisk data reduction scripts, 1.0.0, Zenodo, doi:10.5281/zenodo.7986682
  * Tobin et al. (2012) Tobin, J. J., Hartmann, L., Bergin, E., et al. 2012, ApJ, 748, 16. doi:10.1088/0004-637X/748/1/16
  * Tobin et al. (2010) Tobin, J. J., Hartmann, L., Looney, L. W., et al. 2010, ApJ, 712, 1010. doi:10.1088/0004-637X/712/2/1010
  * Tobin et al. (2020) Tobin, J. J., Sheehan, P. D., Megeath, S. T., et al. 2020, ApJ, 890, 130. doi:10.3847/1538-4357/ab6f64
  * Ulrich (1976) Ulrich, R. K. 1976, ApJ, 210, 377. doi:10.1086/154840
  * Vorobyov & Basu (2015) Vorobyov, E. I. & Basu, S. 2015, ApJ, 805, 115. doi:10.1088/0004-637X/805/2/115
  * Vousden et al. (2016) Vousden, W. D., Farr, W. M., & Mandel, I. 2016, MNRAS, 455, 1919. doi:10.1093/mnras/stv2422
  * Whitworth & Zinnecker (2004) Whitworth, A. P. & Zinnecker, H. 2004, A&A, 427, 299. doi:10.1051/0004-6361:20041131
  * Wilner & Welch (1994) Wilner, D. J., & Welch, W. J. 1994, ApJ, 427, 898
  * Yen et al. (2019) Yen, H.-W., Gu, P.-G., Hirano, N., et al. 2019, ApJ, 880, 69. doi:10.3847/1538-4357/ab29f8
  * Yen et al. (2017) Yen, H.-W., Koch, P. M., Takakuwa, S., et al. 2017, ApJ, 834, 178. doi:10.3847/1538-4357/834/2/178
  * Yen et al. (2014) Yen, H.-W., Takakuwa, S., Ohashi, N., et al. 2014, ApJ, 793, 1. doi:10.1088/0004-637X/793/1/1
  * Yen et al. (2013) Yen, H.-W., Takakuwa, S., Ohashi, N., et al. 2013, ApJ, 772, 22. doi:10.1088/0004-637X/772/1/22
  * Young et al. (2006) Young, K. E., Enoch, M. L., Evans, N. J., et al. 2006, ApJ, 644, 326. doi:10.1086/503327
  * Zucker et al. (2020) Zucker, C., Speagle, J. S., Schlafly, E. F., et al. 2020, A&A, 633, A51. doi:10.1051/0004-6361/201936145
  * Zucker et al. (2019) Zucker, C., Speagle, J. S., Schlafly, E. F., et al. 2019, ApJ, 879, 125. doi:10.3847/1538-4357/ab2388