Early Planet Formation in Embedded Disks (eDisk). XII. Accretion Streamers, Protoplanetary Disk, and Outflow in the Class I Source Oph IRS 63

The 1.3 mm continuum image of the disk around Oph IRS 63 is shown in Figure 1 (left panel). We see a disk that is centrally peaked with a peak intensity of 5.8 mJy beam⁻¹ or 170 K. We fit a 2D Gaussian profile using the imfit function in CASA to measure the basic dust disk properties (see also Ohashi et al. 2023). We obtained a disk position angle (P.A.) = 149° (measured from north through east), a central position at coordinates R.A. = 16: 31: 35.654, decl. = −24: 01: 30.08, and a deconvolved disk aspect ratio of 0.685 that we transform to disk inclination of i_disk = 467. Because the Gaussian fit does not properly account for the extended dust continuum structure, we integrate the flux above 5σ to derive a total flux density of 317 ± 0.5 mJy.

Zoom In Zoom Out Reset image size

To measure the dust disk size, we adopted the now-common definition of the radius enclosing a certain fraction of the total flux (e.g., Sanchis et al. 2021). We used GoFish (Teague 2019) to compute an azimuthally averaged flux (assuming inclination of i = 47° and P.A. = 149°) and defined the radii enclosing 68%, 90%, and 95% of total flux as R_68% = 017, R_90% = 038, and R_95% = 047, respectively. The difference in disk size at different flux cuts is due to the strongly peaked central emission observed for Oph IRS 63.

From Figure 1, we see that beyond the central peak there is at least one diffuse annular structure. Such annular structures around Oph IRS 63 were studied in detail by Segura-Cox et al. (2020), where two ring-like structures were identified at R ∼ 019 and R ∼ 035 with widths of 004 and 009, respectively. In the right panel of Figure 1, we show azimuthally averaged intensity profiles of the continuum emission for different robust parameters. We visually confirm the presence of an inner ring at ∼016 and an outer ring at 034. In this work, we do not intend to perform any detailed analysis of the continuum emission, as a future paper will address the substructures and asymmetries observed in all the sources of the eDisk sample.

We calculate the mass of the dust disk following Equation (1) proposed by Hildebrand (1983), which assumes that the emission is optically thin and comes only from a thermal component. In this equation, S_ν is the flux density calculated by integrating the emission above 5σ, D is the distance to the source of 132 pc, B_ν (T_d ) is the Planck function at the observed frequency (225 GHz), T_d is the dust temperature, and κ_ν is the absorption coefficient adopted as κ_ν = 2.3 cm² per gram of dust (Beckwith et al. 1990):

$$\begin{eqnarray}&&M=\frac{{S}_{\nu }{D}^{2}}{{\kappa }_{\nu }{B}_{\nu }({T}_{d})}.\end{eqnarray} \tag{ 1 }$$

We use two different dust temperatures to measure the mass dust of the disk: (1) the standard T_dust = 20 K value often used in Class II sources (e.g., Pascucci et al. 2016), from which we obtain a total dust disk mass of M_dust ∼ 5 × 10⁻⁴ M_⊙ (or M_dust ∼ 160 M_⊕); (2) the temperature prescription proposed by Tobin et al. (2020) of ${T}_{\mathrm{dust}}=43\,{({L}_{\mathrm{bol}})}^{1/4}$, which for Oph IRS 63 yields a temperature T_dust = 46 K, and therefore a dust mass of M_dust ∼ 2 × 10⁻⁴ M_⊙ (or M_dust ∼ 60 M_⊕). These values should be regarded as lower limits for the total disk mass of solids since the 1.3 mm emission in the disk is optically thick as discussed in Segura-Cox et al. (2020).

In Figure 2, we present a gallery of the detected molecular tracers described in this work. The images correspond to integrated intensity maps (moment 0) of three CO isotopologues, ¹²CO, ¹³CO, and C¹⁸O, in addition to the SO and H₂CO lines. The different tracers reveal a large-scale envelope, a compact disk, spiral-like features, and an outflow. Although some lines trace multiple structures, below, we attempt to present the different components of the protostellar system, focusing on the most important tracers for each case.

Zoom In Zoom Out Reset image size

In some molecular tracers, our observations suffer from cloud contamination at channel maps close to systemic velocity. This feature is most likely caused by gas being filtered out by the ALMA interferometer. Additionally, we see a strong negative absorption at the disk position for all the molecular species (except in SO). This artifact is caused by continuum oversubtraction, resulting from optically thick molecular lines at the line center, and it must be kept in mind when discussing some molecular maps close to the stellar position (see channel maps in Appendix B).

Following the results of the dust continuum, we now focus on the gas emission at or close to the dusty disk position. Figures 2 and 3 show three prominent tracers of gas at the disk position: ¹³CO, C¹⁸O, and SO. The H₂CO line only shows weak emission, which is strongly affected by continuum oversubtraction.

Zoom In Zoom Out Reset image size

Figure 3 shows zoom-in maps of the molecular tracers centered around the disk position. The moment 0 maps of ¹³CO and C¹⁸O show emission mostly elongated along the dust disk major axis, and their moment 1 maps show a clear velocity gradient along this direction. The SO and H₂CO moment 0 maps (shown in Figure 3) show a velocity gradient along the direction of the disk major axis, but they do not show an elongation on the moment 0 maps. Instead, they show a complex pattern with an apparent lack of emission at the center of the map. Although the lack of central emission in H₂CO can be attributed to continuum oversubtraction (see Figure B5 of Appendix B), the ring-like pattern in SO may be real and is further described in Section 3.5.

The velocity gradients along the major axes and the elongated emission (in the same direction as the dust disk major axis) in ¹³CO and C¹⁸O suggest the presence of a rotating disk-like structure around Oph IRS 63. The size of the gas disk and the exact velocity pattern of the rotation, e.g., if it is consistent with Keplerian rotation or not, are further analyzed in Section 4.1.

