Early Planet Formation in Embedded Disks (eDisk) VI:
Kinematic Structures around the Very Low Mass Protostar 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 N₂H⁺ $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 C¹⁸O $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 C¹⁸O 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 ¹²CO and C¹⁸O $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 ¹²CO emission arises from a Keplerian disk; they found offsets of the ¹²CO 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 ¹²CO, ¹³CO, C¹⁸O, 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.

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 2(a) shows the integrated intensity (moment 0) and the mean velocity (moment 1) maps in the ¹²CO 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 ¹²CO 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 ¹²CO 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(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 ¹²CO 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.

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.

Figure 5 shows the moment 0 and 1 maps in the C¹⁸O and ¹³CO 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 C¹⁸O line than in the ¹³CO line with both robust parameters. This is because the ¹³CO emission tends to be optically thicker and more extended than the C¹⁸O emission, and thus the ¹³CO emission is resolved out more severely around the systemic velocity.

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 C¹⁸O and ¹³CO lines. On the other hand, the mean velocity of the SO emission is larger outside the peaks than inside the peaks, unlike the C¹⁸O and ¹³CO 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 C¹⁸O and ¹³CO 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 C¹⁸O and ¹³CO 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 C¹⁸O and ¹³CO lines and also different from the velocity in the eastern cavity wall of the ¹²CO outflow (Figure 2a). This SO component is discussed in more detail in Section 5.2.

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 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 C¹⁸O and ¹³CO moment 0 maps.

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 ¹²CO (blue contours), C¹⁸O (red contours), and ¹³CO (color) lines. The C¹⁸O and ¹³CO 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 ¹²CO 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 C¹⁸O and ¹³CO 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 ¹²CO 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 ¹²CO line traces a purely rotating, i.e., Keplerian disk, while the C¹⁸O and ¹³CO 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 ¹²CO 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 ¹³CO and C¹⁸O emission is not detected sufficiently in this compact region at the high velocities where the ¹²CO emission is detected. The nondetection in the ¹³CO and C¹⁸O lines suggests the ¹²CO emission is not optically thick at the high velocities. In contrast to the high velocities, the ¹²CO emission is more extended and strongly resolved out at the low velocities where the ¹³CO and C¹⁸O emission traces the envelope.

In order to verify whether the ¹²CO 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)²²2https://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$:

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 ¹²CO line traces the Keplerian disk around I16253.

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 ¹²CO $J=2-1$ emission at two velocity channels by assuming Keplerian rotation ($\sim 0.12~{}M_{\sun}$; Hsieh et al., 2019).

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 ¹³CO and C¹⁸O emission is present (see Figure 14 for the comparison of the SO and C¹⁸O 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 ¹²CO PV diagram (Figure 11), implying that the SO emission traces a different part from the ¹²CO 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 C¹⁸O and ¹³CO emission.

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 C¹⁸O and ¹²CO 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 C¹⁸O, SO, and ¹²CO 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 C¹⁸O 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.

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 C¹⁸O (and ¹³CO) 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 ¹²CO 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 C¹⁸O and ¹³CO lines.

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.

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 C¹⁸O and ¹³CO 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 C¹⁸O and ¹³CO 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.

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 ¹²CO 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 ¹²CO $J=2-1$ line using IRAM 30-m and ¹²CO $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 H₂ 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 ¹²CO 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 C¹⁸O 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 C¹⁸O 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.