© ESO, 2014

1. Introduction

Rotationally supported accretion disks (RSDs) are thought to form very early during the star formation process (see, e.g., Bodenheimer 1995). The RSD transports a significant fraction of the mass from the envelope onto the young star and eventually evolves into a protoplanetary disk as the envelope dissipates. The presence of an RSD affects not only the physical structure of the system but also the chemical content and evolution as it promotes the production of more complex chemical species by providing a longer lifespan at mildly elevated temperature (Aikawa et al. 2008, 2011; Visser et al. 2011). Although the standard picture of star formation predicts early formation and evolution of RSDs, theoretical studies suggest that the presence of magnetic fields prevents the formation of RSDs in the early stages of star formation (e.g., Galli & Shu 1993; Chiang et al. 2008; Li et al. 2011; Joos et al. 2012). Unstable flattened disk-like structures are formed in such simulations, and indeed, disks in the embedded phase have been inferred through continuum observations (e.g., Looney et al. 2003; Jørgensen et al. 2009; Enoch et al. 2011; Chiang et al. 2012). However, with only continuum data, it is difficult to distinguish between such a feature and an RSD. Thus, spatially and spectrally resolved observations of the gas are needed to unravel the nature of these embedded disks.

Only a handful of studies have explored the change in the velocity profiles from the envelope to the disk with spectrally resolved molecular lines (Lee 2010; Yen et al. 2013; Murillo et al. 2013). It is essential to differentiate between the disk and the infalling envelope through such analysis. This paper presents IRAM Plateau de Bure Interferometer (PdBI) observations of ¹³CO and C¹⁸O J = 2–1 toward four Class I sources in Taurus (d = 140 pc) at higher angular resolution (~0.8″= 112 AU diameter at 140 pc) and higher sensitivity than previously possible. We aim to study the velocity profile as revealed through these molecular lines and constrain the size of embedded RSDs toward these sources.

The embedded phase can be divided into the following stages (Robitaille et al. 2006): Stage 0 (M_env > M_⋆), Stage I (M_env < M_⋆ but M_disk < M_env) and Stage II (M_env < M_disk), where M_⋆, M_env and M_disk are masses of the central protostar, envelop and RSD respectively. This is not to be confused with Classes, which are observationally derived evolutionary indicators (Lada & Wilking 1984; Andre et al. 1993) that do not necessarily trace evolutionary stage due to geometrical effects (Whitney et al. 2003; Crapsi et al. 2008; Dunham et al. 2010). However, since all Class I sources discussed in this paper also turn out to be true Stage I sources, we use the Class notation throughout this paper for convenience.

RSDs are ubiquitous once most of the envelope has been dissipated away (Class II) with outer radius R_out up to approximately 200 AU (Williams & Cieza 2011, and references therein). On the other hand, the general kinematical structure of deeply embedded protostars on small scales is still not well understood. There is growing evidence of rotating flattened structures around Class I sources (Keene & Masson 1990; Hayashi et al. 1993; Ohashi et al. 1997a,b; Brinch et al. 2007; Lommen et al. 2008; Jørgensen et al. 2009; Lee 2010, 2011; Takakuwa et al. 2012; Yen et al. 2013), but the question remains how and when RSDs form in the early stages of star formation and what their sizes are.

Class I young stellar objects (YSOs) are ideal targets to search for embedded RSDs. At this stage, the envelope mass has substantially decreased such that the embedded RSD dominates the spatially resolved CO emission (Harsono et al. 2013). CO is the second most abundant molecule and is chemically stable, thus the disk emission should be readily detected. Furthermore, Harsono et al. (2013) showed that the size of the RSD can be directly measured by analyzing the velocity profile observed through spatially and spectrally resolved CO observations. The sources targeted here are TMC1 (IRAS 04381+2540), TMR1 (IRAS 04361+2547), L1536 (IRAS 04295+2251), and TMC1A (IRAS 04365+2535). These Class I objects have been previously observed by the Submillimeter Array (SMA) and the Combined Array for Research in Millimeter-wave Astronomy (CARMA) at lower sensitivity and/or resolution (e.g., Jørgensen et al. 2009; Eisner 2012; Yen et al. 2013). Embedded rotating flattened structures have been proposed previously around TMC1 and TMC1A (Ohashi et al. 1997a; Hogerheijde et al. 1998; Brown & Chandler 1999; Yen et al. 2013). Thus, we target these sources to determine the presence of Keplerian disks, constrain their sizes near the end of the main accretion phase on scales down to 100 AU and compare this with the dust structure.

This paper is laid out as follows. Section 2 presents the observations and data reduction. The continuum and line maps are presented in Sect. 3. In Sect. 4, disk masses and sizes are obtained through continuum visibility modelling. Position–velocity diagrams are then analyzed to determine the size of rotationally supported disks. Section 5 discusses the implications of the observed rotational supported disks toward Class 0 and I sources. Section 6 summarizes the conclusions of this paper.

2. Observations

2.1. IRAM PdBI observations

The sources were observed with IRAM PdBI in the two different configurations tabulated in Table A.1. Observations with an on-source time of ~3.5 h per source with baselines from 15 to 445 m (11 kλ to 327 kλ) were obtained toward TMC1A and TMR1 in March 2012. Additional track-sharing observations of L1536 and TMC1 were obtained in March-April 2013 with an on-source time of ~3 h per source and baselines between 20 to 450 m (15 kλ to 331 kλ). The receivers were tuned to 219.98 GHz (1.36 mm) in order to simultaneously observe the ¹³CO (220.3987 GHz) and C¹⁸O (219.5603 GHz) J = 2–1 lines. The narrow correlators (bandwidth of 40 MHz ~54 km s^-1) were centered on each line with a spectral resolution of 0.078 MHz (0.11 km s^-1). In addition, the WideX correlator was used, which covers a 3.6 GHz window at a resolution of 1.95 MHz (2.5−3 km s^-1) with 6 mJy beam^-1 channel^-1 RMS.

The calibration and imaging were performed using the CLIC and MAPPING packages of the IRAM GILDAS software¹. The standard calibration was followed using the calibrators tabulated in Table A.1. The data quality assessment tool flagged out any integrations with significantly deviating amplitude and/or phase and the continuum was subtracted from the line data before imaging. The continuum visibilities were constructed using the WideX correlator centered on 219.98 GHz (1.36 mm). The resulting beam sizes and noise levels for natural and uniform weightings can be found in Table 1. The uniform weighted images will be used for the analysis in the image space.

2.2. WideX detections

The WideX correlators cover frequencies between 218.68–222.27 GHz. Strong molecular lines that are detected toward all of the sources within the WideX correlator include C¹⁸O 2–1 (219.5603541 GHz), SO 5₆–4₅ (219.9488184 GHz), and ¹³CO 2–1 (220.3986841 GHz). The WideX spectra and integrated flux maps are presented in Appendix C but not analyzed here.

3. Results

3.1. Continuum maps

Strong 1.36 mm continuum emission is detected toward all four sources up to ~450 m baselines (~330 kλ, which corresponds to roughly 80 AU in diameter at the distance of Taurus). To estimate the continuum flux and the size of the emitting region, we performed elliptical Gaussian fits to the visibilities (>25 kλ). The best fit parameters are listed in Table 2. Our derived PdBI 1.36 mm fluxes are <30% lower than the single-dish fluxes tabulated in Motte & André (2001). The only exception is L1536 for which 72% of the singe-dish flux is recovered, which indicates a lack of resolved out large-scale envelope around L1536. The deconvolved sizes indicate a large fraction of the emission is within <100 AU diameter, consistent with compact flattened dust structures.