In these close-in maps, we also observe an S-shaped pattern at the transitional zone between blue- and redshifted emission in C¹⁸O and SO. This S shape depicts a change in the direction of the velocity gradient between an outer region (≳1'') and the inner region (≲1'') that aligns with the dust disk's P.A. This change in velocity gradient is a clear example of an infalling envelope connecting with a rotating disk structure (as shown previously in, e.g., Momose et al. 1998; Murillo et al. 2013; Yen et al. 2013).

From the gallery plot of Figure 2, we see that a large envelope structure is detected in ¹³CO and C¹⁸O, which extends to distances greater than 10'' (∼1300 au). In fact, from Figure 4, we see ¹³CO extended asymmetric emission reaching ∼13'' from the stellar position. While most of the extended emission is relatively faint, it is still detected at levels of ∼22 mJy beam⁻¹ (S/N ∼ 5). The large-scale structure of the ¹³CO is highly complex, with several arcs/arms features extending to the northwest, northeast, and southeast. In contrast, toward the southwest direction, there is a lack of extended structure, which coincides with the position where the ¹²CO blueshifted outflow is observed. The large-scale velocity structure of ¹³CO (right panel in Figure 4) shows a velocity gradient mainly in the north–south direction, with most of the redshifted emission in the north and the blueshifted emission in the south.

Zoom In Zoom Out Reset image size

Figure 5 shows the moment 8 (maximum value in the spectrum) and moment 1 maps for the C¹⁸O line. Similar to the ¹³CO line, we see extended emission reaching ∼13'', but in the case of C¹⁸O, the extended emission is mostly detected toward the northern part of the system. In fact, a large arc-like feature is evident from the moment 8 map, which branches off to the northeast and northwest. In addition, a small arc structure (∼5'' in length) is detected in the eastern part at levels of ∼20–40 mJy beam⁻¹.

Zoom In Zoom Out Reset image size

In both CO isotopologues, the northern emission (redshifted) is brighter and more extended, consistent with our interpretation that the northern part is the near side of the system and the southern part is partially hidden behind the large-scale ¹²CO outflow (see Section 5.6). The observed north–south brightness asymmetry could also be caused by the large number of channels missed owing to cloud contamination, which affects more strongly the ¹³CO as more emission is being resolved out.

Finally, when considering the orientation of the ¹³CO and C¹⁸O envelopes, we find that they have velocity gradient directions roughly aligned north–south, as measured from the peak emission in moment 1 and 8 maps. The orientation of these large-scale velocity gradients, however, differs from what is observed at smaller scales (≲1''), where the velocity gradient aligns with the P.A. of the dust disk (i.e., P.A. = 149°, as shown in Figure 3). Similar to the S-shape zero-velocity line in SO, the change in velocity gradient direction between large and small scales in ¹³CO and C¹⁸O is indicative of an infalling envelope (e.g., Sakai et al. 2016).

One of the most striking observational results for Oph IRS 63 is the detection of a spiral-like structure in the SO 6₅ − 5₄ tracer. Figure 6 shows the SO integrated intensity map, with annotations marking the spiral-like structures that connect at the north (S1), east (S2), and south (S3) with the disk. The apparent motion of the three spirals is as if they were rotating clockwise. The S1 and S3 spirals are tightly wound around the disk, while S2 extends outward toward the south. The three spirals connect to a ring-like structure in the inner 05 of the SO emission. The S1 spiral connects with the annular structure at a P.A. of ∼320°, where the gas has its maximum intensity at about 200 mJy beam⁻¹ km s⁻¹. The S2 spiral connects with the ring at a P.A. of ∼90° (almost exactly east), but it does not produce an emission peak. The S3 spiral connects with the ring at a P.A. of ∼170°, where the gas emission reaches about 190 mJy beam⁻¹ km s⁻¹.

Zoom In Zoom Out Reset image size

Figure 3 shows the SO velocity map, which exhibits a dominant velocity gradient along the major axis of the dust disk, similar to that found for ¹³CO and C¹⁸O at distances ≲1'', but it becomes significantly distorted with increasing radius (leading to the S-shape zero-velocity line). First, the S3 spiral coincides with the blueshifted velocity pattern, while the S1 spiral traces an opposite structure with a redshifted velocity. The S2 spiral starts close to the disk at a positive but almost zero-velocity offset and then becomes blueshifted at increasing radial distances. From the channel maps presented in Appendix B (Figure B4), we can see that most of the spiral structures appear close to the systemic velocity at values of ± 1.3 km s⁻¹ with respect to V_sys = 2.8 km s⁻¹. These velocity ranges are similar to the extended envelope structure identified in ¹³CO and C¹⁸O and significantly different from the outflow detected in ¹²CO.

As briefly mentioned above, an annular structure with a radius of ∼05 ± 02 is observed on the moment 0 map, and this ring-like feature peaks close to the edge of the dust disk (see Section 3.1). The center of the SO ring, or more precisely, the center of the emission drop, is offset by approximately 011 relative to the dust continuum emission, as measured through a 2D Gaussian fitting (see Figure 6). Although small, this difference is significant given that we can measure the positions of the SO emission with an accuracy of 0023. ¹⁹

In the literature, different SO transitions can trace various components of a protostar, including jets, low-velocity outflow, warm inner envelope, and disk emission (e.g., Tychoniec et al. 2021). Due to the kinematics and spatial distribution of SO shown, we can rule out an outflow origin for the case of Oph IRS 63. Instead, given the larger rotating envelope structure detected in ¹³CO and the prominent S-shape pattern at the transitional zone between blue- and redshifted emission, we speculate that it comes from the connection between the warm inner envelope and the disk. This idea is further discussed in Section 5.1.

