A Hot and Massive Accretion Disk around the High-Mass Protostar IRAS 20126+4104
We present new spectral line observations of the molecule in the accretion disk around the massive protostar IRAS 20126+4104 with the Submillimeter Array that for the first time measure the disk density, temperature, and rotational velocity with sufficient resolution (, equivalent to ) to assess the gravitational stability of the disk through the Toomre- parameter. Our observations resolve the central region that shows steeper velocity gradients with increasing upper state energy, indicating an increase in the rotational velocity of the hotter gas nearer the star. Such spin-up motions are characteristics of an accretion flow in a rotationally supported disk. We compare the observed data with synthetic image cubes produced by three-dimensional radiative transfer models describing a thin flared disk in Keplerian motion enveloped within the centrifugal radius of an angular-momentum-conserving accretion flow. Given a luminosity of , the optimized model gives a disk mass of and a radius of rotating about a protostar with a disk mass accretion rate of . Our study finds that, in contrast to some theoretical expectations, the disk is hot and stable to fragmentation with at all radii which permits a smooth accretion flow. These results put forward the first constraints on gravitational instabilities in massive protostellar disks, which are closely connected to the formation of companion stars and planetary systems by fragmentation.
What role accretion disks play in the formation of high-mass stars () remains a long-standing question. Circumstellar disks form naturally in the centers of rotating inflows and are a key element in the standard paradigm of the formation of Sun-like stars, providing for the growth of planetary systems. However, it remains debatable whether high-mass stars form in a similar fashion. The accretion rates in massive star formation may be high enough to induce gravitational instabilities and put disk-mediated accretion in doubt. Previous observations reporting disk-like accretion flows around high-mass protostars have not yet assessed the stability of the candidate disks (Chini et al. 2004; Patel et al. 2005; Jiménez-Serra et al. 2007; Cesaroni et al. 2014; Johnston et al. 2015). On theoretical grounds, a very high accretion rate () is required to form stars more massive than . Stars this massive undergo rapid enough Kelvin-Helmholtz contraction that they begin hydrogen burning while still accreting (Palla & Stahler 1993). A continuous resupply of fresh hydrogen is required to allow a growing massive protostar to reach the mass of a B or O star before exhausting its hydrogen fuel and leaving the main sequence (Keto 2003). The gas densities required for such rapid accretion could induce gravitational instabilities in the accretion flow that would make the disk prone to fragmentation, perhaps producing companion objects (Kratter & Matzner 2006). On the other hand, the stability of a disk is determined not only by the gas density, but also the differential shear and the gas temperature which are much affected by the disk environment (Durisen et al. 2007). High stellar luminosity () and shock heating may keep the disk warm enough to allow accretion to proceed steadily (Pickett et al. 2000b; Kratter et al. 2010). Massive disks are also subject to dynamical heating from spiral shocks, and vertical shear at the interface with the envelope (Pickett et al. 2000a; Harsono et al. 2011). The stability of the disk may be assessed by the Toomre- parameter,
which compares the stabilizing effects of the temperature (sound speed ) and shear (angular velocity ) against the clumping tendency toward instability induced by the surface density (). Values of imply a stable disk. Up until now, there has been insufficient observational evidence to answer the question. We present new molecular line observations of the massive protostar IRAS 20126+4104 that directly address the question of the gravitational stability of massive protostellar disks.
IRAS 20126+4104 (hereafter I20126) is a nearby (; Moscadelli et al. 2011), luminous (; Johnston et al. 2011) high-mass protostar with a well collimated bipolar molecular outflow (Zhang et al. 1999; Moscadelli et al. 2005; Su et al. 2007; Hofner et al. 2007) and an accretion disk (Zhang et al. 1998; Cesaroni et al. 2005, 2014; Xu et al. 2012). Previous radio and infrared observations suggest a disk mass of several and a protostellar mass of to (Sridharan et al. 2005; Cesaroni et al. 2005; Keto & Zhang 2010; Johnston et al. 2011). With such a high ratio of the disk to stellar mass, the self-gravity of the disk is significant in the accretion dynamics, and the stability of the disk is questionable. A previous study compared observed infrared emission against that predicted by a model of a disk-mediated accretion flow (Johnston et al. 2011) and found the Toomre- high enough to allow self-consistency with their original assumption of a disk. Our new observations directly measure the temperature, angular velocity, and surface density to enable a direct determination of the Toomre- and the disk stability.
Our observations of the line emissions around I20126 with the Submillimeter Array444The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics, and is funded by the Smithsonian Institution and the Academia Sinica. (SMA; Ho et al. 2004) achieved the highest possible angular resolution (, equivalent to ) sufficient to spatially resolve the accretion flow. The molecule is ideal for the identification of the hot, dense gas in the flow near the star. It requires high gas densities, , for collisional excitation, and high temperatures, , to produce detectable emission (Araya et al. 2005; Chen et al. 2006). The Doppler shifting of the molecular lines measures the rotational velocities and shear within the accretion disk. The emission includes multiple line transitions with different excitation energies allowing a measurement of the gas density from the line brightness and a measurement of the gas temperature from the brightness ratios.
The present study of the high-mass protostar I20126 finds its massive disk hot and stable against gravitational fragmentation even as the accretion proceeds at a high rate. This paper is organized as fellows. Sect. 2 describes details of the observations and data reduction. We present the basic observational results in Sect. 3 and explain details of radiative transfer models for both continuum and line emissions in Sect. 4. We discuss the implication of our results in Sect. 5 and conclude with a brief summary in Sect. 6.
2.1 SMA Observations
We carried out observations with the SMA on 2006 July 7 and 9 in the very extended configration and 2007 August 7 in the extended configuration. The phase tracking center is at (20:14:26.024,+41:13:32.76). The lower sideband covers the frequency range of while the upper sideband covers . The projected baseline is in the range of , which is insensitive to structures larger than (equivalent to ). The system temperature varied from to for the very extended configuration tracks and from to for the extended configuration track. Data inspection as well as bandpass and flux calibrations were done within the IDL superset MIR. The flux scale was derived from observations of the Jovian moon, Callisto, and is estimated to be accurate within 15%. Temporal gains were derived from the calibrator MWC 349 and imaging were performed with the MIRIAD package. Both sidebands were used to generate line-free continuum maps in a multi-frequency synthesis with an effective central frequency of . The continuum image reach an rms of with a beam size of and a position angle (P.A.) of using uniform weighting. The continuum flux density of I20126 is . The spectral data were gridded with a velocity resolution of at . Imaging with the robust parameter equal to zero yields a beam size of () and an rms of . The pixel size is set to be , which yields 68.0 and 90.8 pixels within one beam for the continuum and spectral line images, respectively.
2.2 BIMA Observations
In order to find the spectral index, , of the dust opacity in the disk, we measure the flux density from a continuum image produced by visibilities within the same - range as the SMA observations. The data were obtained with the Berkeley-Illnois-Maryland Association (BIMA) Millimeter Array on January 16 and February 14, 2002 in the A and B configuration, respectively. The system temperature varied from to for the A configuration track and from to for the B configuration track. The calibration and imaging were performed with the MIRIAD software package. Temporal gains and flux density were determined by observation on the calibrator MWC 349. Both sidebands were used to generate the continuum map in a multi-frequency synthesis with an effective central frequency of . The continuum flux density of I20126 is , which is consistent with the flux density reported by Cesaroni et al. (2014). Comparing the flux density between and , we obtain a continuum spectral index of , corresponding to a grain opacity spectral index of in the disk.
3 Observational Results
The dust emission (Fig. 1) shows a rather symmetric distribution with the peak position at (20:14:26.029,+41:13:32.57), which is assumed to be the position of the protostar. Such symmetric morphology already suggests a disk at a moderate inclination, but the character of a rotating disk is revealed through spectral imaging. The observed spectrum covers the components of the transition with upper state energy in the range of and the () emission with (Fig. 2). The lines are collisionally excited and thus sensitive to gas temperature (Araya et al. 2005; Chen et al. 2006). We can compare the Doppler shifting of spectral lines of different excitation energies to identify the spin-up of the accretion flow and the temperature gradient within disk created by the hot protostar. The excitation is indicated by two quantum numbers, and , which specify the total angular momentum of the molecule and its projection along the principal axis. As examples, we showed the components of the transition (Fig. 3a–3d), whose upper states have increasing energies of and are expected to trace emission progressively close to the hot protostar. In addition, the () line shows a similar and slightly more extended emission, likely to trace part of the envelope (Fig. 3e).
The spectral line data are three-dimensional (3-D) with axes of position, position, and velocity and can be displayed in several ways. Plotting the average spectral line velocity (first moment) at each position shows a consistent gradient characteristic of rotation for all five observed lines (Fig. 3a–3e). We can select a single plane in the data cube oriented along this gradient, and plot the line intensity as a function of position and velocity (P-V diagram). These plots (Fig. 4a–4e) reveal progressively steeper velocity gradients with increasing , indicating an increase in the rotational velocity of the higher temperature gas nearer the star. Such spin-up motions are characteristics of an accretion flow that at least partially conserves angular momentum, the limit of which is a rotationally supported or Keplerian disk.
The velocity gradients in the spectral line data also show the inward or radial flow of accretion in the more circular and bluer pattern of the average velocity of the line (Fig. 3a) compared with line (Fig. 3d). This effect arises from self-absorption of the spectral line emission in a radial flow. Along a line of sight through the center, a radial flow splits the line emission into red and blue shifted components from the near and far sides of the core. If the line were optically thin, we would see a symmetrically split profile. However, absorption by colder gas in the outer part of the core selectively absorbs the red-shifted emission, which is closer to its own velocity, while the emission from the far side, which is blue-shifted to a dissimilar velocity, passes through. If the flow were purely radial, this would produce a bullseye pattern of velocities, bluest in the center. In Fig. 3a, we see a combination of this effect along with the rotation. In contrast, a purely rotational flow will show a simpler one-dimensional gradient across the image as shown in Fig. 3d. This indicates that the high temperature gas has spun-up so much that the rotational velocities dominate the average. In I20126, the spectral line velocities and brightnesses are consistent with a more radial flow in a cooler accreting envelope that spins up and flattens to a hot rotating disk as the accretion flow approaches the star.
4 Radiative Transfer Models for Continuum and Line Emissions
The velocity patterns provide a qualitative overview of the accretion flow feeding I20126. We can extract precise measurements of the disk and envelope temperatures, densities, and velocities by comparing the observations with model accretion flows. Following earlier studies (Keto & Zhang 2010; Johnston et al. 2011) that modeled lower angular resolution observations of other spectral lines, e.g. , and the infrared continuum observations of I20126, we constructed a 3-D analytical model describing a thin accretion disk (Pringle 1981) in Keplerian motion enveloped within the centrifugal radius of an angular-momentum-conserving accretion flow (Keto & Zhang 2010; Ulrich 1976). We include the stellar irradiation to heat the flared disk (Kenyon & Hartmann 1987) consistent with the presence of outflow cavities (Qiu et al. 2008; Moscadelli et al. 2011). A bolometric luminosity of is assumed (Johnston et al. 2011). We also adopt the systemic velocity of (Cesaroni et al. 1999).
Under conditions of local thermodynamic equilibrium (LTE), we solved the radiative transfer equation for the intensity, , of the continuum and spectral lines simultaneously to construct synthetic continuum images and 3-D spectral image cube for comparison with the observations. The radiative transfer is evaluated with the source function, i.e. the Planck function, , modulated by the linear sum of opacities for the continuum and spectral lines. At the center frequency of each individual channel, we solve for with
where is the total absorption coefficient given by
For the continuum emission, we have , where is the mass density, and a gas-to-dust mass ratio of is assumed. We use the dust opacity law (Hildebrand 1983) with for the disk (Sec. 2.2) and for the envelope (Johnston et al. 2011). Regarding spectral lines, the line blending is significant for a few molecular lines so we have
where is the total number of lines included in the model. The absorption coefficient of the -th line, , is given by
where is the gas density, the abundance of the molecule, the upper state degeneracy, the upper state energy, the partition function of the molecule, the spontaneous emission rate of the line, the rest frequency of the line, and the model line profile (see Appendix A).
The integration along line of sight is performed using the Runge-Kutta method in steps of . For comparison with the observed spectral line data, synthetic spectral line images are generated through observation simulation, including visibility sampling and image making, with MIRIAD. The Levenberg-Marquardt method was used to optimize models by minimizing value computed with both the continuum and spectral line data. Only channels with emission stronger than in the central region (dark gray histogram in Fig. 2) are considered for the model fitting and the computation. The model spectral imaging includes 19 molecular lines (Table 1): the components of the transition, the components of the transition, and (). The ratio is assumed to be 70 for a galactocentric distance of (Wilson & Rood 1994). The model contains nine adjustable parameters as listed in Table 2. The uncertainty of the parameter is estimated by , where is the -th diagonal term of the covariance matrix, is the number of pixels in the synthesized beam, and gives the 68.3% confidence level for degrees of freedom of nine.
Analogous to the standard model of star formation with disk-mediated accretion, our model consists of a thin disk (Pringle 1981) of mass in Keplerian motion around a stellar mass, , residing within the centrifugal radius, , of an accretion flow with constant specific angular momentum, , in an infalling envelope (Ulrich 1976). To account for a fairly large disk mass, the centrifugal radius is computed with
Density-weighted means are calculated for velocity and temperature in positions where both disk and envelope are present. Due to the prominent outflow cavities (Qiu et al. 2008), stellar irradiation is also included for anticipated heating due to the flared disk geometry (Kenyon & Hartmann 1987). Considering various heating processes, e.g. accretion shocks, that may occur to further raise the disk temperature, we introduce a scaling factor, , to take these effects into account. Due to a fairly large disk mass, corrections for the enclosed disk mass as a function of radius are applied to rotation velocities and viscous heating produced by differential shear in the disk.
We use to denote the envelope radius in spherical coordinates and to denote the disk radius in cylindrical coordinates. Given an envelope mass accretion rate, , we obtain the density distribution of the envelope to be
where is the initial polar angle of the streamline, and and are the polar radius and angle along a streamline. The envelope radius is assumed to be , larger than the structures sensitive with our observations. The velocity field is described by (Ulrich 1976; Mendoza et al. 2004)
where is the Keplerian velocity at . In the inner region , the mass of the disk gradually decreases to zero as one approaches the protostar. Eq. (8) is an approximation for the actual velocity field that responds to the mass distribution of the disk. Since our observed features are mainly attributed to the disk component, we just treat the inner envelope approximately.
The density distribution of the flared disk is described as
where is the disk scale-height with , and is related to the disk mass, , by , where is the inner radius of the disk. The surface density is defined by , and the enclosed disk mass is given by . The Keplerian rotation speed in the disk is hence
Given a disk mass accretion rate, , and the disk surface density, , an inward accretion velocity can be computed from
which is roughly constant through the disk.
We describe the central massive protostar as a zero-age main-sequence star (Schaller et al. 1992) with additional surface heating by gas accreted from the inner edge of the disk, releasing all its free-fall energy (Calvet & Gullbring 1998; Johnston et al. 2011). The very low X-ray luminosity of (Anderson et al. 2011) suggests that the free-fall energy is largely absorbed in the stellar surface. Hence, the accretion luminosity that goes into heating the stellar surface is given by
The model stellar luminosity, , including the luminosity of the protostar, , and the accretion heating, , is
A stellar luminosity of determined from infrared observations (Johnston et al. 2011) is applied to constrain the emerging flux at model stellar surface
where and is determined by the protostellar mass, , following relations of zero-age main-sequence stars.
To obtain the disk temperature distribution, we consider both the accretion heating and stellar irradiation, which is responsible for the vertical thermal gradient in the disk and heating in the outer part. First, we compute the disk temperature at one disk scale-height, , with
where is the temperature scaling factor, the Stefan-Boltzmann constant, the accretion flux (Pringle 1981), and the irradiation flux (Kenyon & Hartmann 1987). The disk temperature at small approaches dominated by accretion luminosity (Pringle 1981) while the outer part is mainly heated by stellar irradiation with (Kenyon & Hartmann 1987). We then obtain the vertical thermal gradient due to disk surface heating by stellar irradiation with
where is a parameter describing the increase in temperature with height from the mid-plane and is set to (Dartois et al. 2003). The temperature profile of the envelope follows the analytical scheme (Kenyon et al. 1993) that gives in the inner optically thick regime for infrared photons and in the optically thin outer part with (Johnston et al. 2011). The innermost dust-free zone has and the boundary is set by dust sublimation temperature of .
The total velocity dispersion of line broadening, , is computed by combining the turbulent broadening of velocity dispersion, , in quadrature with the thermal broadening, , where is the mass of Hydrogen atom. For the envelope, we adopted Larsonâs law to describe the turbulent broadening of velocity dispersion, (Larson 1981). Since this turbulent broadening is smaller than our spectral resolution of in most part of the envelope, we apply a line profile with instrumental broadening, (see Appendix A). In the disk, the Shakura-Sunyaev parameter (Shakura & Sunyaev 1973), is used to describe the kinematic viscosity, . In the thin-disk theory, the kinematic viscosity connects mass accretion rate to surface density through
Hence, we have the parameter described by
which just weakly depends on and with . For non-thermal broadening in the disk, we found the ratio of the turbulent to thermal broadening through the parameter with , where is the magnetic field. Here we assume the kinematic viscosity mainly attributed to turbulence. The magnetic term becomes important if the field strength is greater than a critical value of at . Given the typical field strength of in star-forming cores (Crutcher et al. 2010), the approximation of the parameter for deriving the turbulent broadening is reasonable in I20126. We kept in the disk to make turbulence subsonic when optimizing models.
5 Results and Discussion
The model is specified by nine adjustable parameters, listed in Table 2, that are optimized by minimization. Since our observations are insensitive to , we do not intend to determine its value but assume based on previous infrared studies (Johnston et al. 2011) and findings in disks around low-mass stars (Shu et al. 1994). The best-fit model gives a disk mass of and a centrifugal radius, , of rotating about a protostar with an accretion luminosity of , giving a disk mass accretion rate, , of . The optimization obtains a reduced value of 1.6, which is optimized with pixels including both the continuum and spectral line data, equivalent to independent data points. The derived parameters along with their uncertainties are listed in Table 3. The synthetic continuum image from the best-fit model is shown in Fig. 1b, while the integrated intensity images and the P-V diagrams of the and emissions are shown in Fig. 3f–3j and Fig. 4f–4j, respectively. For choices of in the range of , parameters of the best-fit models vary slightly within 4% of those listed in Table 2 with the exception of , which compensates the change of and varies within 30%.
|Inclination angle††Inclination angle is defined as the angle between the disk axis and the line of sight.,|
|Disk temperature scaling factor,|
|Specific angular momentum,|
|Envelope accretion rate,|
The accretion timescale is estimated by
which gives the time for a mass element in the disk to be accreted, and . The rotation period of a Keplerian disk is
Hence the ratio and reaches a minimum at . As long as at , the accretion timescale is longer than the rotation period at all radii, and the accretion is mediated by a rotationally supported disk. In I20126, the surface density at is , which implies an inward velocity, , of . The accretion timescale at , , is , longer than the rotation period of , so the accretion is mediated by a rotationally supported disk.
As the stellar luminosity is actually a proxy for the stellar mass, our model derives a larger stellar mass, leading to a moderately inclined disk at an angle of rather than an edge-on geometry, i.e. (Cesaroni et al. 2005, 2014). A smaller inclination angle of has also been suggested by numerical simulations with a misaligned magnetic field with respect to the disk rotation axis (Shinnaga et al. 2012). The temperature of the disk is above at all radii and becomes warmer than the enveloping accretion flow beyond . Figure 5 shows the temperature distributions of the disk and the envelope. The mid-plane density of the disk is everywhere larger than and is significantly higher than the critical density for collisional de-excitation for thermalization of the lines, which is about at . The densities in the envelope are lower and thus may not fully satisfy the LTE conditions. Since the observed lines are dominated by the disk component, the approximate treatment of the envelope is not significant.
The derived gas density and mass depend on the molecular abundances, which are assumed to be constant in the entire system, including the disk and envelope. While variations of the molecular abundance cannot be ruled out, the fact that the derived temperature exceeds ice mantle sublimation point readily implies the enhancement of in the gas phase. The derived abundance of is comparable with those found in hot molecular cores (Hernández-Hernández et al. 2014). An abrupt jump of abundance within the domain of interests is therefore not expected.
|Disk accretion rate,||(12)|
|Inward accretion velocity,||(11)|
|Representative value at|
Having determined the properties of the disk by comparison with our observations, we can assess its dynamical stability. We calculate the Toomre- parameter, , and find it larger than everywhere in the disk, which makes the disk stable to fragmentation (Fig. 5). The fractional uncertainty of Toomre- is about 27% through the disk. The Shakura-Sunyaev parameter (Shakura & Sunyaev 1973) is within and decreases slowly to at . The value of is much higher than the typical value of , above which fragmentation occurs in simulations of isolated disks with simple cooling laws (Rice et al. 2005). Yet numerical models appropriate for high-mass protostars tend to find higher (Kratter et al. 2010). A direct comparison between observations and simulations may still be premature because the evolution of the accretion flow in I20126 may not necessarily follow the evolution assumed for the simulations (Kratter et al. 2008).
We present new continuum and spectral line observations of the disk around the high-mass protostar I20126 with the SMA that achieved the highest angular resolution (, equivalent to ) to resolve the accretion flow. The continuum emission shows a fairly symmetric morphology, which suggests the disk at a moderate inclination angle. Observations of and lines resolve the central region, where the kinematics display a clear rotation velocity pattern. Position-velocity diagrams of the lines reveal progressively steeper velocity gradients with increasing upper state energy, indicating an increase in the rotational velocity of the hotter gas nearer the protostar. Such spin-up motions are characteristics of a rotationally supported disk.
We assess the dynamical stability of this massive disk through the Toomre- parameter. To evaluate the Toomre- as a function of radius through the disk, we measure the gas density, temperature, and rotational velocity in the disk by comparing data with synthetic data generated by radiative transfer models analogous to the standard model of star formation with disk-mediated accretion. Given a luminosity of , the optimized model finds a disk mass of and a centrifugal radius of rotating about a protostar with a disk mass accretion rate of . These physical conditions render everywhere in the disk, which makes the disk stable to fragmentation.
Our high angular resolution SMA observations of I20126 provide evidence for a stable massive accretion disk around a high-mass protostar. In contrast to some theoretical expectations of massive disks prone to local instabilities, the disk of I20126 is found to be hot and stable to fragmentation even as the accretion proceeds at a high rate. Such conditions may help to maintain the disk around massive stars and preserve opportunities for developing companions or a planetary system in a later phase of the protostellar evolution.
Appendix A Instrumental Broadening Using a Boxcar Function
Since the turbulence line width in the inner part of the envelope is smaller than the spectral resolution of our observations, it is necessary to account for the instrumental broadening in our models. The model line profile, , is given by the convolution of a boxcar, , set by channel spectral resolution, , and a gaussian function, , of line broadening, , set by thermal broadening, , and nonthermal broadening, , with . The intrinsic gaussian line profile is given by
where is the systemic velocity. The frequency response of one channel is approximated by a boxcar function
Hence, we calculate the line profile with the instrumental broadening by convolving with
where is the error function. Since a simple formula is available for the complementary error function, , one can also rewrite the line profile as
Examples of are shown in Fig. 6 for the case of a spectral resolution of , same as our observations.
- Anderson et al. (2011) Anderson, C. N., Hofner, P., Shepherd, D., & Creech-Eakman, M. 2011, AJ, 142, 158
- Araya et al. (2005) Araya, E., Hofner, P., Kurtz, S., Bronfman, L., & DeDeo, S. 2005, ApJS, 157, 279
- Calvet & Gullbring (1998) Calvet, N., & Gullbring, E. 1998, ApJ, 509, 802
- Cesaroni et al. (1999) Cesaroni, R., Felli, M., Jenness, T., et al. 1999, A&A, 345, 949
- Cesaroni et al. (2014) Cesaroni, R., Galli, D., Neri, R., & Walmsley, C. M. 2014, A&A, 566, 73
- Cesaroni et al. (2005) Cesaroni, R., Neri, R., Olmi, L., et al. 2005, A&A, 434, 1039
- Chen et al. (2006) Chen, H.-R., Welch, W. J., Wilner, D. J., & Sutton, E. C. 2006, ApJ, 639, 975
- Chini et al. (2004) Chini, R. S., Hoffmeister, V., Kimeswenger, S., et al. 2004, Nature, 429, 155
- Crutcher et al. (2010) Crutcher, R. M., Wandelt, B., Heiles, C., Falgarone, E., & Troland, T. H. 2010, ApJ, 725, 466
- Dartois et al. (2003) Dartois, E., Dutrey, A., & Guilloteau, S. 2003, A&A, 399, 773
- Durisen et al. (2007) Durisen, R. H., Boss, A. P., Mayer, L., et al. 2007, Protostars and Planets V, 607
- Edris et al. (2005) Edris, K. A., Fuller, G. A., Cohen, R. J., & Etoka, S. 2005, A&A, 434, 213
- Harsono et al. (2011) Harsono, D., Alexander, R. D., & Levin, Y. 2011, MNRAS, -1, 141
- Hernández-Hernández et al. (2014) Hernández-Hernández, V., Zapata, L., Kurtz, S., & Garay, G. 2014, ApJ, 786, 38
- Hildebrand (1983) Hildebrand, R. H. 1983, QJRAS, 24, 267
- Ho et al. (2004) Ho, P. T. P., Moran, J. M., & Lo, K. Y. 2004, ApJ, 616, L1
- Hofner et al. (2007) Hofner, P., Cesaroni, R., Olmi, L., et al. 2007, A&A, 465, 197
- Jiménez-Serra et al. (2007) Jiménez-Serra, I., Martín-Pintado, J., Rodríguez-Franco, A., et al. 2007, ApJ, 661, L187
- Johnston et al. (2011) Johnston, K. G., Keto, E., Robitaille, T. P., & Wood, K. 2011, MNRAS, 415, 2953
- Johnston et al. (2015) Johnston, K. G., Robitaille, T. P., Beuther, H., et al. 2015, ApJ, 813, L19
- Kenyon et al. (1993) Kenyon, S. J., Calvet, N., & Hartmann, L. 1993, ApJ, 414, 676
- Kenyon & Hartmann (1987) Kenyon, S. J., & Hartmann, L. 1987, ApJ, 323, 714
- Keto (2003) Keto, E. 2003, ApJ, 599, 1196
- Keto & Zhang (2010) Keto, E., & Zhang, Q. 2010, MNRAS, 406, 102
- Kratter & Matzner (2006) Kratter, K. M., & Matzner, C. D. 2006, MNRAS, 373, 1563
- Kratter et al. (2008) Kratter, K. M., Matzner, C. D., & Krumholz, M. R. 2008, ApJ, 681, 375
- Kratter et al. (2010) Kratter, K. M., Matzner, C. D., Krumholz, M. R., & Klein, R. I. 2010, ApJ, 708, 1585
- Larson (1981) Larson, R. B. 1981, MNRAS, 194, 809
- Mendoza et al. (2004) Mendoza, S., Cantó, J., & Raga, A. C. 2004, Revista Mexicana de Astronomía y Astrofísica, 40, 147
- Moscadelli et al. (2005) Moscadelli, L., Cesaroni, R., & Rioja, M. J. 2005, A&A, 438, 889
- Moscadelli et al. (2011) Moscadelli, L., Cesaroni, R., Rioja, M. J., Dodson, R., & Reid, M. J. 2011, A&A, 526, 66
- Palla & Stahler (1993) Palla, F., & Stahler, S. W. 1993, ApJ, 418, 414
- Patel et al. (2005) Patel, N. A., Curiel, S., Sridharan, T. K., et al. 2005, Nature, 437, 109
- Pickett et al. (2000a) Pickett, B. K., Cassen, P., Durisen, R. H., & Link, R. 2000a, ApJ, 529, 1034
- Pickett et al. (2000b) Pickett, B. K., Durisen, R. H., Cassen, P., & Mejía, A. C. 2000b, ApJ, 540, L95
- Pringle (1981) Pringle, J. E. 1981, ARA&A, 19, 137
- Qiu et al. (2008) Qiu, K., Zhang, Q., Megeath, S. T., et al. 2008, ApJ, 685, 1005
- Rice et al. (2005) Rice, W. K. M., Lodato, G., & Armitage, P. J. 2005, MNRAS, 364, L56
- Schaller et al. (1992) Schaller, G., Schaerer, D., Meynet, G., & Maeder, A. 1992, A&AS, 96, 269
- Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
- Shinnaga et al. (2012) Shinnaga, H., Novak, G., Vaillancourt, J. E., et al. 2012, ApJ, 750, L29
- Shu et al. (1994) Shu, F. H., Najita, J., Ostriker, E., et al. 1994, ApJ, 429, 781
- Sridharan et al. (2005) Sridharan, T. K., Williams, S. J., & Fuller, G. A. 2005, ApJ, 631, L73
- Su et al. (2007) Su, Y.-N., Liu, S.-Y., Chen, H.-R., Zhang, Q., & Cesaroni, R. 2007, ApJ, 671, 571
- Trinidad et al. (2005) Trinidad, M. A., Curiel, S., Migenes, V., et al. 2005, AJ, 130, 2206
- Ulrich (1976) Ulrich, R. K. 1976, ApJ, 210, 377
- Wilson & Rood (1994) Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191
- Xu et al. (2012) Xu, J.-L., Wang, J.-J., & Ning, C.-C. 2012, ApJ, 744, 175
- Zhang et al. (1998) Zhang, Q., Hunter, T. R., & Sridharan, T. K. 1998, ApJ, 505, L151
- Zhang et al. (1999) Zhang, Q., Hunter, T. R., Sridharan, T. K., & Cesaroni, R. 1999, ApJ, 527, L117