Figure 1 presents the uniform weighted continuum maps with red circles indicating the continuum positions (Table 2). The total flux of TMC1A is ~30% lower than in Yen et al. (2013) within a 4.0″ × 3.5″ beam, which indicates that some extended structure is resolved-out in our PdBI observations. The peak fluxes of our cleaned maps agree to within 15% with those reported by Eisner (2012) except toward L1536, which is a factor of two lower in our maps. However, our image toward L1536 (Fig. 2) shows that there are two peaks whose sum is within 15% of the peak flux reported in Eisner (2012). The position in Table 2 is centered on the Western peak and the two peaks are separated by ~70 AU. The analysis in the next section (Sect. 4) is with respect to the Western peak. The maps indicate the presence of elongated flattened structure perpendicular to the outflow direction found in the literature (TMC1: 0°; TMR1: 165°; TMC1A: 155°; L1536: None², Ohashi et al. 1997a; Hogerheijde et al. 1998; Brown & Chandler 1999). Such elongated dust emission is indicative of flattened disk-like structures.

3.2. Molecular lines: maps and morphologies

The ¹³CO and C¹⁸O 2–1 lines are detected toward all sources. Figure 3 presents the spectra integrated within a 1″ box around the continuum position. The integrated line flux densities within a 5″ box are tabulated in Table 3. The molecular line emission is strongest toward TMC1A and weakest toward L1536, while most of the emission around the systemic velocity is resolved out. This can be seen in the ¹³CO channel map toward TMC1A (Figs. 3 and 4, see Appendix A for channel maps toward other sources) where most of the emission is absent near υ_lsr = 6.6 km s^-1. In contrast to ¹³CO, the C¹⁸O line is detected around the systemic velocities (±2 km s^-1) as shown in Fig. 3 except toward TMC1A, which shows emission within the velocity range comparable to that of ¹³CO.

Single-dish observations of ¹³CO J = 3–2 (except for L1536, which is the J = 2–1 transition) and C¹⁸O J = 2–1 are also shown in Fig. 3 (right). The comparison indicates that a large part of the emission at the systemic velocity originating from the large-scale envelope is filtered out in our PdBI observations. The PdBI spectra within a 1″ box around each source are also broader than the single-dish spectra. Thus, our PdBI observations are probing the kinematics within the inner few hundreds of AU.

The zeroth moment (i.e., velocity-integrated flux density) maps, are shown in Fig. 5 toward all four sources. The morphologies of the ¹³CO and C¹⁸O maps are different toward 3 out of the 4 sources due to C¹⁸O emission being weaker. Only toward TMC1A, which has the strongest C¹⁸O emission, do the C¹⁸O and ¹³CO show a similar morphologies. In general, the zeroth moment maps indicate the presence of flattened gas structures that are perpendicular to the outflow direction except for TMR1. In most cases, the C¹⁸O emission is weak and compact toward the continuum position as revealed by our higher S/N spectrally unresolved spectra taken with the WideX correlator (bottom panels of Fig. 5).

Moment one (i.e., intensity-weighted average velocity) maps are presented in Fig. 6 and show the presence of velocity gradients in the inner few hundreds of AU. These velocity gradients are perpendicular to the large-scale outflow. Large-scale rotation around TMC1A was reported previously by Ohashi et al. (1997a) extending up to 2000 AU. In addition, Hogerheijde et al. (1998) also detected rotation around TMC1. The high sensitivity and spatial resolution of our observations allow us to map the kinematics on much smaller scales and therefore unravel the nature of these gradients.

4. Analysis

4.1. Continuum: dust disk versus envelope

The advantage of an interferometric observation is that one can disentangle the compact disk-like emission from the large-scale envelope. Thus, an estimate of the disk mass can be obtained using both the single-dish and the spatially resolved continuum observations (visibility amplitudes). Figure 7 shows the circularly averaged visibility amplitudes as a function of the projected baselines in kλ (assuming spherically symmetric emission). All sources show strong emission on long baselines that is near constant with increasing baseline. The continuum visibilities of TMR1 and TMC1A are consistent with the scaled (F_ν ∝ ν^α) 1.1 mm SMA observations (red circles, Jørgensen et al. 2009) with α = 2.5. An excess of continuum emission at the longest baselines compared with the expected noise is seen toward all sources, which corresponds to unresolved continuum structure at scales smaller than 320 kλ (50 AU radius). Moreover, rapidly increasing amplitudes at short baselines are typically associated with the presence of a large-scale envelope. This combination of features in the visibilities is a typical signature of embedded disks (e.g., Looney et al. 2003; Jørgensen et al. 2005; Enoch et al. 2009).

The first important step is quantifying the large-scale envelope contribution to the uv amplitudes. Using radiative transfer models, Jørgensen et al. (2009) estimated that a spherical envelope contributes at most 4% at 50 kλ at 1.1 mm, which translates to 2% at 1.36 mm. By using the single-dish 850 μm (S_15″) and the 50 kλ fluxes (S_50kλ), the envelope, S_env, and disk, S_disk, flux densities can be estimated by solving the following system of equations (Eqs. (2) and (3) in Jørgensen et al. 2009):

$\begin{array}{ccc}& & {\mathit{S}}_{\mathrm{50}\mathrm{k}\mathit{\lambda }}\mathrm{=}{\mathit{S}}_{\mathrm{disk}}\mathrm{+}\mathit{c}\mathrm{\times }{\mathit{S}}_{\mathrm{env}}\\ & & {\mathit{S}}_{\mathrm{15}″}\mathrm{=}{\mathit{S}}_{\mathrm{disk}}\mathrm{(}\mathrm{1.36}\mathit{/}\mathrm{0.85}{\mathrm{)}}^{\mathit{\alpha }}\mathrm{+}{\mathit{S}}_{\mathrm{env}}\mathit{.}\end{array}$

We have adopted α ≃ 2.5 to translate the 850 μm disk flux to 1.36 mm (Jørgensen et al. 2007) and c = 0.02 as the fractional envelope contribution. The derived disk masses within 100 AU diameter are tabulated in Table 2 using the OH5 (Ossenkopf & Henning 1994) opacities at 1.36 mm and a dust temperature of 30 K. Our values for TMR1 and TMC1A are consistent with Jørgensen et al. (2009). The envelope masses listed in Table 2 were estimated using Eq. (1) of Jørgensen et al. (2009). Three of the sources (TMC1, TMR1, TMC1A) have similar envelope masses of ~0.1  M_⊙, which are higher than the disk masses. In contrast, at most a tenuous envelope seems to be present around L1536. Using the disk to envelope mass ratio as an evolutionary tracer, L1536 is the most evolved source in our sample, which is consistent with its high bolometric temperature. Thus, these sources are indeed near the end of the main accretion phase where a significant compact component is present in comparison to their large-scale envelope.

In interferometric observations toward YSOs, the lack of short baselines coverage may result in underestimation of the large-scale emission and, consequently, overestimating the compact structure. Kristensen et al. (2012) calculated 1D spherically symmetric models derived using the 1D continuum radiative transfer code DUSTY (Ivezic & Elitzur 1997) for a number of sources, including TMC1A, TMR1, and TMC1. The model is described by a power law density structure with an index p (n = n₀(r/r₀)^−p) where n is the number density and with an inner boundary r₀ at T_dust = 250 K. The parameters are determined through fitting the observed spectral energy distribution (SED) and the 450 and 850 μm dust emission profiles. These parameters are tabulated in Table B.1. The models in Kristensen et al. (2012) do not take into account the possible contribution by the disk to the submillimeter fluxes and the envelope masses are consequently higher than ours by up to a factor of 2. Using these models, the missing flux can be estimated by modelling the large-scale continuum emission as it would be observed by the interferometer (blue line in Fig. 7). We can therefore place better constraints on the compact flattened structure.