A different type of arc or spiral-like structures are found in the C¹⁸O moment maps (shown in Figure 5). These are large-scale features (>5''), have lower pitch angles than the SO (less tightly wound), and can definitely be associated with the envelope structure. Three to four of these structures can be identified from the moment maps of C¹⁸O (see Figure 5), but the most conspicuous ones are the northern and eastern ones extending ∼13'' and ∼5'', respectively. Section 4 explores the idea that these structures result from an infalling rotating envelope that feeds the disk with material from the collapsing cloud.

In Figure 7, we plot the ¹²CO emission around Oph IRS 63, revealing a large-scale bipolar outflow structure. The outflow extends to the northeast and southwest direction until the edge of the map, where the sensitivity steadily decreases owing to the poorer primary beam response. Panels (a) and (b) of Figure 7 display a strong and blueshifted outflow emission toward the southwest. Conversely, the northeast outflow is redshifted and weaker, appearing only in a few channel maps. A tentative detection of the inner ∼1'' of this outflow was recently reported in Antilen et al. (2023).

Zoom In Zoom Out Reset image size

The blueshifted side of the outflow exhibits shell-like structures, which can be roughly categorized into narrow and broad components. The narrow component is clearly detected up to velocities of −15 km s⁻¹, and it forms a closed loop-like structure with a maximum radial size of ∼8'' at velocities of −3.0 to −7.0 km s⁻¹. In the moment 0 map (see Figure 2), the narrow blueshifted component displays a characteristic V-shape pattern. We estimated the outflow's P.A. to be ∼240° by averaging the P.A.s of the two legs of the V-shape structure.

The broad component is detected at approximate velocities of −3 to −9.0 km s⁻¹. The latter always fully encompasses the narrower component, and its opening angle decreases with increasing negative velocity. A similar configuration, in which a broader outflow fully contains a narrower component, has been observed in ALMA observations of DO Tau (Fernández-López et al. 2020), DK Cha (Harada et al. 2023), and GSS 30 IRS3 (A. Santamaria-Miranda & eDisk team 2023, in preparation). Near-infrared observations of such a configuration have also been reported, for example, for HL Tau (Takami et al. 2007). Furthermore, the narrower blueshifted component in Oph IRS 63 (and marginally also the broader component) is split into even finer structures as shown in panel (c) of Figure 7. A future study will determine whether the complex shell-like structure is a kinematic feature produced by different episodic events (e.g., Zhang et al. 2019; Fernández-López et al. 2020) or whether a single outflow driven by magnetized winds and their interplay with the environment can produce the observed shells (Shang et al. 2020, 2023).

The redshifted outflow is less prominent and more asymmetric and can be more easily identified in panel (c) of Figure 7, where selected channels are used to emphasize the large-scale ¹²CO emission. Unlike the blueshifted side, the redshifted outflow only shows a broad (large opening angle) component that appears to extend broadly toward the east direction. Although purely speculative, interaction with an inhomogeneous ambient medium may have modified the appearance of the redshifted outflow direction. Alternatively, an outflow precession mechanism may be at play (e.g., Vazzano et al. 2021).

In panel (d) of Figure 7, we overlaid the redshifted component of the outflow with the northern C¹⁸O arc-like emission, as these two appear spatially and morphologically close to each other. The northern C¹⁸O emission appears to be at the outer edge of the ¹²CO outflow and follows its shape closely. Although we only access a two-dimensional projection of the complex system, the C¹⁸O gas likely lies outside the cavity opened by the outflow. If this is the case, then when the outflow opened the cavity observed in ¹²CO, it must have pushed the ambient gas around, creating a localized higher-density and higher-temperature region that we observe as this C¹⁸O arc-like emission. Thieme et al. (2022) recently reported a similar streamer and outflow configuration for the Lupus 3-MMS source.

Lastly, the outflow morphology is related to the orientation of the disk. We can infer that the northeastern edge is the near edge of the disk based on the redshift of the northeast outflow component relative to the source's systemic velocity. Additionally, using the outflow P.A. of ∼240°, we confirm that the blueshifted outflow is perpendicular to the disk direction. However, the redshifted outflow axis appears to be tilted with respect to the disk rotational axis.

Protostellar systems are typically identified as young stars surrounded by a disk and an envelope. The presence of an envelope around a young star distinguishes it from a more evolved Class II or T Tauri star. To determine whether the compact (<2'') disk-like structure is a rotating Keplerian disk and whether the large (∼10'') structure is an infalling envelope, we perform a kinematic analysis using the Spectral Line Analysis/Modeling (SLAM) code (Aso & Sai 2023).

SLAM fits power-law functions to rotational velocity profiles. This method takes a position–velocity (PV) diagram as input and fits a power-law function (described by Equation (2)) to the data points. In this equation, V_m (x) represents the velocity of the model as a function of the impact parameter x (i.e., the distance coordinate in the PV diagram), and p_in is the power-law index of the inner part of the emission that changes to p_in + dp at a breaking radius R_b and breaking velocity V_b . V_sys is the systemic velocity of the system. SLAM selects data points from the input PV in two ways: the ridge mode traces the peak of the emission at different channels, while the edge mode traces the outermost emission contour. SLAM performs the fitting process in a Markov Chain Monte Carlo (MCMC) fashion, and it returns the best-fit parameters along with their formal uncertainties:

$$\begin{eqnarray}{V}_{m}(x)=\left\{\begin{array}{ll}\mathrm{sgn}(x){V}_{b}{\left(\displaystyle \frac{| x| }{{R}_{b}}\right)}^{-{p}_{\mathrm{in}}}+{V}_{\mathrm{sys}} & | x| \leqslant {R}_{b}\\ \mathrm{sgn}(x){V}_{b}{\left(\displaystyle \frac{| x| }{{R}_{b}}\right)}^{-({p}_{\mathrm{in}}+{dp})}+{V}_{\mathrm{sys}} & | x| \gt {R}_{b}\end{array}\right..\end{eqnarray} \tag{ 2 }$$

It is worth mentioning that a Keplerian rotating structure has a rotational profile that falls off as the square root of the radius (i.e., V ∝ r^−0.5), while an infalling envelope has a steeper drop in the rotational velocity with a profile following V ∝ r⁻¹ (e.g., Momose et al. 1998; Yen et al. 2013).

4.1.1. CO Emission

First, we obtained PV cuts along the major axis of the disk for the ¹³CO and C¹⁸O isotopologues (5 pixels wide or ∼05) and show them in Figure 8. In both PV cuts, we observe differential rotation, emission asymmetries caused by envelope contribution, and interferometric filtering. The most prominent signatures of infalling motion can be seen in the ¹³CO PV diagram, as it displays a characteristic emission excess in the upper left non-Keplerian quadrant (quadrant II). This emission has been described in several works, using simple infall models (e.g., Momose et al. 1998; van 't Hoff et al. 2018).

Zoom In Zoom Out Reset image size

We fit the full functional form of Equation (2), which consists of two power laws, to the ¹³CO emission. We did this because we observe emission above 5σ extending to distances >200 au along the major axis. By using Equation (2), we were able to derive the change in the power-law index in the outer parts dp and the braking radius R_b . For C¹⁸O, we only fit a single power law (the first portion of Equation (2)), as the extended emission along the major axis is weaker. In both cases, we assumed a distance of d = 132 pc and a disk inclination of i = 47° and excluded the filtered-out velocity channels from the fit. In Figure 8, we overlaid the best-fit model from SLAM to the PV diagrams along the major axis of ¹³CO and C¹⁸O. The two sets of lines (interior and exterior) in quadrants I and III represent two fitting methods of SLAM: the ridge (circles) and the edge (triangles) modes, respectively.

For both CO isotopologues, we find almost Keplerian inner power-law indices of p_in ∼ 0.6 ± 0.04 using the edge and ridge methods (second column in Table 2). We take the median value of the different methods as the most likely inner power-law index. The stellar mass obtained from SLAM relates to this inner power-law index, and by forcing Keplerian rotation, we derive a mass of M_⋆ = 0.5 ± 0.2 M_⊙ (fifth column in Table 2). Our derived dynamical mass is approximately consistent with the one derived by Brinch & Jørgensen (2013) of M_⋆ = 0.4 M_⊙ after rescaling by our chosen source distance (d) and disk inclination angle (${M}_{\star }\propto \tfrac{d}{{\sin }^{2}({i}_{\mathrm{disk}})}$).

Table 2. Best-fit Parameters from SLAM

Line (Method)	pin	pin + dp	Rb (au)	M⋆ (M⊙)	Vsys (km s–1)
13CO (edge)	0.59 ± 0.01	1.6 ± 0.1	264 ± 3	0.73 ± 0.02	2.80 ± 0.01
13CO (ridge)	0.67 ± 0.02	0.9 ± 0.1	265 ± 22	0.22 ± 0.03	2.82 ± 0.01
C18O (edge)	0.58 ± 0.01	⋯	⋯	0.66 ± 0.01	2.82 ± 0.01
C18O (ridge)	0.60 ± 0.01	⋯	⋯	0.33 ± 0.01	2.81 ± 0.1
SO (edge)	0.48 ± 0.19	1.13 ± 0.23	97 ± 12	0.76 ± 0.19	3.05 ± 0.02
SO (ridge)	0.83 ± 0.06	1.30 ± 0.26	120 ± 39	0.14 ± 0.14	3.08 ± 0.01

Note. The indicated uncertainties are the 14th and 68th statistical uncertainties from the MCMC analysis from SLAM. The adopted disk and stellar properties are obtained by averaging the results from the edge and ridge methods.

Download table as:  ASCIITypeset image

We derive the radius of the rotating gas disk from the fit to the ¹³CO emission, and the radius is defined as the point where the inner and outer power laws change in functional form (R_b ), i.e., where the quasi-Keplerian disk meets the envelope emission. We obtain that such a radius is at R_b = 265 ± 9 au using ridge and edge methods. The calculation of the difference of the power-law index from the inner to the outer parts dp yields significantly different values between the edge and ridge methods. The p_in + dp values measured for the edge and ridge methods are 1.6 and 0.9, respectively, which gives an average of p_in + dp = 1.2 ± 0.4, where the uncertainty is dominated by the standard deviation. A p_in + dp ∼ 1 means that the ¹³CO emission at >260 au (>2'') follows approximately an infalling velocity profile.

4.1.2. SO Emission

Two molecular tracers presented in Section 3.5 exhibit a spiral or arc-like structure. C¹⁸O displays a large-scale arc structure detected at the northern edge of the map, which connects to the disk on the north and northwest sides. There is also a smaller-scale C¹⁸O arm toward the east of the disk (see Figure 5). The SO emission, on the other hand, exhibits three spiral-like structures. As discussed later in Section 5.1, these spiral features might reflect temperature enhancements possibly induced by the impact of the ambient cloud materials on the disk surface. Thus, we would like to test whether their morphology can be reproduced by a collapsing rotating cloud that accretes to a central gravitational object. For this, we follow the analytical solutions proposed by Mendoza et al. (2009). Similar streamer modeling work for other sources has been presented in, e.g., Yen et al. (2017), Yen et al. (2019), Pineda et al. (2020), Thieme et al. (2022), Valdivia-Mena et al. (2022), and Kido et al. (2023).

In brief, the Mendoza et al. (2009) model is described by a gravitational potential produced by the central star and a rotating gas cloud with a finite radius r₀, where particles can collapse starting with a radial velocity v₀. We use spherical coordinates to trace the motion of the infalling particles, where r is the radial coordinate, θ is the polar angle (where the zenith is defined by the disk's rotational axis), and ϕ is the azimuthal angle (which moves in the plane of the disk) (see Figure 1 in Mendoza et al. 2004). For simplicity, we assume that the rotational axis of the cloud and the disk are aligned (i = 47°, P.A. = 149°, the near side is 59° east of north).

The Mendoza et al. (2009) model relaxes the assumption of the classical Ulrich (1976) model, and it differentiates, as the former assumes that the cloud has a finite size and that accreting particles might start its collapse with a nonzero initial radial velocity. These differences translate into trajectories described by conical sections rather than parabolic motions. The assumptions involved in this model are as follows: (1) The pressure gradient of the flow is negligible, (2) the particles at the boundary of the spherical cloud rotate as a rigid body, (3) self-gravity is negligible, and thus the massive central source dominates the motion. The first assumption is valid if the flow is supersonic and the viscosity and radiation effects are small enough. The second one is necessary to give the material a distribution of angular momentum. The final requirement is not difficult to overcome, considering that the central protostar is more massive than the envelope and disk.

To model C¹⁸O, we defined the inner edge of the cloud to be of size R_cloud ∼ 5300 au (20 times the disk's radius), fixed the mass of the protostar to M_⋆ = 0.5 M_⊙, and set the initial specific angular momentum of the particles to ${h}_{0}\,=\sqrt{{R}_{\mathrm{disk}}\,{{GM}}_{\star }}\sim 1250$ au km s⁻¹, where G is the gravitational constant and R_disk = 265 au is the size of the disk. For this molecule, we named the three most prominent arc-like features as East, N. East, and N. West, as shown in the right panel of Figure 9. We manually adjusted the θ and ϕ angles to match the positions and velocities of the three arc-like features. As shown in Figure 9, the streamer model reproduces well the projected positions, as well as the velocity structure along the trajectories of the C¹⁸O line.

Zoom In Zoom Out Reset image size

For SO, we set a smaller disk radius of R_disk = 110 au (following the results from SLAM on SO), and to better match the tightness of the spirals, we had to define a much smaller initial position of the streaming particles of only 275–1320 au to match different spirals. Then, we varied θ and ϕ to visually match the three spirals' features detected in SO. In this way, we were able to reproduce the spatial characteristics of the SO spirals as shown in Figure 10. The velocity profiles of the three streamers (white dashed lines) follow well the emission seen in the PV diagrams, meaning that their velocity structures of the SO spirals can be reproduced by infalling streamer models (right panel of Figure 10).

Zoom In Zoom Out Reset image size

For both C¹⁸O and SO lines, we marked the radial position (with a cyan horizontal line) where the infalling profile should turn into a Keplerian-dominated profile as shown in our SLAM analysis (Section 4.1). Physically this means that the infalling streamers do not reach the midplane of the disk, but instead they likely landed at the disk's surface at distances of several tens of au above the midplane level.

We acknowledge that since we have not performed a mathematical fitting procedure, the values shown in Table 3 might not be the best or even unique solutions to match the observations. In fact, for each model trajectory, small changes in the θ and ϕ angles also produce similarly valid results. These small ranges of solutions could, for example, represent that these streamers have an intrinsic three-dimensional width.

We have previously mentioned that we were unable to confirm whether the ring-like structure is a real ring or the product of continuum oversubtraction. These claims come mainly from a radiative transfer analysis of a three-layer model of dust and gas, where the dust continuum is concentrated in the disk midplane in the form of a thin layer as depicted in Figure 11 of Yen et al. (2019). In such a case, the radiative transfer equation can be written as

$$\begin{eqnarray}\begin{array}{rcl}{I}_{\nu } & = & [{B}_{\nu }({T}_{g})(1+{e}^{-{\tau }_{d}-{\tau }_{g}})-{B}_{\nu }({T}_{d})(1-{e}^{-{\tau }_{d}})]\\ & & (1-{e}^{-{\tau }_{g}}),\end{array}\end{eqnarray} \tag{ 3 }$$

where the subscripts d and g stand for dust and gas, respectively, τ is the optical depth, and B_ν is the Plank function at temperature T. Additionally, when the dust is optically thick (τ_d ≫ 1), as near the stellar position (where T _d ∼ 170 K), it can be written as

$$\begin{eqnarray}&&{I}_{\nu }=[{B}_{\nu }({T}_{g})-{B}_{\nu }({T}_{d})](1-{e}^{-{\tau }_{g}}).\end{eqnarray} \tag{ 4 }$$

A natural conclusion from Equation (4) is that the temperature of SO must be higher than the temperature of the dust everywhere in the disk, as negative SO emission is not observed in any channel map (I_ν > 0). Then, two scenarios are possible to explain the lower SO emission at the disk center. We could have a real ring structure, which would imply a very small SO optical depth toward the center τ_g (i.e., a small (1- ${e}^{-{\tau }_{g}}$) factor). Alternatively, we could have abundant SO gas in the center, i.e., (1 – ${e}^{-{\tau }_{g}}$) ∼ 1, but the apparent emission drop happens as a result of SO being marginally hotter than the dust T_g ≳ T_d . This second scenario would be consistent with the continuum subtraction case, i.e., there is no real ring. Since we cannot distinguish between these two possibilities, we will refer to this feature as a ring. However, we caution the reader that more observations are needed to confirm that this is not a purely continuum subtraction artifact.

An additional support to the ring being a real structure comes from the off-centering between the continuum emission and SO emission depression. We speculate that if the dust were the cause of a central SO continuum oversubtraction, then the SO emission depression should coincide with the dust maximum emission, which is not the case (see Section 3.5).

An SO ring tracing heated material has been proposed to explain ALMA observations of the Class I protostar L1489 IRS (Yen et al. 2014; Yamato et al. 2023), the Class 0 protostar IRAS 16253–2429 (Aso et al. 2023), and the embedded source L1527 IRS (IRAS 04368+2557; Ohashi et al. 2014; Sakai et al. 2014). While these structures are often inferred from PV diagram analysis, a bright ring-like structure in SO can be directly observed in the disk around Oph IRS 63 owing to its moderate inclination. The observed SO enhancement could result from OH reactions after the liberation of S via sputtering (Hartquist et al. 1980), or SO could result directly via sublimation from dust grains when temperatures are ≳40–60 K (e.g., Hasegawa & Herbst 1993). In either case, the temperature increase needed to enhance SO might have originated at the landing zone of the infalling material, i.e., the SO ring could trace an impact zone of the infalling envelope on the disk, as predicted, for example, in the slow shock models of Aota et al. (2015). This scenario would naturally connect the origin of the SO ring with the observed spiral structures. Further evidence of this connection comes from the fact that the outer edge of the SO ring approximately marks a transition between the infalling spiral structures and the inner 100 au that rotate closer to Keplerian velocity. Such kinematic connection has been proposed, for example, for L1527 IRS through the centrifugal radius, where a large portion of the kinetic energy of the infalling material is transformed to rotational energy, allowing the gas to orbit the central source (e.g., Sakai et al. 2014).

Furthermore, some of the features observed in SO have been predicted through magnetohydrodynamical (MHD) models. Infalling and rotating spiral structures have been reported when there is a misalignment between the mean magnetic field and rotation axis of a collapsing core (e.g., Li et al. 2013; Väisälä et al. 2019; Wang et al. 2022; Tu et al. 2023). These MHD simulations have even found a connection between shock regions and spirals of streaming material, as well as formation of ring-like structures due to the circularization of spiraling inflows of matter (Wang et al. 2022).

Since an infalling envelope model reproduces the arc-like structures of C¹⁸O, as well as the spirals in SO, we would like to calculate how much mass these streamers transfer from the inner part of the envelope to the disk. To estimate this, we use the C¹⁸O (2–1) emission, as it is optically thin, and thus we can calculate the amount of mass traced by it. We first confirm that C¹⁸O (2–1) is optically thin using the results obtained by van 't Hoff et al. (2018). They found that if the ratio of the brightness temperature of ¹³CO (2–1) to C¹⁸O (2–1) is higher than 1.5, then C¹⁸O (2–1) emission is optically thin. These results were tested using a range of gas temperatures (20–100 K) and column densities (10¹³–10¹⁸ cm⁻²) assuming standard isotope ratios of [¹²C]/[¹³C] = 77 and [¹⁶O]/[¹⁸O] = 560 (Wilson & Rood 1994). For the gas around Oph IRS 63, we find that the ratio between the peak brightness temperature of the ¹³CO and C¹⁸O lines is between 2 and 3 in the region where most of the CO gas is detected. This result confirms that C¹⁸O is optically thin, and because the ratio is <4, it also confirms that ¹³CO is optically thick (see van 't Hoff et al. 2018).

Since ¹³CO is optically thick, we can use it to calculate the excitation temperature of the gas according to Equation (5). The excitation temperature (T_ex) coincides with the kinetic temperature when the gas is in local thermodynamic equilibrium (LTE). In this equation, h, ν, and k refer to Plank's constant, the frequency of the emission, and the Boltzmann constant, respectively. T¹³ is the peak of the brightness temperature inferred from the ¹³CO (2–1) line, and $J({T}_{\mathrm{bg}})=\tfrac{h\nu /k}{{e}^{h\nu /{{kT}}_{\mathrm{bg}}}-1}$ is the intensity expressed in units of temperature for the microwave background temperature of 2.7 K:

$$\begin{eqnarray}&&{T}_{\mathrm{ex}}=\displaystyle \frac{h\nu /k}{\mathrm{ln}\left(1+\tfrac{h\nu /k}{{T}^{13}+J({T}_{\mathrm{bg}})}\right)}.\end{eqnarray} \tag{ 5 }$$

From this equation, we measure a peak excitation temperature of 21 K in the envelope. Once the excitation temperature is obtained, we then use the C¹⁸O emission to calculate the gas mass of the streamer according to Equations (6) and (7) (following, e.g., Zapata et al. 2014; López-Vázquez et al. 2020). In these equations N(J) is the column density of the J level obtained from Boltzmann's equation, B_e is the rotational constant of the molecule, Z is the partition function of C¹⁸O, m(H₂) is the molecular hydrogen mass, ΔΩ is the angular size of the emission, and d is the distance to Oph IRS 63:

$$\begin{eqnarray}&&{N}_{\mathrm{total}}=N(J)\displaystyle \frac{Z}{(2J+1)}{e}^{\tfrac{{{hB}}_{e}J(J+1)}{{{kT}}_{\mathrm{ex}}}}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&M({{\rm{H}}}_{2})=m({{\rm{H}}}_{2})\displaystyle \frac{X({{\rm{H}}}_{2})}{X({{\rm{C}}}^{18}O)}{\rm{\Delta }}{\rm{\Omega }}{d}^{2}{N}_{\mathrm{total}}.\end{eqnarray} \tag{ 7 }$$

We consider a region that encompasses the extended C¹⁸O emission but excludes the disk's emission, which is concentrated in the inner 265 au. Using a standard abundance of 3 × 10⁻⁷ (Frerking et al. 1982) for C¹⁸O and a B_e value taken from the JPL database (Pickett et al. 1998), we obtain a total C¹⁸O mass of M_envelope ∼ 0.04 M_⊙. The main uncertainty in obtaining this mass comes from the assumed abundance of C¹⁸O with respect to H₂. Additionally, our mass estimate is likely a lower-limit value owing to the flux missed by the interferometer. This is already suggested by the fact that our value is about half of the total envelope mass derived from Brinch & Jørgensen (2013), which was 0.07 M_⊙.

Finally, we estimate a mass infalling rate from the envelope to the disk of ${\dot{M}}_{\mathrm{infall}}\sim 4\times {10}^{-6}$ M_⊙ yr⁻¹. To do this, we divided the mass contained in the infalling envelope by the infalling timescale (${\dot{M}}_{\mathrm{infall}}\sim {M}_{\mathrm{envelope}}/{t}_{\mathrm{infall}}$), where t_infall is the time it takes a parcel of gas at the outermost edge of the C¹⁸O emission to reach the disk. As mentioned in previous works that calculate such mass infall rates (e.g., Ohashi et al. 1997; Takakuwa et al. 2013; Yen et al. 2019), this quantity is just an order-of-magnitude estimate; it could be higher if we consider that part of the large-scale emission can be filtered out by the interferometer, or it could be lower if the infalling velocity is slower than freefall.

As with most young stellar sources, we expect some material to be funneled from the inner parts of the disk to the central (proto)star. This mass transfer is typically traced via hydrogen recombination lines, such as Hα in the optical and Paβ or Brγ in the near-infrared (e.g., Muzerolle et al. 1998; Hartmann et al. 2016). Due to the embedded nature of Oph IRS 63, optical tracers cannot be used to measure the mass accretion rate. Thus, we use the Brγ line at 2.166 μm obtained from the SpeX and iSHELL near-infrared spectra (see Appendix A). These spectra are pseudo-flux-calibrated, and a region around the 2.166 μm Brγ line of the iSHELL spectrum is presented in Figure 11.

Zoom In Zoom Out Reset image size

We deredden the infrared spectrum using our A_v = 24 mag derived from the procedure described in Appendix A. We integrate the Brγ emission in the spectra to calculate the line flux (F_line) and the line luminosity (L_line) using Equation (8):

$$\begin{eqnarray}&&{L}_{\mathrm{line}}=4\pi {d}^{2}{F}_{\mathrm{line}}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\mathrm{log}\left(\displaystyle \frac{{L}_{\mathrm{acc}}}{{L}_{\odot }}\right)=a\mathrm{log}\left(\displaystyle \frac{{L}_{\mathrm{line}}}{{L}_{\odot }}\right)+b\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{acc}}\sim \left(1-\displaystyle \frac{{R}_{\star }}{{R}_{\mathrm{in}}}\right)\displaystyle \frac{{L}_{\mathrm{acc}}{R}_{\star }}{{{GM}}_{\star }}.\end{eqnarray} \tag{ 10 }$$

We convert the line luminosity to accretion luminosity following the prescription in Alcalá et al. (2017) (Equation (9)), where a = 1.19 ± 0.1 and b = 4.02 ± 0.51 for the Brγ line. Finally, we transform this line measurement into a stellar mass accretion rate via Equation (10), assuming R_in = 5R_⋆ (Hartmann & Stauffer 1989) and R_⋆ = 3.0 R_⊙ (this assumes a gravity of $\mathrm{log}g=3.2$, which is consistent with values found around Class 0/I protostars; Doppmann et al. 2005; Greene et al. 2018; Flores et al., submitted). Putting all these numbers together, we measure a (proto)stellar accretion luminosity of L_acc = 0.21 L_⊙ and, consequently, a mass accretion rate of ${\dot{M}}_{\mathrm{acc}}=5\times {10}^{-8}$ M_⊙ yr⁻¹, for an A_v of 24 mag. To estimate the uncertainties of this stellar accretion measurement, we varied the A_v from 20 to 40 mag of extinction and the stellar radius between 2.0 and 3.7 R_⊙ (i.e., $\mathrm{log}g=3.0-3.5$) and obtain (proto)stellar mass accretion rates from ${\dot{M}}_{\mathrm{acc}}=2\times {10}^{-7}$ M_⊙ yr⁻¹ to 2 × 10⁻⁸ M_⊙ yr⁻¹. Although this method has traditionally been applied to T Tauri stars or Class II sources (Alcalá et al. 2017; Manara et al. 2019), it has recently been demonstrated to also work for Class I and flat-spectrum sources (e.g., Fiorellino et al. 2021, 2023).

An independent check on the mass accretion rate for Oph IRS 63 can be obtained following the standard assumption that the bolometric luminosity is the sum of the stellar and the accretion components:

$$\begin{eqnarray}&&{L}_{\mathrm{bol}}={L}_{\star }+{L}_{\mathrm{acc}}.\end{eqnarray} \tag{ 11 }$$

For Oph IRS 63, the bolometric luminosity L_bol is 1.3 L_⊙ (Ohashi et al. 2023), while the stellar luminosity (L_⋆) using R_⋆ = 3.0 R_⊙ and T_eff = 3500 K can be calculated as

$$\begin{eqnarray}&&{L}_{\star }=4\pi \,\sigma \,{R}_{\star }^{2}\,{T}_{\mathrm{eff}}^{4}=1.21\,{L}_{\odot }.\end{eqnarray} \tag{ 12 }$$

This means that the accretion luminosity is of the order of L_acc ∼ 0.1 L_⊙, which corresponds to a mass accretion rate of ${\dot{M}}_{\mathrm{acc}}\sim 3\times {10}^{-8}$ M_⊙ yr⁻¹ according to Equation (10). This value agrees with the lower range of the mass accretion rate derived from the Brγ line.

We observe a significant discrepancy between the mass infall rate ${\dot{M}}_{\mathrm{infall}}\sim 4\times {10}^{-6}$ M_⊙ yr⁻¹ and the stellar accretion rate ${\dot{M}}_{\mathrm{acc}}=(2\mbox{--}20)\times {10}^{-8}$ M_⊙ yr⁻¹, i.e., between the mass infalling from the envelope to the disk and the mass transfer from the disk to the (proto)star.

If these measurements represent their respective "accretion process" and are not just transients, then the natural consequence is a mass buildup in the disk. This is because more mass is infalling into the disk than leaving it, assuming that the mass outflow rate is ∼10% of the mass accretion rate (Ellerbroek et al. 2013).

During the embedded phase, envelope accretion and viscous processes naturally increase a disk's mass and size as the angular momentum is redistributed. However, when a disk becomes too massive, it can become gravitationally unstable, forming asymmetries and spiral features (e.g., Kuffmeier et al. 2018) and sometimes fragmenting (e.g., Vorobyov & Basu 2006; Das & Basu 2022). The Toomre Q parameter (Toomre 1964) is a standard metric used to determine whether a disk is gravitationally stable (Q ≳ 1) or unstable (Q ≲ 1). We calculate the Toomre Q parameter for the disk around Oph IRS 63 using

$$\begin{eqnarray}&&Q\equiv \displaystyle \frac{{c}_{s}{\rm{\Omega }}}{\pi G{\rm{\Sigma }}}.\end{eqnarray} \tag{ 13 }$$

We obtain the angular speed of the gas disk from ${\rm{\Omega }}=\sqrt{{{GM}}_{\star }/{R}_{\mathrm{disk}}^{3}}$ and approximate the disk surface density by ${\rm{\Sigma }}\sim {M}_{\mathrm{disk}}/\pi {R}_{\mathrm{disk}}^{2}$. The gas sound speed (c_s ) is calculated using ${c}_{s}=\sqrt{{k}_{{\rm{B}}}T/\mu {m}_{{\rm{H}}}}$, where k_B is Boltzmann's constant, m_H is the mass of a hydrogen atom, and μ is the mean molecular weight, which is often taken to be μ = 2.3 in disks (e.g., Armitage 2019). As pointed out in Section 3.1, we can use two disk temperatures: T = 20 K or T = 46 K, which yield different dust disk masses. We calculate the total gas disk mass assuming a standard gas-to-dust ratio of g/d = 100 and obtain M_disk ∼ 0.05 M_⊙ and M_disk ∼ 0.02 M_⊙, respectively. We obtain a Toomre Q parameter ranging from Q ∼ 2 to 8 from these two temperatures. It is worth mentioning that the main uncertainties in deriving these values are that the gas-to-dust ratio could be different from 100 and the dust mass is likely a lower limit.

Although it is not possible to predict whether these infall and mass accretion rates will remain at the current levels, it can be calculated that it would take approximately ∼10⁴ yr for the disk to achieve a mass of M_disk ∼ 0.1 M_⊙. If this happens, the Toomre Q parameter would be close to unity, making the disk gravitationally unstable, and possibly triggering episodic accretion bursts.

Comparing the extension of gas and dust in protoplanetary disks is central to the question of radial drift (Facchini et al. 2019; Trapman et al. 2020) and dust settling (Villenave et al. 2020), and ultimately to the planet formation process itself (Armitage 2019). Measuring gas and dust radial sizes can inform us whether grain growth has begun and how dynamically coupled these two components are.

In Section 3.1, we measured the dust disk size as the radius enclosing 95% of the flux and found R_95% = 047 or R_95% = 62 au (at d = 132 pc). Then, in Section 4.1, we use the kinematics of ¹³CO and C¹⁸O to determine that the disk around Oph IRS 63 follows an almost Keplerian rotation up to 265 ± 9 au. By dividing the gas and dust disk sizes, we find a gas-to-dust size ratio of ${R}_{\mathrm{gas}}/{R}_{95 \% }^{\mathrm{dust}}\unicode{x00303}4.3$. This gas-to-dust size ratio might not be directly comparable to studies in Class II disks because, in such cases, both the dust and gas disks' radii are measured based on a fixed amount of flux enclosed in a certain radius (e.g., 90% or 95%). However, due to the embedded nature of our source, we could not perform such a measurement for the gas and instead used a different definition than for the dust. In any case, for the sake of comparison, the gas-to-dust size ratio of the disk around Oph IRS 63 is slightly larger than the average value obtained in the more evolved T Tauri stars of 2.9 ± 1.2 (Long et al. 2022) and compared to other values found in the literature (e.g., Ansdell et al. 2018; Trapman et al. 2020).

Using thermochemical models, Trapman et al. (2019) and Sanchis et al. (2021) found that a gas-to-dust size ratio of R_gas/R_dust > 4 is a clear sign of dust evolution and radial drift in disks. They also showed that in cases where R_gas/R_dust < 4 optical depth effects alone could reproduce these observational signatures. In the latter case, identifying dust evolution from R_gas/R_dust requires detailed modeling of the disk, including the total CO content. For the case of Oph IRS 63, the gas-to-dust size ratio is slightly above the critical value of 4. Thus, it is possible that radial drift effects alone could account for the size difference. However, as mentioned above, the results obtained for Oph IRS 63 might not be directly comparable to studies in Class II sources. Future studies with detailed modeling of the dust and gas around Class I sources will be necessary to clarify whether the radial size difference observed for Oph IRS 63 implies dust growth and radial drift.

Figure 12 shows a schematic view of the large-scale emission of the system. The envelope is represented by the gray rotating flattened structure observed in ¹³CO at radial offset ≥2''. The red arrows indicate the portion of the envelope receding, while the blue ones mark emission toward us. The envelope should be at least as large as the field of view (FOV) of our observations (∼25'') based on the streamer structure seen toward the northeast in C¹⁸O, the extended emission seen at low intensity in ¹³CO, and the extended emission observed close to the systemic velocity of ¹²CO (more specifically in channel maps with velocities of 0.96, 1.6, 3.5, and 4.14 km s⁻¹; see Figure B1 of Appendix B).

Zoom In Zoom Out Reset image size

The almost Keplerian rotating disk structure seen in the observations with angular scale ≤2'' is represented in colors with an oval shape at the center of the image. Assuming that the orientation of the envelope is similar to the disk (i_disk = 47° and P.A. = 149^◦), and given the outflow directions, we expect the southwest part of the envelope to lie behind the blueshifted outflow, while the northeast part should be in front (from our perspective) of the redshifted outflow. This geometry could explain the lack of envelope emission toward the west and southwest seen in the ¹³CO peak intensity and moment 0 map of Figure 4. Additionally, this geometry partially explains why the redshifted outflow is weaker and has a simpler structure than the blueshifted one owing to self-absorption and emission being more resolved out. The figure also depicts the blue- and redshifted wide-angle outflows, along with the narrow blueshifted outflow emission.