4.2. Constraining the dust disk parameters

The goal of this section is to estimate the size and mass of the compact flattened dust structure or the disk using a number of more sophisticated methods than Eqs. (1) and (2) to compare with the estimates from the gas in the next section. We define the disk component to be any compact flattened structure that deviates from the expected continuum emission due to the spherical envelope model. In the following sections, the disk sizes will be estimated by fitting the continuum visibilities with the inclusion of emission from the large-scale structure described above. As the interferometric observation recovers a large fraction of the single-dish flux toward L1536, its large-scale envelope contribution is assumed to be insignificant, and therefore not included. The best-fit parameters are found through χ² minimization on the binned amplitudes. The range of values of R_disk, M_disk, i and PA shown in Table 4 correspond to models with χ² between ${\mathit{\chi }}_{\min }^{\mathrm{2}}$ and ${\mathit{\chi }}_{\min }^{\mathrm{2}}\mathrm{+}\mathrm{15}$.

4.2.1. Uniform disk

The first estimate on the disk size and mass is obtained by fitting a uniform disk model assuming optically thin dust emission (see Fig. 8, Eisner 2012). Uniform disks are described by a constant intensity I within a diameter θ_ud whose flux is given by $\mathit{F}\mathrm{=}{}^{\mathrm{\int }}\mathit{I}\cos \mathit{\theta }\mathrm{d\Omega }\mathrm{=}\mathit{I\pi }{{}^{\mathrm{\right(}}\frac{{\mathit{\theta }}_{\mathrm{ud}}}{\mathrm{2}}^{\mathrm{\left)}}}^{\mathrm{2}}$. The visibility amplitudes are given by: $\mathit{V}\mathrm{\left(}\mathit{u,v}\mathrm{\right)}\mathrm{=}\mathit{F}\mathrm{\times }\mathrm{2}\frac{{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}\mathit{\pi \theta }{\mathit{r}}_{\mathit{uv}}\mathrm{\right)}}{\mathit{\pi \theta }{\mathit{r}}_{\mathrm{uv}}}\mathit{,}$(3)where J₁ is the Bessel function of order 1 and r_uv is the projected baseline in terms of λ. Deprojection of the baselines follow the formula given in Berger & Segransan (2007): $\begin{array}{ccc}{\mathit{u}}_{\mathrm{PA}}& \mathrm{=}& \mathit{u}\cos \mathrm{PA}\mathrm{+}\mathit{v}\sin \mathrm{PA}\mathit{,}\\ {\mathit{v}}_{\mathrm{PA}}& \mathrm{=}& \mathrm{-}\mathit{u}\sin \mathrm{PA}\mathrm{+}\mathit{v}\cos \mathrm{PA}\mathit{,}\\ {\mathit{r}}_{\mathit{uv,i}}& \mathrm{=}& \sqrt{\mathrm{(}{\mathit{u}}_{\mathrm{PA}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathrm{PA}}^{\mathrm{2}}\cos \mathrm{(}\mathit{i}{\mathrm{)}}^{\mathrm{2}}\mathrm{)}}\mathit{,}\end{array}$where PA is the position angle (East of North) and i is the inclination.

The best fit parameters are tabulated in Table 4 with disk masses estimated using a dust temperature of 30 K and OH5 opacities (Table 5 of Ossenkopf & Henning 1994) within $\frac{{\mathit{\theta }}_{\mathrm{ud}}}{\mathrm{2}}$ radius. The disk radii vary between 7 to 171 AU. The smallest disk size is found toward TMR1, which suggests that most of its continuum emission is due to the large-scale envelope. Moreover, disk sizes around TMR1 and TMC1A are 65%–90% lower than those reported by Eisner (2012, Table 3) since we have included the spherical envelope contribution. This illustrates the importance of including the large-scale contribution in analyzing the compact structure.

4.2.2. Gaussian disk

The next step is to use a Gaussian intensity distribution, which is a slightly more realistic intensity model that represents the mm emission due to an embedded disk. The visibility amplitudes are described by the following equation: $\mathit{V}\mathrm{\left(}\mathit{u,v}\mathrm{\right)}\mathrm{=}\mathit{F}\mathrm{\times }\exp \left(\mathrm{-}\frac{\mathit{\pi }{\mathit{\theta }}_{\mathrm{G}}{\mathit{r}}_{\mathit{uv}}}{\mathrm{4}\log \mathrm{2}}\right)\mathit{,}$(7)where θ_G is the FWHM of the Gaussian distribution and F is the total flux.

The difference between the Gaussian fit presented in Table 2 and this section is the inclusion of the large-scale emission. In Table 2, the Gaussian fit gives an estimate of the size of emission in the observed total image, while this section accounts for the simulated large-scale envelope emission in order to constrain the compact flattened structure. The free-parameters are similar to the uniform disk except that the size of the emitting region is defined by full-width half maximum (FWHM). For a Gaussian distribution, most of the radiation (95%) is emitted from within 2σ_G ~ 0.85 × FWHM, thus we define the disk radius to be 0.42 × FWHM. The best-fit parameters are tabulated in Table 4 and the corresponding emission is included in Fig. 7. The resulting disk sizes are slightly lower than those derived using a uniform disk.

The visibility amplitudes in Fig. 7 show the presence of an unresolved point source component at long baselines indicated by the arrow in the first panel. Jørgensen et al. (2005) noted that a three components model (large-scale envelope + Gaussian disk + point-source) fits the continuum visibilities. By using such models, the disk sizes generally increase by 20–30 AU with the addition of an unresolved point flux, which is comparable to the uniform disk models. Such models were first introduced by Mundy et al. (1996) since the mm emission seems to be more centrally peaked than a single Gaussian. The unresolved point flux is more likely due to the unresolved disk structure since the free-free emission contribution is expected to be low toward Class I sources (Hogerheijde et al. 1998).

4.2.3. Power-law disk

The next step in sophistication is to fit the continuum visibilities with a power-law disk structure. The difference is that the first two disk models are based on the expected intensity profile of the disk, while the power-law disk models use a more realistic disk structure. Similar structures were previously used by Lay et al. (1997, see also Lay et al. 1994; Mundy et al. 1996; Dutrey et al. 2007; Malbet et al. 2005. We adopt a disk density structure described by Σ_disk = Σ_50  AU(R/50  AU)^−p where Σ_50  AU is the surface density at 50 AU with a temperature structure given as T_disk = 1500(R/0.1  AU)^−q. The inner radius is fixed at 0.1 AU with a dust sublimation temperature of 1500 K, while the disk outer radius is defined at 15 K. The free parameters are p, q, Σ_50  AU, i and PA where i is the inclination with respect of the plane of the sky in degrees (0° is face-on) and PA is the position angle (East of North). The visibilities are constructed by considering the flux from a thin ring given by: $\mathrm{d}\mathit{S}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }\cos \mathit{i}}{{\mathit{d}}^{\mathrm{2}}}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}\mathit{T}\mathrm{\right)}\left(\mathrm{1}\mathrm{-}{\exp }^{\mathrm{-}\mathit{\tau }}\right)\mathit{R}\mathrm{d}\mathit{R,}$(8)where d is the distance, i is the inclination (0° is face-on), $\mathit{\tau }\mathrm{=}\frac{\mathrm{\Sigma }\left(\mathit{R}\right)\mathit{\kappa }}{\cos \mathit{i}}$ is the optical depth and B_ν(T) is the Planck function. For the dust emissivity (κ), we adopt the OH5 values. Since the model is axisymmetric, the visibilities are given by a one-dimensional Hankel transform which is given by the analytical expression: $\mathit{V}\mathrm{(}\mathit{u,}{\mathit{v}}^{\mathrm{)}}\mathrm{=}\mathrm{d}\mathit{S}\mathrm{\times }{\mathit{J}}_{\mathrm{0}}\mathrm{\left(}\mathrm{2}\mathit{\pi R}{\mathit{r}}_{\mathit{uv}}\mathrm{\right)}\mathit{,}$(9)where J₀ is the Bessel function of order 0. The total amplitude at a given r_uv is integrated over the whole disk.

The first advantage of such modelling is the treatment of the optical depth which gives better mass approximations. Secondly, the unresolved point-flux is already included in the model by simulating the emission from the unresolved disk component. Figure 8 shows the normalized visibilities for a fixed p = 0.9 and a number of R_disk. The best-fit parameters are tabulated in Table 4.

The derived disk masses obtained by integrating out to ~50 AU radius are similar to those from the 50 kλ point (Table 2). This is expected as the 1.3 mm continuum emission (κ_OH5 = 0.83 cm² g^-1) is on average nearly optically thin (Jørgensen et al. 2007). The optical depth can be up to 0.4 at radii <10 AU for the best-fit power-law disk models with i > 70°. A wide range of “best-fit” masses is found toward TMR1, which is mainly due to the uncertainties in the inclination.

More recently, Eisner (2012) presented embedded disk models derived from continuum data by fitting the I band image, SED and 1.3 mm visibilities. Our continuum disk radii toward TMC1A and TMR1 are consistent with their results. However, the high quality of our data allow us to rule out disks with R_out > 100 AU toward TMR1. On the other hand, our disk radius toward L1536 is a factor of 3 higher than that reported in Eisner (2012) due to the difference in the treatment of the large-scale envelope.

In summary, different methods have been used to constrain the disk parameters from thermal dust emission. The first two models, Gaussian and uniform disk models, are based on the intensity profile, while the power-law disk models use a simplified disk structure with given density and temperature profiles. The intensity based models find the smallest disk sizes, while the power-law disk models predict up to a factor of two larger disk sizes. On the other hand, the disk masses from these fitting methods are similar. The advantage of the power-law disk modeling is the determination of the density and temperature structures including the unresolved inner part of the disk. The continuum images provide no kinematic information, however, so the important question is whether these flattened structures are rotationally supported disks and if so, whether the continuum sizes agree with those of RSDs.

4.3. Line analysis: Keplerian or not?

The nature of the velocity gradient can already be inferred from the moment maps. First, the ¹³CO integrated flux maps shown in Fig. 5 indicate the presence of flattened gas structures perpendicular to the outflow direction. Second, the velocity gradient is also oriented perpendicular to the outflow within hundreds of AU, which is a similar size scale to the compact dust structure measured in Sect. 4.2. On the other hand, the velocity gradient does not show a straight transition from the red to the blue shifted component as expected from a Keplerian disk. Taking TMC1A as an example, the blue-shifted emission starts at the North side and skews inward toward the continuum position. Such a skewed moment 1 map is expected for a Keplerian disk embedded in a rotating infalling envelope (Sargent & Beckwith 1987; Brinch et al. 2008). Thus, the moment maps indicate the presence of a rotating infalling envelope leading up to a rotationally dominated structure.

In image space, position velocity (PV) diagrams are generally used to determine the nature of the velocity gradient. Analysis of such diagrams starts by dividing it into four quadrants around the υ_lsr and source position. A rotating structure occupies two of the four quadrants that are symmetric around the center, while an infalling structure will show emission in all quadrants (e.g., Ohashi et al. 1997a; Brinch et al. 2007). In addition, one can infer the presence of outflow contamination by identifying a velocity gradient in PV space along the outflow direction (e.g., Cabrit et al. 1996).

Figure 9 presents PV diagrams at a direction perpendicular to the outflow direction; for L1536, the direction of elongation of the continuum was taken. In general, the PV diagrams are consistent with Keplerian profiles, except for TMR1. The PV diagram toward TMR1 seems to indicate that the ¹³CO line is either dominated by infalling material or the source is oriented toward us, which is difficult to disentangle. However, recent ¹²CO 3–2 and 6–5 maps presented in Yıldız et al., (in prep.) suggest that the disk is more likely oriented face-on (i < 15°). Furthermore, a velocity gradient is present in the ¹³CO PV analysis both along and perpendicular to the outflow direction indicating minor outflow contamination for this source.

The focus of this paper is to differentiate between the infalling gas and the rotationally supported disk. Such analysis in image space is not sensitive to the point where the infalling rotating material (υ ∝ r^-1) enters the disk and becomes Keplerian (υ ∝ r^-0.5). Thus, additional analysis is needed to differentiate between the two cases.

As pointed out by Lin et al. (1994), infalling gas that conserves its angular momentum exhibits a steeper velocity profile (υ ∝ r^-1) than free-falling gas (υ ∝ r^-0.5). In the case of spectrally resolved optically thin lines, the peak position of the emission corresponds to the maximum possible position of the emitting gas (Sargent & Beckwith 1987; Harsono et al. 2013). On the other hand, molecular emission closer to the systemic velocity is optically thick and, consequently, the inferred positions are only lower limits. Thus, with such a method, one can differentiate between the Keplerian disk and the infalling rotating gas. Moreover, Harsono et al. (2013) suggested that the point where the two velocity profiles meet corresponds to the size of the Keplerian disk, R_k. In the following sections, we will attempt to constrain the size of the Keplerian disk using this method, which we will term the uv-space PV diagram.

The peak positions are determined by fitting the velocity resolved visibilities with Gaussian functions (Lommen et al. 2008; Jørgensen et al. 2009). It can be seen from the channel maps (Fig. 4 and Appendix A) that a Gaussian brightness distribution is a good approximation in determining the peak positions of the high velocity channels.

To characterize the profile of the velocity gradient, the peak positions are projected along the velocity gradient. This is done by using the following transformation $\begin{array}{ccc}{\mathit{x}}_{\mathrm{PA}}& \mathrm{=}& \mathit{x}\cos \mathrm{\left(}\mathrm{PA}\mathrm{\right)}\mathrm{+}\mathit{y}\sin \mathrm{\left(}\mathrm{PA}\mathrm{\right)}\mathit{,}\\ {\mathit{y}}_{\mathrm{PA}}& \mathrm{=}& \mathrm{-}\mathit{x}\sin \mathrm{\left(}\mathrm{PA}\mathrm{\right)}\mathrm{+}\mathit{y}\cos \mathrm{\left(}\mathrm{PA}\mathrm{\right)}\mathit{,}\end{array}$where x,y are the peak positions and x_θ,y_θ are the rotated positions. A velocity profile (υ ∝ r^−η) is then fitted to a subset of the red- and blue-shifted peaks to determine the best velocity profile.

4.3.1. TMC1A

TMC1A shows the greatest promise so far for an embedded Keplerian disk around a Class I protostar . The moment 1 maps (Fig. 6) indicate a gradual change from blue- to red-shifted gas in the inner few hundreds of AU, which is consistent with a Keplerian disk. Figure 10 presents the C¹⁸O velocity profiles using the nominal υ_lsr and PA perpendicular to the outflow direction along with small deviations from those values, which correspond to the uncertainty in outflow direction. In the absence of foreground clouds, one expects that the blue and red-shifted peaks overlap for a rotating system, which is the case if one adopts υ_lsr =  + 6.4 km s^-1. The C¹⁸O and C¹⁷O observations presented in Jørgensen et al. (2002) indicate a υ_lsr =  + 6.6 km s^-1, however υ_lsr =  + 6.4 km s^-1 is consistent with the N₂H⁺ observations toward the nearby L1534 core (Caselli et al. 2002). For comparison, the ¹³CO velocity profile with υ_lsr =  + 6.4 km s^-1 and PA perpendicular to the outflow direction is shown in Fig 11. Both of the ¹³CO and C¹⁸O lines show a velocity profile close to υ ∝ r^-0.5, which indicates that the gas lines trace similar rotationally supported structures. The typical uncertainty in the power-law slope is ± 0.2. In the region where the red- and blue-shifted peaks overlap, a clear break from r^~−0.5 to a steeper slope of r^-1 can be seen at ~100 AU. This break may shift inward at other PAs and υ_lsr to 80 AU. However, it is clear that the flattened structure outside of 100 AU flattened structure is dominated by infalling rotating motions. Thus, the embedded RSD toward TMC1A has a radius between 80–100 AU.

To derive the stellar mass associated with the rotationally dominated region, simultaneous fitting to the PV diagrams in both image and uv space has been performed. The best-fit stellar mass associated with the observed rotation toward TMC1A is 0.53$\begin{array}{c}\mathrm{+}\mathrm{0.2}\\ -0.1\end{array}$M_⊙ at i = 55° ± 10°.

Additional analysis is performed for TMC1A with a simple radiative transfer model to confirm the extent of the Keplerian disk following a method similar to that used by Murillo et al. (2013) for the VLA 1623 Class 0 source. Such fitting is only performed for this source to investigate whether they give similar results. A flat disk model given by: $\begin{array}{ccc}\mathit{T}& \mathrm{=}& {\mathit{T}}_{\mathrm{0}}{\mathrm{(}\mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{0}}\mathrm{)}}^{\mathrm{-}\mathit{q}}\mathit{,}\\ \mathit{N}& \mathrm{=}& {\mathit{N}}_{\mathrm{0}}{\mathrm{(}\mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{0}}\mathrm{)}}^{\mathrm{-}\mathit{p}}\mathit{,}\end{array}$where N is the gas column density with T₀ and N₀ as the reference gas temperature and column densities, respectively, at the reference radius, r₀. Similar to the power-law disk models in Sect. 4.2.3, we have used T₀ = 1500 K at r₀ of 0.1 AU. The images are then convolved with the clean beams in Table 1. We have assumed i = 55° and PA = 65°. For such orientation, we find that a Keplerian structure within the inner 100 AU is needed, surrounded by a rotating infalling envelope. The radius of the Keplerian structure is similar to the radius derived from previous methods (~100 AU), however, a stellar mass of 0.8 ± 0.3  M_⊙ is needed to reproduce the observed molecular lines. This is consistent with Yen et al. (2013) calculated at i = 30°, but a factor of 1.5 higher than the value obtained earlier (although the two values agree within their respective uncertainties).

4.3.2. TMC1

The velocity profile for TMC1 projected along PA = 90° is shown in Fig. 11 for ¹³CO. The moment maps toward this source indicate the presence of a compact flattened structure around the continuum position. Similar to TMC1A, the PV diagram also indicates a rotationally dominated structure. However, in contrast to TMC1A, the ¹³CO line is best described by an r^-1 profile indicating that it traces infalling rotating gas, while the C¹⁸O line is fitted well with a r^-0.5 slope as shown in Fig. 12 indicative of an RSD. Such discrepancy can occur if the two lines are tracing different regions. The moment zero maps in Fig. 5 already gave some indication that the ¹³CO and C¹⁸O emissions do not arise from similar structure. Furthermore, the disk mass toward TMC1 as inferred through the continuum is the lowest in our sample, hence it is also possible that the embedded RSD is not massive enough to significantly contribute to both the ¹³CO and C¹⁸O emission. The relatively more optically thick ¹³CO traces the inner envelope and thus is dominated by its kinematics, while the rotationally supported structure significantly contributes to the C¹⁸O line. A smooth extrapolation of the envelope down to 100–150 AU scales predicts C¹⁸O emission <18% of the observe emission at >1.5 km s^-1, suggesting that the majority of the emission at these velocities originates in a disk (see also models by Harsono et al. 2013, Fig. 8). Hence, using the C¹⁸O line, the fit indicates the presence of an RSD extending between 100−150 AU in radius. The best-fit stellar mass is 0.54$\begin{array}{c}\mathrm{+}\mathrm{0.2}\\ -0.1\end{array}$M_⊙ at i = 55° with typical error of ±0.2 km s^-1 on the determination of υ_lsr. The stellar mass is a factor of 1.5 lower than that previously derived by Hogerheijde et al. (1998) as indicated by the black lines in Fig. 9. However, a stellar mass of 0.4 M_⊙ also reasonably fits the ¹³CO PV diagram. In addition, the discrepancy between the ¹³CO and C¹⁸O velocity profiles may suggest that the most stable RSD is less than the C¹⁸O break point. Thus, we will use 100 AU as the size of the RSD toward TMC1.

4.3.3. TMR1

The determination of the velocity profile toward TMR1 is difficult. Adopting υ_lsr =  + 6.3 km s^-1 (Kristensen et al. 2012; Yıldız et al. 2013), the red-shifted peaks do not coincide with the blue-shifted peaks as expected from a symmetric rotating infalling structure in the absence of foreground clouds. The ¹³CO velocity profile is uncertain with a slope between 0.5 to 1 (Fig. 11). However, a velocity gradient is present in the PV analysis along and perpendicular to the outflow direction, which suggest outflow contamination toward TMR1. With the 5σ threshold used for this analysis, only 3 peaks are available in the C¹⁸O line, which are located along the blue outflow lobe as inferred from the ¹²CO J = 3–2 and J = 6–5 maps in Yıldız et al. (in prep.). The on-source C¹⁸O line is detected within the WideX correlator, which is double peaked within 100 AU diameter (Fig. 5). Thus, the non-detection of the turn-over from r^-1 to r^-0.5 indicates that the disk is much smaller than the synthesized beam (<0.7″ = 98 AU diameter). Also, the RSD has to be less massive than TMC1 and TMC1A due to the lack of on-source C¹⁸O emission and rotational motion. Thus, the compact flattened dust structure at >50 kλ toward this source cannot be fully associated with an RSD. The molecular emission from TMR1 must be mostly due to the infalling rotating envelope and the entrained outflow gas similar to the conclusion of Jørgensen et al. (2009).

4.3.4. L1536

Figure 11 shows the uv-space PV diagram with systemic velocity of υ_lsr =  +5.3 km s^-1. The peak positions were shifted to the center of the two peaks seen in the continuum (Fig. 2) and the C¹⁸O integrated map (Fig. 5). The dominant motion at low velocity offsets (Δυ =  ±2 km s^-1) is infall, while a flatter slope with a power-law close to r^-0.5 best describes the high velocity components. If the peak positions were not shifted, only the blue-shifted peaks follow the r^-0.5 profile starting from 70–80 AU. The flat velocity profile is based on 4 points and it does not show a very clear Keplerian disk profile such as seen toward TMC1A. From the high velocity peaks, the velocity profile seems to indicate a break at ~80 AU.

4.4. Summary and stellar masses

Three out of the four Class I sources show clear indications of embedded RSDs. From the turn-over radius (from r^-1 to r^-0.5), the radii of the Keplerian disks must be between 100–150 AU in TMC1, 80 AU toward L1536 and 100 AU toward TMC1A. Considering that the turnover radius is the outer radius of the embedded Keplerian disk, the stellar masses are derived from that radius and tabulated in Table 5. In the case of L1536, we chose the red-shifted peak (~3 km s^-1) at 70 AU as the turnover radius to determine the combined stellar masses. Thus, by obtaining the masses of the different components (envelope, disk, and star), it is now possible to place the sources in an evolutionary stage for comparison with theoretical models. The three sources are indeed in Stage I of the embedded phase, where the star has accreted most of its final mass.

The uv-space PV diagram assumes that the disk rotation is perpendicular to the outflow direction as inferred from large ¹²CO maps, foreground clouds do not contaminate the molecular lines, and the source systemic velocity is well determined. Small systemic velocity shifts can indeed shift both the red-shifted and blue-shifted velocity profiles from r^-0.5 to a steeper slope. However, for the firm detection of an embedded RSD toward TMC1A, this does not significantly affect the break radius. Further observations at higher resolution and sensitivity are certainly needed to confirm the extent of the Keplerian disks toward all of these sources. However, we argue that even with the current data it is possible to differentiate the RSD from the infalling rotating envelope with spatially and spectrally resolved molecular line observations.

4.5. Disk gas masses

After determining the extent of the embedded Keplerian disks, their gas masses can be estimated from the C¹⁸O integrated intensities. The J=2–1 line is expected to be thermalized at the typical densities (n_H2 > 10⁷ cm^-3) of the inner envelope due to the low critical density. The gas mass assuming no significant freeze-out is then given by the following equation (Scoville et al. 1986; Hogerheijde et al. 1998; Momose et al. 1998): ${\mathit{M}}_{\mathrm{gas}}\mathrm{=}\mathrm{5.45}\mathrm{\times }{\mathrm{10}}^{-4}\frac{{\mathit{T}}_{\mathrm{ex}}\mathrm{+}\mathrm{0.93}}{\exp \mathrm{(}\mathrm{-}{\mathit{E}}_{\mathrm{u}}\mathit{/}\mathit{k}{\mathit{T}}_{\mathrm{ex}}\mathrm{)}}\frac{\mathit{\tau }}{\mathrm{1}\mathrm{-}{\exp }^{\mathrm{-}\mathit{\tau }}}\mathrm{\int }{\mathit{S}}_{\mathit{\upsilon }}\mathrm{d}\mathit{\upsilon }{\mathit{M}}_{\mathrm{\odot }}\mathit{,}$(14)where h and k are natural constants, T_ex is the excitation temperature, τ is the line optical depth, E_u is the upper energy level in K and ${}^{\mathrm{\int }}\mathit{S}_{\mathit{\nu }}\mathrm{d}\mathit{\upsilon }$ is the integrated flux densities in Jy km s^-1 using a distance of 140 pc. The gas mass estimates of the disks inferred by integrating over a region similar to the extent of R_K in Table 5 are listed in Table 6 with τ = 0.5 and T_ex = 40 K. The adopted excitation temperature comes from the expected C¹⁸O rotational temperature within a 1″ beam if a disk dominates such emission (Harsono et al. 2013). From the observed $\frac{{\mathrm{13}}^{}\mathrm{CO}}{{\mathrm{C}}^{\mathrm{18}}\mathrm{O}}$ flux ratios, the ${\mathit{\tau }}_{{\mathrm{13}}^{}\mathrm{CO}}$ is estimated to be ≤4.0 with isotopic ratios of ¹²C/¹³C = 70 and ¹⁶O/¹⁸O = 550 (Wilson & Rood 1994). This implies that the C¹⁸O is almost optically thin (τ_C¹⁸O ≤ 0.5) and only a small correction is needed under the assumption that ¹³CO and C¹⁸O trace the same structure. Uncertainties due to the combined T_ex and τ result in a factor of ≤2 uncertainty in the derived disk masses. The gas disk masses toward TMC1A, TMC1 and TMR1 are at least a factor of 2 higher than the mass derived from the continuum data (assuming gas to dust ratio of 100). Since the velocity gradient is prominent in the moment one maps, there must be enough C¹⁸O in the gas phase within the RSD in order to significantly contribute to the lines. This is consistent with the evolutionary models of Visser et al. (2009) and Harsono et al. (2013) where only a low fraction of CO is frozen out within the embedded disks unless R_K > 100 AU due to the higher luminosities. Only toward L1536, which is the most evolved source in our sample, is the gas mass a factor of 5 lower than the mass derived from the continuum. This would be consistent with CO starting to freeze-out within the disk near the end of the main accretion phase (Visser et al. 2009). Thus, the small differences between gas and dust masses indicate relatively low CO freeze-out in embedded disks.

5. Discussion

5.1. Disk structure comparison: dust versus gas

In this paper, we have determined the sizes of the disk-like structures around Class I embedded YSOs using both the continuum and the gas lines. The continuum analysis focuses on the extent of continuum emission excluding the large-scale emission, while the gas line analysis focuses on placing constraint on the size of the RSD. The continuum analysis utilizes both intensity based disk models and a simple disk structure (power-law disk model). The intensity based disk models give a large spread of disk radii toward all sources. On the other other hand, the best-fit radii of power-law disks are larger than the Keplerian disk toward TMR1 and L1536, while the values are comparable toward TMC1A and TMC1. Not all of the disk structure derived from continuum observations may be associated to an RSD.

The main caveat in deriving the flattened disk structure from continuum visibilities is the estimate of the large-scale envelope emission. It has been previously shown that the large-scale envelope can deviate from spherical symmetry (e.g., Looney et al. 2007; Tobin et al. 2010). Disk structures may change if one adopts a 2D flattened envelope structure due to the mass distribution at 50–300 AU scale. However, in this paper, we focused on the size of the flattened structure where any deviation from the spherically symmetry model gives an estimate of its size. For the purpose of this paper, a spherical envelope model is used to estimate the large-scale envelope contribution since it is simple, fits the observed visibilities at short baselines, and does not require additional parameters.

5.1.1. Comparison to disk sizes and masses from SED modelling

A number of 2D envelope and disk models have been published previously constrained solely using continuum data. More recently, Eisner (2012) used the combination of high resolution near-IR images, the SED and millimeter continuum visibilities to constrain the structures around TMR1, TMC1A and L1536. The size of his disk is defined by the centrifugal radius, which is the radius at which the material distribution becomes more flattened (Ulrich 1976). His inferred sizes are consistent with the extent of RSDs in our sample (TMR1: 30–450 AU; L1536: 30–100 AU; TMC1A: 100 AU), but our higher resolution data narrow this range. For the case of TMR1, we can definitely rule out any disk sizes >100 AU.

Others have used a similar definition of the centrifugal radius using the envelope structure given by Ulrich (1976) with and without disk component (e.g., Eisner et al. 2005; Gramajo et al. 2007; Furlan et al. 2008). Thus, the value of R_c gives an approximate outer radius of the RSD. These models are then fitted to the SED and high resolution near-IR images without interferometric data. In general, the values of R_c derived from such models are lower (at least by factor of 2) than the extent of the RSDs as indicated by our molecular line observations. By trying to fit the models to the near-IR images which are dominated by the hot dust within the disk, they put more weight to the inner region. Thus, we argue that spatially resolved millimeter data are necessary in order to place better constraints on the compact flattened structure.

The disk masses derived in our work from both 50 kλ flux and uv modelling are consistent with each other. These values are typically higher by a factor of 2 and a factor of 8 for TMC1A³ than the disk masses reported by Eisner (2012) and Gramajo et al. (2010, only toward TMR1). The discrepancy may be due to the difference in the combination of the assumed inclination, the adopted envelope structure, and definition of the flattened structure.

5.2. Constraints on disk formation

A number of RSDs toward other embedded YSOs have been reported in the literature, which are summarized in Table 7. The values are obtained from analysis that included interferometric molecular line observations. Keplerian radii, R_K, between 80 and 300 AU are found toward these sources with disk masses in the range 0.004–0.055 M_⊙ around 0.2–2.5 M_⊙ stars. Most of the observed RSDs are found toward Stage I sources (M_⋆ > M_env), while RSDs are claimed toward three Stage 0/I sources. From the sample, large (≥100 AU) RSDs are generally found toward sources where M_⋆/M_tot > 0.65 with M_tot = M_⋆ + M_env + M_disk.

How do these disks compare to disk formation models? Semi-analytical disk formation models by Visser et al. (2009), which are based on Terebey et al. (1984) and Cassen & Moosman (1981), predict that >100 AU RSDs are expected by the end of the accretion phase (M_env < M_disk, see also Young & Evans 2005). The disk sizes depend on the initial envelope temperature or sound speed and the rotation rate. The evolution of R_K within these models is shown in Fig. 13 for the same initial conditions as Harsono et al. (2013) with comparable initial rotation rates to those measured by Goodman et al. (1993). Remarkably, the observed RSDs in Stage I fall within the expected sizes from these semi-analytical models. On the other hand, most of the Class 0 disks fall outside of these models. From this comparison, the observed RSDs in Stage I can be well represented by models with Ω = 10^-14 Hz and c_s = 0.19−0.26 km s^-1.

In semi-analytical models, the disk forms as a consequence of angular momentum conservation. However, recent numerical simulations of collapse and disk formation with the presence of magnetic fields do not always form a large RSD (>100 AU). Galli et al. (2006) showed that an RSD does not form in the ideal magnetohydrodynamics (MHD) case due to efficient magnetic breaking (e.g., Galli & Shu 1993; Mellon & Li 2008). The problem can be alleviated if the magnetic field is weak or with the inclusion of non-ideal MHD, although the results with non-ideal MHD are not entirely clear (Mellon & Li 2008, 2009; Dapp & Basu 2010; Krasnopolsky et al. 2011; Li et al. 2011; Braiding & Wardle 2012). Once the RSD forms, it can easily evolve to >400 AU (Dapp et al. 2012; Vorobyov 2010; Joos et al. 2012). The question centers on whether RSDs are already present in the Class 0 phase. This depends on the strength and orientation of the magnetic field, as well as the initial rotation rate of the envelope.

5.2.1. On the magnetic field

The dynamical importance of the magnetic field can be numerically determined through the value of the effective mass-to-flux ratio, λ_eff, which is the ratio of the gravitational and magnetic energies within the cloud. A value of 1 means that the cloud is strongly magnetized while value of >10 characterizes weakly magnetized clouds. Troland & Crutcher (2008) presented a recent survey of magnetic field strengths toward molecular cloud cores. They find an average λ_eff ~ 2–3 or an average magnetic field strength of 16 μG, which argues for the importance of magnetic field during the dynamical evolution of pre- and protostellar cores (see also Crutcher 2012).

The magnetic field strength within Taurus is found to vary toward different cores. For example, TMC1A is within the L1534 core whose line of sight magnetic field is weak (<1.2  μG, Crutcher et al. 2010), while the TMC1 dark cloud has considerable magnetic field strength of the order of 9.1 ± 2.2  μG. Since λ is inversely proportional to the magnetic field strength, the λ_eff toward these two sources are >10 and ~4 around TMC1A and TMC1, respectively. It is interesting to note that a large and massive RSD (R = 100 AU and M_disk/M_⋆ ~ 0.1) is found toward TMC1A, while the RSDs toward both TMC1 and TMR1, which are located closer to the TMC1 dark cloud, are small and less massive. Considering that M_disk/M_env is lower toward TMC1 and TMR1 than toward TMC1A, this also illustrates that the magnetic field may play an important role in determining the disk properties.

How likely are these disks formed out of weakly magnetized cores? Molecular cloud cores are found to have lower than average magnetic field strength (λ_eff > 3, Crutcher et al. 2010), than found on cloud scales. On the other hand, the fraction of cores with λ_eff > 5 is low. Consequently, it is unlikely that TMC1A and other observed RSDs in Stage I all formed out of weakly magnetized cores. Thus, simulations need to be able to form RSDs out of a moderately magnetized core that evolve up to sizes of 100–200 AU by Stage I. In addition, we can argue that the magnetic breaking efficiency drops as the envelope dissipates, which allows the rapid growth of RSDs near the end of the main accretion phase as shown in Fig 13. On the other hand, there may be a low fraction of cores which are weakly magnetized, which facilitate formation of large disks as early as Stage 0 phase (Krumholz et al. 2013; Murillo et al. 2013).

5.2.2. On the large-scale rotation

Another important parameter is the large-scale rotation around these sources, which is a parameter used in both simulations and semi-analytical models. Generally, a cloud rotation rate between 10^-14 Hz to a few 10^-13 Hz is used in numerical simulations (e.g., Yorke & Bodenheimer 1999; Li et al. 2011). The average rotation rate measured by Goodman et al. (1993) is ~10^-14 Hz, which is consistent with the initial rotation rate of the semi-analytical model that best describes a large fraction of the observed RSDs. The time at which Stage I starts within that model is 3 × 10⁵ years. For comparison, around the same computational time, Yorke & Bodenheimer (1999) form a > 500 AU RSD for their low-mass case (J). Similarly, Vorobyov (2010) also found >200 AU RSDs close to the start of the Class II phase. We find no evidence of such large RSDs at similar evolutionary stages in our observations.

Another aspect of the large-scale rotation that is worth investigating is its direction with respect to that of the disk rotation. Models assume that the disk rotation direction is the same as the rotation of its core. To compare the directions, we have used the velocity gradients reported by Caselli et al. (2002) and Goodman et al. (1993) toward the dark cores in Taurus. Velocity gradients are detected toward L1534 (TMC1A), L1536 and TMC-1C, which is located to the North East of TMR1 and TMC1. Table 7 lists the velocity gradients of the RSDs toward the embedded objects and their envelope/core rotation⁴. It is interesting to note that the large-scale rotation detected around or nearby the RSDs have a different velocity gradient direction. Similar misalignments in rotation have been pointed out by Takakuwa et al. (2012) toward L1551 NE and by Tobin et al. (2013) for L1527. Although the positions of the detected molecular emission in Caselli et al. (2002) do not necessarily coincide with the embedded objects due to chemical effects, the systematic velocities of the cores are similar. Misalignment at 1000 AU scales was also concluded by Brinch et al. (2007) toward L1489 where they found that the Keplerian disk and the rotating envelope are at an angle of 30° with respect to each other. If these objects formed out of these cores in the distant past, then either the collapse process changed the angular momentum of the envelope and, consequently, the disk rotation direction or it affects the angular momentum distribution of nearby cores through feedback. Interestingly, the former may suggest the importance of the non-ideal MHD effect called the Hall effect during disk formation, which can produce counter-rotating disks (Li et al. 2011; Braiding & Wardle 2012).

Finally, the simulations and models above require the initial core to be rotating in order to form RSDs. However, by analyzing the angular momentum distribution of magnetized cores, Dzib et al. (2010) suggested that the NH₃ and N₂H⁺ measurements may overestimate the angular momentum by a factor of ~8−10. This poses further problems to disk formation in previous numerical studies of collapse and disk formation. On the other hand, the inclusion of the Hall effect does not require a rotating envelope core to form a rotating structure (e.g., Krasnopolsky et al. 2011; Li et al. 2011). Further studies are required to compare the expected observables from rotating and non-rotating models.

5.3. Disk stability

Self gravitating disks can be important during the early stages of star formation in regulating the accretion process onto the protostar, which can be variable (e.g., Vorobyov 2010). Dunham & Vorobyov (2012) show that such episodic accretion events can reproduce the observed population of YSOs. A self-gravitating disk is expected when M_disk/M_star ≥ 0.1 (e.g., Lodato & Rice 2004; Boley et al. 2006; Cossins et al. 2009). However, most of the observed embedded RSDs around Class I sources have disk masses such that M_disk/M_star ≤ 0.1. It is unlikely for such disks to become self-gravitating, but they may have been in the past. Interestingly, a few embedded disk seem to be self-gravitating since their M_disk/M_star ≥ 0.1 such as toward TMC1A, IRS63, L1527 and NGC1333 IRAS4A2 (See Table 7).

The low fraction of large and massive disks during the embedded phase suggests that disk instabilities may have taken place in the past and limited their current observed numbers (i.e. low fluxes at long baselines). In this scenario, the embedded disks will have a faster evolution and, thus, spend most of their lifetime with small radii. Alternatively, disks experiencing high infall rates from envelope to the disk as expected during the embedded phase may also have higher accretion rates (Harsono et al. 2011). Such a scenario argues for disk instabilities to be present even for M_disk/M_⋆ < 0.1. In both cases, the dominant motion of the compact flattened structure during the embedded phase is infall rather than rotation. There are indeed sources where higher sensitivity and spatial resolution observations such as with Atacama Large Millimeter/submillimeter Array (ALMA) can detect features associated with disk instabilities (e.g., Cossins et al. 2010; Forgan et al. 2012; Douglas et al. 2013).

6. Conclusions

We present spatially and spectrally resolved observations down to a radius of 56 AU around four Class I YSOs in Taurus. The C¹⁸O and ¹³CO 2–1 lines are used to differentiate the infalling rotating envelope from the rotationally supported disk. Analysis of the dust and gas lines were performed directly in uv space to avoid any artefacts introduced during the inversion process. The main results of this paper are:

• Dust disk sizes and masses can be derived from the continuumby using a power-law spherical envelope model to account forthe large-scale envelope emission. Intensity based disk models(Uniform and Gaussian) give similar disk sizes to each otherwhich are lower by typically 25%–90% than realistic disk mod-els (power-law disk). On the other hand, they give similar diskmasses to more realistic models, which are consistent with thedisk mass derived from the 50 kλ flux point. Inclusion of the envelopein the analysis is important to obtain reliable disk masses. Theobservationally derived masses indicate that there is still a signif-icant envelope present toward three of the four sources: TM1A,TMC1, and TMR1.

• Three of the four sources (TMC1, L1536, TMC1A) host embedded rotationally supported disks (RSDs) derived from line data. By fitting velocity profiles (υ ∝ r^-0.5 or ∝ r^-1.0) to the red- and blue-shifted peaks, the RSDs are found to have outer radii of ~80–100 AU. In addition, the large-scale structure (>100 AU) is dominated by the infalling rotating envelope (υ ∝ r^-1). The derived stellar masses of these sources are of order 0.4–0.8 M_⊙, consistent with previous values. Consequently, these objects are indeed Stage I young stellar objects with M_⋆/(M_env + M_disk + M_⋆) > 0.7.

• Disk radii derived from the power-law dust disk models toward TMC1 and TMC1A are consistent with sizes of their RSDs. However, the dust disk radii toward TMR1 and L1536 are not the same as the extent of their RSDs. Thus, we emphasize that spatially and spectrally resolved gas lines observations are required to study the nature of flattened structures toward embedded YSOs.

• Semi-analytical models with Ω = 10^-14 Hz and c_s = 0.19−0.26 km s^-1describe most of the observed RSDs in Stage I. The observed RSDs argue for inefficient magnetic breaking near the end of the main accretion phase (M_⋆/M_tot > 0.65). More theoretical studies are needed to understand how and when RSDs form under such conditions.

• Comparison between disk masses derived from the continuum and C¹⁸O integrated line intensities (Table 6) suggests that relatively little of CO is frozen out within the embedded disk.

The current constraints on disk formation rely on a small number of observed RSDs around Class I sources. Certainly one needs to constrain the size of RSDs near the end of the main accretion phase to understand how late Class II disks form. On the other hand, the observed RSDs in the Class I phase do not answer the question when and how RSDs form. Such constraints require the detection and characterization of RSDs during Class 0 phase, which are currently lacking. From the sample of sources in Yen et al. (2013) and our results, it is clear that the infalling rotating flattened structure is present at >100 AU. Thus, extending molecular line observations down to <50 AU radius with ALMA toward Class 0 YSOs will be crucial for differentiating between the different scenarios for disk formation.

Online material

Appendix A: Observational data

The observational log can be found in Table A.1 and the channel maps of all the observed lines toward the four sources are shown in Figs. A.1–A.6.

Appendix B: Large-scale structure

The parameters of the large-scale envelope structure used to estimate the continuum emission at short baselines are in Table B.1.

Appendix C: WideX data

Spectra toward all sources integrated within a 3″ box around their continuum positions are shown in Fig. C.1. The integrated flux density maps toward C¹⁸O, SO 5₆−4₅ and SO₂ 11₁₁₁−10₀₀₁₀ (221.96521 GHz) are shown in Fig. C.3. The SO emission toward all sources is much more extended compared to the C¹⁸O emission, while the SO₂ emission is compact.

Appendix D: L1536

Figure D.1 shows the ¹²CO J = 3–2 spectral map obtained with JCMT, where there is no indication of classical bipolar outflows are present.

---

Acknowledgments

We thank the anonymous referee for the constructive comments. We are grateful to the IRAM staff, in particular to Chin-Shin Chang, Tessel van der Laan and Jérôme Pety, for their help with the observations, reduction and modelling of the data. We thank Umut Yıldız for reducing the single-dish data and making them available. This work is supported by the Netherlands Research School for Astronomy (NOVA) and by the Space Research Organization Netherlands (SRON). Astrochemistry in Leiden is supported by the Netherlands Research School for Astronomy (NOVA), by a Spinoza grant and grant 614.001.008 from the Netherlands Organisation for Scientific Research (NWO), and by the European Community’s Seventh Framework Programme FP7/2007–2013 under grant agreement 238258 (LASSIE). This research was supported by a Lundbeck Foundation Group Leader Fellowship to Jes Jørgensen. Research at Centre for Star and Planet Formation is funded by the Danish National Research Foundation and the University of Copenhagen’s programme of excellence. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency.

References

All Tables

All Figures