EECC for Star-Forming Cloud Core L1517B

Starless Cloud Core L1517B in Envelope Expansion with Core Collapse

Tian-Ming Fu1 ,Yang Gao1 2 , Yu-Qing Lou1 3 4
1affiliation: Department of Physics and Tsinghua Center for Astrophysics (THCA), Tsinghua University, Beijing 100084, China.
2affiliation: Center for Combustion Energy and Department of Thermal Engineering, Tsinghua University, Beijing 100084, China.
3affiliation: Department of Astronomy and Astrophysics, the University of Chicago, 5460 South Ellis Avenue, Chicago, IL 60637, USA.
4affiliation: National Astronomical Observatories, Chinese Academy of Science, A20, Datun Road, Beijing 100012, China.

Various spectral emission lines from star-forming molecular cloud core L1517B manifest red asymmetric double-peaked profiles with stronger red peaks and weaker blue peaks, in contrast to the oft-observed blue-skewed molecular spectral line profiles with blue peaks stronger than red peaks. Invoking a spherically symmetric general polytropic hydrodynamic shock model for the envelope expansion with core collapse (EECC) phase, we show the radial flow velocity, mass density and temperature structures of self-similar evolution for L1517B in a dynamically consistent manner. By prescribing simple radial profiles of abundance distribution for pertinent molecules, we perform molecular excitation and radiative transfer calculations using the publicly available RATRAN code set for the spherically symmetric case. Emphatically, spectral profiles of line emissions from the same molecules but for different line transitions as well as spectra of closely pertinent isotopologues strongly constrain the self-similar hydrodynamics of a cloud core with prescribed abundances. Our computational results show that the EECC model reproduces molecular spectral line profiles in sensible agreement with observational data of Institut de Radioastronomie Millimétrique (IRAM), Five College Radio Astronomical Observatory (FCRAO) and Effelsberg 100 m telescopes for L1517B. We also report spatially resolved observations of optically thick line HCO using the Purple Mountain Observatory (PMO) 13.7 m telescope at Delingha in China and the relevant fitting results. Hyperfine line structures of NH and NH transitions are also fitted to consistently reveal the dynamics of central core collapse. As a consistent model check, radial profiles of 1.2 mm and 850 m dust continua observed by IRAM 30 m telescope and the Submillimeter Common-User Bolometer Array (SCUBA), respectively, are also fitted numerically using the same EECC model that produces the molecular line profiles. L1517B is likely undergoing an EECC shock phase. For future observational tests, we also predict several molecular line profiles with spatial distributions, radial profile of sub-millimeter continuum at wavelength 450m, as well as the radial profiles of the column density and visual extinction for L1517B.

Subject headings:
dust, extinction — hydrodynamics — ISM: clouds — line: profiles — radio continuum: ISM — shock waves
slugcomment: Received 2010…; accepted 2011…; published…

1. Introduction

Observations over past two decades provide effective diagnostics to examine dynamic properties of molecular clouds where star formations are taking place in relatively early phases. These observations, including molecular spectral line profiles (e.g. Benson & Myers, 1989; Tafalla et al., 2002; Aguti et al., 2007), radial profiles of (sub)millimeter radio continua (e.g. Motte, André & Neri, 1998; Shirley et al., 2000; Stanke et al., 2006; Nutter & Ward-Thompson, 2007), and infrared dust extinctions (e.g. Alves, Lada & Lada, 2001), can be utilized separately and/or in combination to reveal several key hydrodynamic cloud features, such as core collapses, envelope expansions, travelling shocks and turbulence in the early evolution of star-forming cloud cores.

On a relatively independent track, the theoretical model framework based on self-similar hydrodynamic collapses has achieved substantial development since late 1960s (e.g. Bodenheimer & Sweigart, 1968; Larson, 1969; Penston, 1969; Shu, 1977; Hunter, 1977), and different variations of these theoretical model applications have been advanced in recent years (e.g. Terebey, Shu & Cassen, 1984; Lou & Shen, 2004; Wang & Lou, 2008; Stahler & Yen, 2009). For grossly spherical molecular clouds, the gas flow velocity, mass density and thermal temperature structures of star-forming molecular cloud cores predicted by theoretical models can be tested by several complementary high-resolution observations of selected source candidates. Attempts of this kind have been pursued all along by comparing theoretical cloud model results with the observed molecular spectral line profiles (e.g. Zhou et al., 1993; Testi & Sargent, 1998; Tafalla et al., 2004, 2006), (sub)millimeter radio continuum emissions (e.g. Shirley et al., 2002; Motte & André, 2001; Harvey, Wilner & Myers, 2003), and infrared dust extinctions (e.g. Alves, Lada & Lada, 2001; Kandori et al., 2005).

Molecular spectral line profiles are important to study physical structures and kinematics for various phases in star formation processes (e.g. Dyson & Williams, 1997). They need to be physically interpreted and can be properly utilized to probe gas densities and temperatures through excitations of molecules, and detect cloud turbulence as well as systematic bulk flow motions through their line widths and Doppler shifts etc. Observations of spectral lines more frequently exhibit blue skewed double-peak molecular profiles with stronger blue peaks and weaker red peaks (i.e. the so-called blue profiles, see e.g. Gregersen et al., 1997, 2000; Lee, Myers & Tafalla, 1999). These blue asymmetric spectral line profiles of molecular transitions are usually interpreted as signatures of core collapses with a static outer envelope and a temperature variation involved (e.g. Zhou et al., 1993; Myers, Evans & Ohashi, 2000; Gao, Lou & Wu, 2009). However, more and more high-resolution observations reveal distinct red asymmetric double-peak profiles with weaker blue peaks and stronger red peaks (i.e. the so-called red profiles, see e.g. Mardones et al., 1997; Thompson & White, 2004; Fuller, Williams & Sridharan, 2005; Velusamy et al., 2008), which cannot be understood within the scenario of complete collapse models. The presence of ‘red profiles’ for molecular transitions most likely indicates a collapsing central core surrounded by an expanding outer envelope (Lou & Shen, 2004; Thompson & White, 2004; Gao & Lou, 2010; Lou & Gao, 2011), though there have been alternative origins proposed (see Gao & Lou 2010 and Lou & Gao 2011).

The theoretical framework of self-similar hydrodynamics for envelope expansion with core collapse (EECC) was first advanced by Lou & Shen (2004) for an isothermal gas cloud dynamics with spherical symmetry. This isothermal EECC model was later extended to a more general polytropic equation of state (EoS) with the specific entropy conservation along streamlines (Wang & Lou 2007; Wang & Lou, 2008).

We advance the following scenario for the possible emergence of EECC dynamic phase in star-forming molecular clouds. Both theoretical (e.g. Broderick et al., 2002; Stahler & Yen, 2010) and observational (e.g. Lee & Myers, 2011) works reveal that a static core would show expansive or oscillatory motions once being perturbed on large scales. While a molecular cloud is undergoing an expansive or oscillatory mode on large scales, its core may begin to collapse or contract due to nonlinear instabilities under self-gravity. With the enhancement of such central infalling dynamics induced by self-gravity, the collapsed or contracted region expands outward in a dynamic manner. Eventually, the entire cloud starts to collapse once the frontier of the collapsed region reaches the outer cloud edge. Thus, our EECC phase usually exists in the early stage of protostar formation, consistent with observations regarding the internal dynamic motions of starless cores (Lee & Myers, 2011).

Radiative transfer calculations based on such general polytropic EECC shock model have been systematically studied by Gao & Lou (2010) and by Lou & Gao (2011) for the source FeSt 1-457 to emphatically demonstrate that the widely observed ‘red profile’ signatures in molecular line spectral profiles from low-mass star-forming clouds could be sensibly viewed as a diagnostic evidence revealing self-similar EECC hydrodynamic shock processes therein. To test the EECC interpretation of ‘red profiles’, we investigate low-mass star-forming cloud cores with quasi-spherical configurations. In this paper, we focus on a candidate source of starless cloud core L1517B for comprehensive data fitting and comparison.

To achieve high spatial resolutions, star-forming clouds with relatively simple internal structures yet without apparent envelope disturbances are chosen from the nearby Taurus complex. After a survey of known cases for star-forming cloud cores with red-skewed molecular line profiles in the Taurus complex, we choose the cloud core L1517B for a more comprehensive investigation as it possibly involves global expansions. Central spectral line profile observations of molecular transitions HCO(), HCO(), HCO(), HCO(), CS() and CS() from this candidate source clearly show red asymmetric spectral line profiles (e.g. Tafalla et al., 2004, 2006), which likely reveal the existence of an EECC shock dynamic phase. Moreover, millimeter and sub-millimeter continuum mappings serve as independent constraints on radial density and temperature profiles of L1517B (e.g. Tafalla et al., 2004; Kirk, Ward-Thompson & André, 2005) for the same EECC shock model used for fitting molecular line profiles.

Our motivation is to demonstrate specifically that the theoretical explanation of red profiles based on general polytropic EECC solutions with collapses, expansions and shocks (Gao & Lou, 2010), appears grossly consistent with available observations for starless cloud core L1517B. To further constrain and verify the EECC shock model, we present the result of a 5 5 (with steps) spatially resolved HCO molecular line profile observation using the 13.7 m telescope of Purple Mountain Observatory (PMO) at Delingha in Qinghai Province of China. We also report the model fitting results of these spatially resolved spectra on the basis of the same self-similar EECC shock model and show that the EECC shock dynamic phase will influence molecular line profiles observed at different lines of sight (LOSs) from the cloud, rather than simply affects the central spectral line profiles alone. We also predict spatially resolved profiles for several other molecular emission lines, sub-millimeter continuum at 450, and the column density distribution of L1517B for future observational tests.

This paper is structured as follows. The background information and our motivation are introduced here in Section 1. In Section 2, we specify model parameters and physical properties for L1517B using the framework of general polytropic EECC hydrodynamic shock model (Wang & Lou, 2008) without magnetic fields. Model fittings of central spectra from several molecular line transitions as well as spatially resolved spectral line profiles for the HCO transition of L1517B are shown and analyzed in Section 3. We also venture to predict other spectral line profiles with spatial resolutions. In Section 4, we show millimeter and sub-millimeter continua as further proof-tests to the same EECC shock dynamic model. Column number density and visual extinctions are predicted in Section 4.3 and we conclude in Section 5 with speculations. Details of analysis are summarized in Appendices A and B for a convenient reference.

2. Properties of Starless Cloud Core L1517B

2.1. Earlier Studies on Cloud Core L1517B

The candidate source L1517B as a star-forming cloud core is located in the Taurus complex, at an estimated distance of pc from us (Elias, 1978). According to the radio continuum mappings and molecular line intensity mappings (Tafalla et al., 2002, 2006), cloud core L1517B has a radius of (corresponding to AU at a distance of  pc) and a center identified at the right ascension and the declination . Millimeter and sub-millimeter radio continuum data of L1517B have been acquired at 1.2 mm wavelength using the Institut de Radioastronomie Millimétrique (IRAM) 30m telescope by Tafalla et al. (2002, 2004), and at wavelengths of 850 m and 450 m using the Sub-millimeter Common-User Bolometer Array (SCUBA) on the James Clerk Maxwell Telescope (JCMT) by Kirk, Ward-Thompson & André (2005). Model analyses of the 1.2 mm continuum have been done in two ways in the literature. One involves a presumed form of the number density radial profile (where cm is the central number density, is the radius, is the reference radius and is the scaling index) plus a constant dust temperature K (e.g. Tafalla et al., 2004, 2006). The other is a static isothermal Bonnor-Ebert sphere (Bonnor, 1956; Ebert, 1955) by Tafalla et al. (2002, 2004). But the invoked Bonnor-Ebert spheres for clouds appear to be dynamically unstable given the inferred physical parameters from observations.

These previous analyses offer relevant information for L1517B, with a constant temperature of  K and a number density of in the inner region. Information has also been drawn from the 850 m continuum observation of L1517B (e.g. Kirk et al. 2005), including a temperature estimate of  K, a full width at half maximum (FWHM) diameter of  pc (by a 2D Gaussian fit to the 850 m data), a central volume number density of  cm and a total cloud mass of within the 150-arcsec aperture. Radio continuum mappings on shorter wavelengths (850 m and 450 m in Kirk, Ward-Thompson & André, 2005) reveal more substructures (thus not that spherically symmetric as compared to the mapping of 1.2 mm continuum emissions) and may be used to constrain theoretical models (see subsection 4.2). Polarization observations for 850 m and 450 m emission continua are performed using SCUBA by Kirk, Ward-Thompson & Crutcher (2006), from which a magnetic field strength of G in this cloud core is inferred. For the moment, magnetic field is not included in our current model formulation.

Various molecular transition lines from L1517B have been observed using the IRAM 30 m telescope by Tafalla et al. (2004, 2006), and spectral line profile fittings based on a static gas sphere close to the center with an empirical density profile, a constant temperature and an ad hoc outer () flow velocity gradient are performed therein. They also argued for other scenarios, such as rotation and pure contraction with certain cooling mechanisms near the center to generate such asymmetric profiles (e.g. Tafalla et al., 2004). There are inconsistencies in such empirical approach from the theoretical perspective as also discussed in Gao, Lou & Wu (2009). Nevertheless, these model fittings might provide gross information for L1517B. In Tafalla et al. (2004), they declare the existence of internal motions of the order of in addition to turbulence, which we shall see in the following analysis might be related to the cloud systematic infall and expansion motions predicted by our general polytropic EECC shock hydrodynamic model. Another important result is the finding of molecular abundance drops towards the cloud center (e.g. Tafalla et al., 2006), which contributes to molecular line profile structures and will also be adopted in spectral line profile fittings of our model analysis (see subsection 3.4 for details).

2.2. Self-similar EECC Shock Model Parameters

Scaling (km s) (yrs)

Parameters  aafootnotemark:  bbfootnotemark:  ccfootnotemark:
Values 1.2 10.02 4.58
11footnotetext: is the polytropic index of the EoS of a general polytropic gas.22footnotetext: Two integration constants and are mass and velocity parameters in the asymptotic solution for [see eqns. (A11) in Appendix A] adopted as the ‘boundary condition’.33footnotetext: is the dimensionless upstream location of an outgoing shock.
Table 1Six independent physical scaling parameters and dimensionless self-similar EECC model parameters

We adopt the general polytropic self-similar EECC hydrodynamic shock model of Wang & Lou (2008) to simultaneously fit both the (sub)millimeter radio continuum emissions and several available molecular line profiles, especially those profiles with red skewed peaks and with spatial resolutions. The general polytropic model can semi-analytically and numerically describe self-similar behaviors of the hydrodynamic evolution, i.e., a time-evolving molecular cloud preserves its basic structural profiles (e.g. the density, temperature, velocity and pressure profiles remain similar to their initial profiles), of a quasi-spherical polytropic gas under self-gravity with the specific entropy conserved along streamlines (see Appendix A for basic nonlinear hydrodynamic partial differential equations (PDEs), the self-similar transformation, notations, and definitions therein).111To explore the origin of “red profiles”, we actually reduce the number of free model parameters in our data fittings by requiring a constant specific entropy everywhere (i.e. is introduced), corresponding to the conventional polytropic cloud core. This constraint is not necessary in the theoretical framework of our general polytropic model (Wang & Lou 2008). A gas thermal temperature of K (e.g. Tafalla et al., 2002; Kirk, Ward-Thompson & André, 2005) and a typical cloud radius of  pc Kirk et al. (2005) are considered for setting pertinent physical scalings in our self-similar EECC shock model. As L1517B appears in an early evolutionary phase of protostar formation, we may estimate a characteristic infall age of about yrs for this cloud core according to Kirk et al. (2005) and Myers (2005). We also adopt an outer radius of AU for our cloud core model of L1517B, by referring to the observed radius (e.g. Tafalla et al., 2002) at an estimated distance of  pc (e.g. Elias, 1978).

Using these empirical information of L1517B, we choose the self-similar model parameters and the physical scalings of time as well as sound parameter as listed in Tables 1 and 2, of which the sound parameter is derived from the following empirical length scale


The number density scale of L1517B cloud core for hydrogen molecule H is taken as [eq (A6) in Appendix A]


where we adopt the mean molecular weight as used by Harvey et al. (2003). This number density scale is one order of magnitude smaller than the empirical values for the central number density (e.g. Tafalla et al., 2002, Kirk et al. 2005), as the density radial profile of our dynamic model increases rapidly towards the core center (see Fig. 1).

Scaling   (km s) a afootnotemark:
 b bfootnotemark:  c cfootnotemark:  d dfootnotemark:
0.8 0 0.039
11footnotetext: The downstream sound parameter .22footnotetext: Parameter is determined by the polytropic index and the scaling index adopted in the self-similar EECC hydrodynamic shock model; this corresponds to a constant specific entropy everywhere.33footnotetext: The reduced central point mass is obtained by solving the self-similar hydrodynamic equations numerically and by matching the central free-fall solution (see Appendix A for details).44footnotetext: Parameters (regarded as independent in Table 1), , , , and are the dimensionless reduced upstream and downstream shock locations, velocities and densities respectively, and are obtained from hydrodynamic shock jump conditions [see eqns. (B2)(B4)] in the self-similar shock model in Appendix B.
Table 2Derived physical scaling and dimensionless self-similar EECC hydrodynamic shock model parameters

With these specified scalings, physical variables including radius, radial flow velocity, number density and thermal temperature of L1517B can be expressed in terms of the dimensionless self-similar variable and reduced dependent variables of using eqns (A5)(A7), viz.


The enclosed mass and the central mass accretion rate [see eqns. (A6) and (A8) in Appendix A] are




where subscripts 1 and 2 refer to physical variables on the immediate upstream and downstream sides of the outgoing shock front, respectively. As required by the mass conservation across a shock surface, we confirm that by a direct numerical computation with model parameters. To simplify physical units of variables, we have already substituted values of and as listed in Tables 1 and 2 into eq. (A8).

Figure 1.— Structures of cloud core L1517B. From top to bottom are radial profiles of radial flow velocity in unit of km s (positive values for radial infall), number density in unit of cm and temperature in unit of Kelvin (K), respectively. The abscissa is radius in unit of AU in a logarithmic scale. Infall radius and shock radius expand at speeds of km s and km s, respectively. Variables exhibit discontinuities across the shock front. Other parameters for this general polytropic self-similar EECC shock solution are summarized in Table 3. The radial density profile of our dynamic polytropic sphere with expanding envelope and free-fall collapsing core thus consists of two parts connected by an expanding stagnation surface, of which the inner part can be described as [eq. (A12)] and the outer part as [eq. (A11)] (see the middle panel here). This broken power law radial density profile is a typical feature of starless cores (e.g. Caselli et al., 2002).
Variables  a afootnotemark: ( yr) b bfootnotemark: (AU) (AU) c cfootnotemark: (km s) (km s) d dfootnotemark:

Note. – These physical parameters are grossly consistent with those estimated from SCUBA observations by (Kirk, Ward-Thompson & André, 2005), except for the conspicuous feature of the variable temperature profile obtained self-consistently from our self-similar EECC hydrodynamic shock model, which differs from the static and isothermal assumptions in earlier models.
(a) and are the central point mass and the total cloud mass within a radius AU.
(b) denotes the central mass accretion rate which decreases with time in our EECC shock model (see equation A8).
(c) is the expanding boundary separating the core collapse and envelope expansion regions and stands for the shock radius. (d) and are the upstream and downstream radial velocities across the outgoing shock front, with negative values for local inflows

Table 3Physical properties of star-forming cloud core L1517B derived from our polytropic hydrodynamic EECC shock model

All dynamic model parameters in Table 1 are chosen in the procedure of data fitting of molecular spectral lines and dust emission continua from observations as described in Sections 3 and 4. In particular, red profiles, e.g. HCO() and HCO transitions, appear suggestive of an expanding envelope in L1517B; this motivates us to invoke a self-similar polytropic EECC shock solution. Meanwhile, in order to fit the molecular spectral line profiles of these emissions, we introduce an outgoing shock in the EECC model to avoid the outer envelope of this molecular cloud expanding too fast. This shock emerges as an expanding flow rushes into the envelope (see the top panel of Fig. 1). We have systematically explored a wide range of model parameters (including polytropic index , asymptotic behavior characterized by mass and velocity parameters and , and location or speed of the outgoing shock) to identify the best-fit EECC shock model. The finally chosen EECC shock model with parameters in Table 1 does self-consistently and simultaneously fit (sub)millimeter emission continua and molecular spectral line profiles of L1517B. Velocity, density and temperature radial profiles at the present epoch of the star-forming L1517B described by the chosen polytropic EECC shock model are displayed in Fig. 1. By the very dynamic nature, these radial profiles in Fig. 1 qualitatively differ from those of both fitting models adopted by Tafalla et al. (2004, 2006) and Kirk et al. (2005). In principle, all these model fits should be constrained simultaneously by other available observations.

As fitted by the underlying hydrodynamic shock model, L1517B appears to involve a self-similar core collapse with an envelope expansion at a typical outflowing speed of km s in the radial range of AU AU, and a concurrent core collapse within the infall radius AU with a typical infall speed of km s. These infall and expansion speeds are in order-of-magnitude agreement with the km s internal speed estimated in Tafalla et al. (2004). There is an outgoing shock at AU with a travelling speed of km s. The central point mass which represents the mass of the pre-protostar is and the total cloud mass within a radius of 21000 AU is . Thus, there is not yet a star (in the sense of thermal nuclear burning) at the center of L1517B since the inferred mass of the pre-protostar is much smaller than the threshold mass value () to initiate nuclear reaction. Besides, according to the standard equation for the accretion luminosity (e.g. equation 2 in Kenyon & Hartmann, 1995), we can estimate an accretion luminosity of about for L1517B, which is likely below the current infrared (IR) detection limit. Moreover, according to Caselli et al. (2002), starless cores are less massive (for example, a statistical mean mass is estimated as , while the total mass of L1517B inferred from our model is ) than cloud cores with nuclear burning stars (e.g., in contrast). Therefore, L1517B can be consistently considered starless, rather than a core with a star, as noted before. Our inferred total cloud mass is twice that derived from dust continuum flux densities by Kirk et al. (2005), but should still be classified as low-mass (, cf. McKee & Ostriker, 2007) star-forming cloud cores. The ratio of pre-protostar mass to the cloud mass is very low, at , which implies that L1517B is in an early phase of protostar formation. The current central mass accretion rate of L1517B obtained from eq. (12) is yr, giving an evolution timescale of yr. From the polytropic model with , the mass accretion rate decreases with increasing time (see eq. A8), which appears to be a typical feature in low-mass star formation systems (e.g. Schmeja & Klessen, 2004; Evans et al., 2009). Thus the realistic infall age of the cloud should be less than the derived evolution timescale . Derived physical properties of the pre-protostellar core L1517B are summarized in Table 3 for a convenient reference.

3. Molecular Spectral Line Profiles

Molecular line profile observations reveal rich details for physical properties of protostellar cores, including important information of the cloud kinematics (e.g. Pavlyuchenkov et al., 2008). We perform spectral line radiative transfer (LRT) calculations for molecular line spectra using the publicly available numerical code RATRAN (Hogerheijde & van der Tak, 2000, see website, which deals with both LRT and non-local thermal equilibrium (non-LTE) excitations of molecular energy levels based on Monte Carlo method. This Monte Carlo code can handle both optically thin and thick lines, barring very large optical depths (e.g. ) caused by convergence problems. We apply the one-dimensional (1D) version of RATRAN code for spherical geometry. We adopt the EECC shock hydrodynamic model with parameters specified in Section 2 for molecular spectral line profile calculations. Here, the radial flow velocity , number density and gas temperature are all derived self-consistently from this EECC model. The dust temperature in the cloud core is assumed to be equal to the gas temperature as expected for high-density cloud cores where gas molecules and dusts have sufficiently frequent collisional exchanges (e.g. Goldsmith & Langer, 1978; Goldsmith, 2001). When running the RATRAN code, we divide the spherical cloud core into 16 uneven shells with enough accuracy for calculations222Actually, 12 and 14 uneven shells have also been tested separately in running the RATRAN code and one can sense the gradual convergence to the final results with little variance. based on the principle that physical quantities in each shell do not vary significantly. We input the number density , gas and dust temperature , radial flow velocity derived self-consistently from our EECC shock hydrodynamic model for radiative transfer calculations. Besides, the temperature of background radiation is the cosmic microwave background (CMB) temperature K. All molecular line transition data are obtained from the Leiden Atomic and Molecular Database (LAMDA, see Schöier et al., 2005), and various molecular abundance ratios with respect to H molecules are presented in Table 4 and more details of abundance profiles can be found in subsection 3.4.

3.1. Red Skewed Double-Peak Line Profiles

Figure 2.— Central molecular spectral line profiles of L1517B with salient red skewed characteristics. The ordinate stands for the brightness temperature in Kelvin and the abscissa represents the LOS velocity component. Histograms are observational data taken from Tafalla et al. (2004, 2006) while solid curves represent RATRAN LRT computational results based on our EECC shock model. Relevant abundances with respect to H are listed in Table 4. Our spectral profile fitting for the molecular transition line HCO shows slightly stronger red peak than observed. However, spatially resolved observation of this line transition reveals that the profiles in some spatial positions [see e.g. () in Fig. 5] exhibit somewhat stronger red peaks than our calculation in the central position due to non-spherical effects in L1517B. For molecular line transition CS(), the intensity of the red peak is more or less the same between the observation data and our fitting, while the fitting result of Tafalla et al. (2004) is lower than the observation data. These red asymmetric spectral line profiles reveal the presence of global expansions for the cloud core envelope.

Numerical results (based on our EECC model) of central molecular spectral line profiles in comparison with relevant observational data from Tafalla et al. (2004, 2006) are displayed in Figs. 2, 3 and 4. Emission lines of HCO and HCO are averaged over a beam area of about to match with Five College Radio Astronomical Observatory (FCRAO) 13.7 m telescope observations in 2001 April (e.g. Tafalla et al., 2006). Spectral profiles for the following eight molecular transition lines, namely HCO, HCO(), HCO(), CS(), CS(), DCO, SO() and SO(), are averaged over a typical beam area of about to match with the data acquired using the IRAM 30 m telescope in Spain between 1999 October and 2002 November (Tafalla et al., 2002, 2004, 2006). Noting that the observational data of the molecular transition line CS() published in Tafalla et al. (2004) actually differ from those of Tafalla et al. (2002), we here take the observational data of Tafalla et al. (2004) in our model analysis.333In our fittings of spectral profiles from emission lines of CS, we note that the observed intensity of the red peak of CS() in Tafalla et al. (2002) appears too low to simultaneously fit the spectral profiles of CS() and CS() transition lines. We then checked the observations of Tafalla et al. (2004) and found that Tafalla et al. (2004) have made a correction of the relevant data, and the new data appears more consistent with our model fitting calculations. The difference between the CS data of Tafalla et al. (2002) and Tafalla et al. (2004) may be caused by different radio telescopes between FCRAO 14m observations (Tafalla et al. 2002) and IRAM 30m observations (Tafalla et al. 2004). Following Tafalla et al. (2004), we produce the NH and NH spectra by averaging over a beam area of and , respectively. As the NH and lines were observed simultaneously using the 100 m telescope of the Max Planck Institute for Radio Astronomy (MPIRA) at Effelsberg near Bonn between 1998 October and 2001 May (Tafalla et al., 2002), we made an average over a beam area to simulate relevant results. A constant intrinsic line broadening with Doppler b-parameter (defined as 1/e the half-width of a line profile in executing the RATRAN code)444Both the thermal () and turbulent () components that contribute to the line broadening have been included in RATRAN code calculations. is estimated by with consistently obtained from each shell and is the mean particle weight, and (Doppler b-parameter). The total line broadening . km s (i.e. with a FWHM of km s for a Gaussian line profile) and a cloud receding velocity of km s are used for all the data fitting procedure.

Molecular Line Frequency Abundance
Transitions (GHz)
HCO(J) 89.188523  (a) (a)footnotemark:
HCO(J) 267.557619  (a) (a)footnotemark:
HCO(J) 140.839502  (a) (a)footnotemark:
HCO(J) 150.498334  (a) (a)footnotemark:
CS(J) 97.980953  (a) (a)footnotemark:
CS(J) 146.969026  (a) (a)footnotemark:
HCO(J) 86.754288  (a) (a)footnotemark:
DCO(J) 216.112582  (a) (a)footnotemark:
SO(NJ) 99.299890  (a) (a)footnotemark:
SO(NJ) 138.178670  (a) (a)footnotemark:
NH(JFF) 93.176258  (b) (b)footnotemark:
NH(JFF) 279.511811  (b) (b)footnotemark:
NH(JKFF) 23.694501  (c) (c)footnotemark:
NH(JKFF) 23.722680  (c) (c)footnotemark:

Note. – Frequencies of molecular transitions are obtained from the LAMDA database (e.g. Schöier et al., 2005). 11footnotetext: Molecules with central depletion hole. for , and for . For simplification and for keeping a minimum number of parameters in our model construction, we have assumed the same radius of abundance hole, AU, for all molecular species in the fitting procedure. 22footnotetext: Constant abundance without depletion hole. 33footnotetext: The central abundance enhancement with Para-, where is the H number density at radius . We adopt which is the same as Tafalla et al. (2004) for comparisons.

Table 4 Rest-frame molecular line transition frequencies and molecular abundances in reference to H molecules

Optical depths play consequential roles in producing asymmetric spectroscopic signatures in molecular line profiles. Optically thick molecular line transitions of HCO and CS display deeper self-absorption dips than optically thin emissions from molecular line transitions of HCO, which have almost no central dips and show only stronger red shoulders (see Fig. 2). Moreover, the HCO line transition manifests less apparent self-absorption dip in comparison with HCO transition. This contrast is likely due to lower energy level populations on J=3 and J=2 in such a cold interstellar medium (ISM) environment of K, leading to a lower optical depth and source function of the J line transition in HCO (Gao & Lou, 2010).

Emissions from the same molecule but for different energy level transitions strongly constrain the underlying polytropic hydrodynamic cloud core model, because all the dynamic and thermal parameters as well as molecular abundance profiles should remain identical in the LRT calculations using the RATRAN code. Fig. 2 shows the spectral profiles of HCO, HCO and CS each for two distinct transitions from central L1517B. These molecular line spectral profiles present explicit red skewed double-peak signatures (i.e. red profiles), consistent with the plausible existence of self-similar EECC shock phase (Lou & Shen 2004; Shen & Lou 2004; Thompson & White, 2004; Gao & Lou, 2010).555A schematic explanation for the origin of red profiles in molecular cloud cores can be found in Lou & Gao (2011). Such a cloud envelope expansion with core collapse may arise from either large-scale cloud radial oscillations (e.g. Aguti et al., 2007) or molecular outflows driven by protostar embedded in the core (e.g. Thompson & White, 2004). It has been claimed that the presence of molecular outflows is always associated with the evidence of non-Gaussian CO line wings (e.g. Thompson et al., 2004). However, spectra for CO(), CO(), CO() and CO() towards L1517B (see figure 8 of Tafalla et al., 2004) do not manifest such evidence of outflows. Besides, no spatially separated “molecular” outflow lobes has yet been detected in this star forming cloud core. We thus propose that the global envelope expansion with central core collapse revealed in L1517B may originate from damped acoustic radial pulsations on large spatial and temporal scales in molecular cloud cores (Lou & Shen 2004; Lada et al., 2003; Keto et al., 2006; Gao & Lou, 2010; Lou & Gao, 2011).666 More explanations for the coexistence of expansion and collapse are elaborated in introduction and discussion sections.

3.2. Other Relevant Molecular Transition Lines

As consistent checks and further constraints, we also compute central spectra of isotopologues, namely HCO and DCO, of the formyl cation HCO and two distinct line transitions of sulfur monoxide SO using the same EECC shock model in comparison with observation data of Tafalla et al. (2006). Fittings of isotopologue emissions with constant abundance ratios to each other777We simply multiply a constant ratio to the molecular abundance of HCO for its isotopologues in our model profile fittings. provide an effective way to investigate the influence of abundance patterns on spectral profiles. As already noted in subsection 3.1, emissions from identical molecules with different energy level transitions serve as effective evidences to verify our EECC shock model. In extensive numerical explorations, we adopt molecular abundances and isotopic ratios which are somewhat different from those used in Tafalla et al. (2006), e.g. the SO abundance with respect to hydrogen molecule H differs from that adopted by Tafalla et al. (2006) due to the variation in the chosen radius of central depletion hole.

There is no conspicuous appearance of red skewed double-peak profiles in molecular line emissions in Fig. 3; instead, single-peak lines are observed in these optically thin lines as expected. As discussed in subsection 3.1, the absence of red asymmetries and self-absorption dips in spectral line profiles is attributed to the decrease of molecular level populations (which in turn leads to the optical depth decrease). This point is highlighted by comparing the spectral line profile of the optically thin HCO in Fig. 3 with that of the optically thick HCO in Fig. 2 [n.b. the optical depth of HCO is that of HCO]. Therefore, the combined examination of optically thick and thin molecular lines offers an important diagnosis to constrain large-scale structures of low-mass star forming cloud cores (Gao & Lou, 2010; Lou & Gao, 2011).

Compared with observations, our model results in Fig. 3, especially those of DCO and SO, tend to be broader in line profile widths. Actually, model lines of these molecules could be narrowed by reducing their depletion hole sizes888The reduction of the depletion hole (still larger than the collapsing core) would bring extra emission contributions from the denser and hotter inner region, so we need to decrease the total abundance of these molecular species to better fit observations. Due to the lower abundance, emissions from the outer expanding envelope which mainly contribute to the edges of the line spectra are weakened, and this leads to narrower spectral line profiles. . For illustrations, we have decreased the depletion hole radius from to for both DCO and SO. Meanwhile, their abundances are reduced to and , respectively, and then both molecular lines become narrower (see dashed curves in Fig. 3). Physical explanations for the reductions of depletion holes are as follows: The difference between the abundance hole radius of the DCO and that of the HCO results from a central increase in the deuterium fractionation caused by the CO depletion, which partly compensates the DCO freeze out at the inner core (e.g. Tafalla et al., 2006). The non-Carbon-bearing molecule SO does not suffer from the depletion of Carbon components, so that its depletion hole may differ from those of Carbon-bearing molecules. Here, SO provides complementary information to emissions from Carbon-bearing molecules for the protostar-forming cloud core.

Figure 3.— Molecular line profiles of central spectra for four transitions HCO(), DCO(), SO() and SO(). Histograms are observational data taken from Tafalla et al. (2006) while solid curves represent LRT fitting calculations based on the same self-similar EECC shock model. Solid curves represent the model results with the universal depletion hole size for all species. Constant isotopic ratios of C/C and H/D (see Table 4 for details) are adopted in our model fitting analysis. Dashed curves displayed in the DCO and SO line profiles correspond to the simulation results for a smaller depletion hole with , and lower molecular abundances as described in the text. The absence of red asymmetry and self-absorption dip in these molecular spectral line profiles is attributed to the decrease of optical depths because of the lower energy level populations for these molecular line transitions.

3.3. Molecules Without Central Depletion Holes

Figure 4.— Molecular line profiles from the central part of L1517B without depletion holes. Multi-peaks present in both molecular transitions are caused by hyperfine (hf) splittings rather than by the cloud core dynamic collapse. Histograms are observational data taken from Tafalla et al. (2002, 2004). As shown in Table 4, a constant NH (top panels) abundance and a central abundance enhancement pattern for NH (bottom panels) have been adopted in order to fit the observational data.

Molecules with central depletion holes offer important diagnostics to probe the dynamic structure of the outer molecular envelope of L1517B. Meanwhile, molecular tracers of the central core, such as NH and NH, offer essential information to examine the collapsing core, as these molecules do not freeze out onto dust grains at typical core densities and they are present in the gas phase throughout the core (e.g. Tafalla et al., 2004). As shown in Table 4, a constant NH abundance and a variable abundance profile with a central enhancement for NH have been adopted in our RATRAN calculations.

Unlike other molecular transitions discussed in this paper, rotational levels of NH are split into multiple hyperfine (hf) components by the two nitrogen atoms (e.g. Caselli et al., 2002; Tafalla et al., 2002). Such hf splittings make the radiative transfer computation much more complicated. Fortunately, available NH molecular data file containing hf structures is provided in LAMDA, which simplifies our calculations considerably. In all, 15 different transitions of NH (with 7 distinct frequencies due to the overlap of some transitions) for and 29 different transitions (with 26 distinct frequencies) for are actually involved in their spectra. We first compute populations of these hf sub-levels and related radiative transfer processes with regard to each hf transition, and then superpose all the individual hf lines to compare with observations (see Fig. 4). No obvious line broadening deviation from the observed lines has been found in our numerical fittings (either in or ). While a significant contracting central region with inward motions of the order of km s has been predicted by our EECC model, this region is very small as compared to the entire collapsing core ( in volume) and is almost negligible compared to the entire cloud ( in volume). Therefore, this region contributes little to the broadening of molecular line profiles. Such a fairly “high” speed of collapsing region is possible in a cloud since gas materials immediately above the surface of the central dense core (or pre-protostar) approach a state of free-fall.

The ammonia molecule consists of the so-called ortho and para species, which coexist almost independently of each other because normal radiative and collisional transitions do not change spin orientations (e.g. Ho & Townes, 1983). The observed NH and inversion lines arise from para-NH (e.g. Tafalla et al., 2002), so we only focus on this particular species in our modelling. Like NH, ammonia NH also has hf splittings due to the interaction between the electrical quadrupole moment of the nitrogen nucleus and the electric field of electrons. Apart from this major contribution, other weaker interactions, including the I J (where I and J are the nitrogen spin and the total angular momentum of ammonia, respectively) and I J (where I is the sum of hydrogen spins) magnetic interactions, as well as H-N and H-H spin-spin interactions, further split the transition components on the order of (e.g. Ho & Townes, 1983). Unfortunately, molecular data with hf splittings for radiative transfer calculations are not available in LAMDA. We follow the approach of Tafalla et al. (2002) by separating the calculation into two steps, namely, the solution of the level excitation and the prediction of the emergent spectrum. Using the RATRAN code, we first adopt the ammonia data file without hf splittings from LAMDA to obtain the combined population of each rotational energy level. We then compute the population of each hf sublevel999Information about hf sub-levels are taken from the NASA-JPL molecular database. Only the electrical quadrupole hf energy level structure has been considered in the data files. Therefore, 5 distinct components of the transition frequency for both and are taken into account in our radiative transfer modelling. by assuming that the sublevels are populated according to their statistical weights. Once populations for the hf sublevels are derived, we take the full hf structures into account and predict the emergent and spectra (as shown in Fig. 4) by integrating the radiative transfer equation along the LOS and summing up each hf component101010We compute the Einstein A coefficients for different hf line transitions based on equation (9) of Pickett et al. (1998).. The validity of this approach was discussed by Tafalla et al. (2002). We briefly comment here: according to quantum statistics, the population ratio between two different levels is given by , where are the statistical weights (degeneracies) for level (=1, 2) and is the energy difference between the two energy levels. As long as , i.e. the energy difference is far less than the energy scale involving the excitation temperature (which is usually satisfied when we consider the population of hf sub-levels for molecules in starless cores), we can then approximately take the population ratio as simply . The excitation temperature and the kinetic temperature can be equal under thermodynamic equilibrium, but this is not necessarily so in general.

Figure 5.— We show here the grid-map for the observed molecular transition HCO spectra with a step for L1517B. The emission intensities are expressed in main beam temperature (using Kelvin as the unit) with elevation corrections where is the elevation angle of the Delingha telescope, and the efficiency of the Delingha 13.7 m telescope is . This grip map is centered at , based on the central position of the 1.2 mm continuum source of L1517B as observed by Tafalla et al. (2002). The LOS velocity resolution is km s, and the root-mean-square error in the antenna temperature is K. The legend panel to the upper right corner illustrates the central panel as an example. The spectral profiles are not perfectly spherically symmetric with the eastern spectra slightly shifting towards the blue. The difference between the eastern and western spectra might be due to the rotation of L1517B about its north-south axis. However, double-peaked red profiles are evident in HCO lines across the cloud globule and an average over the same distance from the center reveals stronger red peaks in all radial positions (see Fig. 8). In this sense, the presence of red profiles is indicative of cloud core expansive motions (see also Aguti et al., 2007).

Ammonia NH spectral data are usually used to probe the central gas kinetic temperature of star-forming molecular cores (e.g. Tafalla et al., 2004). The database of dense molecular cores mapped in the and transition lines of NH was compiled and presented by Jijina, Myers & Adams (1999), who conclude that the temperature distribution of molecular cores in Taurus complex is very narrow around K. In our polytropic EECC shock model, however, the gas temperature gradually rises up towards the center and reaches K in the collapsing core of L1517B as shown in Fig. 1. According to Tafalla et al. (2004) and Ho & Townes (1983), the gas temperature of K inferred from NH3 observations relies on the presumed (static) density profile, which is an empirical asymptotic power-law envelope with a “flat” central region. It appears that such a constant K temperature reasonably fits the gas kinetic temperature for L1517B derived from NH data analysis (e.g. figure 4 of Tafalla et al., 2004). We note, however, that the gas kinetic temperature thus determined is inferred directly from the rotational temperature of NH (e.g. Walmsley & Ungerechts, 1983; Tafalla et al., 2004) which relates to the relative brightness temperature/population between the and radiative transitions (e.g. Ho & Townes, 1983); therefore, the gas kinetic temperature thus derived only represents an “averaged” temperature along the LOS, not necessarily indicating the actual temperature at a certain spatial point. In other words, it is possible that a radially variable temperature profile could also reproduce the NH data as shown by our dynamic model fitting analysis (see also Galli et al. 2002 for variable temperature profiles in molecular cloud cores). Noting that the temperature profile of our model reaches K only in the very inner region ( 1000 AU from the center111111The gas kinetic temperature for L1517B derived from NH observations starts from 900 AU (see figure 4 of Tafalla et al., 2004); therefore, the “hot” region in our model occupies the inner core and contributes only a very small part to the LOS averaged gas kinetic temperature.) while remaining around K over a wide radial range in the outer portion (see Fig. 1), we would expect that our radially variable temperature profile may mimic the constant LOS temperature distribution of the NH analysis for L1517B. Besides, as the temperature distribution of L1517B in figure 4 of Tafalla et al. (2004) represents LOS averages spatially smoothed with 40”, these averages over the telescope beam resolution would further flatten the temperature distribution and make the whole profile closer to a roughly constant K distribution. Moreover, as our EECC model shown in Fig. 1 as a whole reasonably reproduces both numerous molecular line emissions (including those of NH and NH) and (sub)millimeter dust continuum observations of L1517B (see Table 4 and Section 4 for details), we conclude that this hydrodynamic self-similar EECC model is a viable physical description of cloud core L1517B. Our temperature variation would only introduce minor corrections, because the inner “hot” region is fairly small compared with the entire molecular cloud (). In other words, a variable temperature profile might appear somewhat flattened in this type of inferences by LOS NH observations. Noting that our radial temperature profile drops below K in the outer envelope portion, it is possible to produce a 10 K “averaged” rotational temperature from our EECC model. Jijina et al. (1999) have adopted the - (here is the gas kinetic temperature) from Walmsley & Ungerechts (1983) to convert the inferred rotational temperature into gas temperature, which may cause further uncertainty in the determination of cloud temperature since the - relation of Walmsley & Ungerechts (1983) may not be that accurate (e.g. Tafalla et al., 2004).

3.4. Influence of Molecular Abundance Distributions

The spatial distribution of molecular abundances plays a crucial role in producing molecular spectral line profiles and thus has received considerable attention in the literature (e.g. Herbst & Klemperer, 1973; Rawlings & Yates, 2001; Tsamis et al., 2008). Also, the depletion of molecular species by adhesion onto cold dust grain surfaces has been discovered in central regions of star-forming cloud cores (e.g. Tafalla et al., 2002, 2004; Walmsley et al., 2004). We have selected molecular abundance ratios (relative to the number of H molecules) by referring to Tafalla et al. (2004, 2006) yet with variations to probe molecular distributions in L1517B. Central molecular depletion holes with very much lower abundance ratios [i.e. of the ratio in the outer cloud layers as in Tafalla et al. (e.g. 2006)] are essential in reproducing absolute intensities and relative strengths of ‘blue’ and ‘red’ peaks in molecular spectral line profiles. It is possible to treat the radii of abundance holes, which can vary for different molecular species in general, as extra independent parameters as in Tafalla et al. (2006). However, in order to focus on exploring the dynamic process of L1517B and to reduce the number of independent parameters, we have assumed a common size of depletion hole for all molecules, except for NH and NH molecular transitions. Our chosen abundance distributions and the radius of the depletion hole used in the data fitting are summarized in Table 4 for reference. Based on LRT RATRAN calculations, intensities of both ‘blue’ and ‘red’ peaks would decrease with the increase of ; this is also intuitively sensible. In fact, the intensity of the stronger ‘red’ peak appears to be much more sensitive to the variation of than that of the weaker ‘blue’ peak does. This may allow us to estimate the radius of the depletion hole in L1517B by carefully calibrating the ‘blue’ to ‘red’ peak intensity ratios in pertinent molecular spectral line profiles.

3.5. Spatially Resolved Spectral Line Profiles

Figure 6.— A direct comparison between the central spectra, which have been generated from the average of typically 5 spectra within a radius from the core center, of the molecular transition HCO from L1517B observed by Delingha 13.7 m millimeter-wave radio telescope in mid-April 2010 (solid histogram) and by FCRAO 13.7 m telescope in April 2001 (dashed histogram). The ordinate stands for the main beam brightness temperature converted using the beam efficiencies (0.62 for Delingha telescope and 0.55 for FCRAO telescope respectively) from the respective antenna temperatures . The Delingha observation has a velocity resolution of 0.04 km s as compared with the velocity resolution between 0.03 km s and 0.07 km s of the FCRAO observation, which used the QUARRY array receiver in frequency switching mode together with the facility correlator. With longer integration times for the observation of each spatially resolved regions, the result from the Delingha observation has relatively higher signal to noise ratio (S/N) compared with the FCRAO data.
Figure 7.— The HCO grid-map of four spectral profiles at a distance of around the central spectral profile plot for L1517B (see Fig. 8 for the average of these four spectra). Emission intensities are expressed in the main beam brightness temperature . The net integration time for each spectral profile with the impact parameter b is 50 minutes, and is 80 minutes for the spectral profile at the center. While these spacing spectra overlap previous observations with the beam size of , it is still valuable to examine global dynamic properties of L1517B from these spectra, as the effect of possible rotation about the north-south axis in the inner cloud region appears weaker than in the outer region. Compared with Fig. 5, apparent red profiles in nearly all these inner spectra (n.b. the red asymmetry is not particularly obvious in the eastern spectrum) strongly indicate that this red asymmetry is an intrinsic characteristic of L1517B caused by its internal motion rather than by rotation.

Central spectral profile model fittings of various molecular line transitions provide an important and effective approach to probe both dynamical and thermal structures of pre-protostellar cores. However, only presenting signatures from the cloud center, central spectral profile fittings of molecular lines might still have certain ambiguity and uncertainty in determining all pertinent physical properties of star-forming cloud cores. Mapping molecular transition line profiles with high enough spatial resolutions, on the other hand, would further constrain the large-scale physical conditions within molecular clouds (e.g. Aguti et al., 2007; Gao & Lou, 2010; Lou & Gao, 2011). Here, we report a position-switch observation ( grids) of spectral profiles for the HCO line transition, which is the most prominent diagnostics for dynamic motions within L1517B, using the 13.7 m telescope at Delingha in Qinghai Province of China and present our model fitting results of spatially resolved spectral profiles based on the same underlying EECC shock model. Meanwhile, we predict other pertinent red skewed spectral profiles of molecular lines with spatial resolutions for future observational tests.

3.5.1 Data Acquisition and Comparison

We recently observed L1517B for HCO at 89.186769 GHz using the 13.7 m millimeter-wave radio telescope of PMO at Delingha during April 1222, 2010. The telescope is m above sea level with an extremely dry air. One SIS receiver operating at GHz was used in the observation together with a FFTS digital spectrometer with 16384 channels and a working bandwidth of 200 MHz, which gives a velocity resolution of 0.04 km s at the frequency of 89 GHz. The position-switch mode was chosen, with the pointing and tracking accuracy within the range of . The telescope beam size is approximately 50 at this frequency, and we adopted a mapping step size of 50 in observation. The typical system temperature is K during our observations. The standard chopper wheel calibration was used during the observation runs to get the antenna temperature , which has been corrected for the atmospheric absorption and telescope elevations. The beam efficiency of 0.62 is used to convert the telescope antenna temperatures into main beam brightness temperatures . The net integration time for most positions is 60 minutes, and is 80 minutes for the central position, resulting a root-mean-square (rms) noise of K for the measured antenna temperature .

The spectral profile data are processed with the standard CLASS software. Linear baselines are removed from all spectra and an elevation correction ( is the elevation angle) is adjusted.

We also made an average of five spectra within a 50 radius from the cloud core center (i.e. the central spectrum together with four spectra that are away from the center) to directly compare with the result of Tafalla et al. (2006) obtained by the FCRAO 13.7 m telescope in April 2001. From the data comparison of Fig. 6, we see that these two observational results are generally consistent with each other with the notable exception of integrated intensities around the dip between the two peaks. This difference may be tolerated within uncertainties of the observations and might be due to slight discrepancies in telescope calibrations or velocity resolutions. This comparison of data from two independent observations 9 years apart confirm the validity of both observational data from FCRAO 2001 April and Delingha 2010 April.

To examine whether these “red-profiles” are caused by rotation, we have observed four spectral profiles at a distance of 20 from the center (Fig. 7) since the rotational effect is much weaker in the inner region than in the outer region. Conspicuous red asymmetric profile characteristics in nearly all these spectra exclude the possibility that a pure rotation produces these “red-profiles”.

3.5.2 Model Fittings of Observational Data

Simultaneous comparisons between model results and observational data for molecular line profiles with distinct impact parameters (distance of LOS from the core center) is another powerful way to examine the overall dynamical and thermal structures of star-forming cloud cores. As our model is spherically symmetric while the spatially resolved observation (Fig. 5) exhibits non-spherical effects, we make an average over grids of data with the same radial distance from the center. We carried out observations of four spectra spaced by from the core center under the same observational conditions (yet with 50 minute integration time each) to check EECC shock model (Fig. 7). All the dynamic parameters as well as the abundance distribution and the intrinsic broadening in these simulations are consistent with those used in fitting the central spectrum of HCO observed by the FCRAO 13.7m telescope (see Section 2).

From the model fitting results shown in Fig. 8, we see that our computed molecular line profiles based on the general polytropic EECC shock hydrodynamic model coincide with the observational data. In comparison with observations, molecular line profiles produced from our numerical model calculations show somewhat stronger red peaks and weaker blue peaks at the very center (Fig. 8). This likely originates from the spatial asymmetry of the spectral line profiles, which might be caused by the possible rotation of L1517B about its north-south axis. The small blue shifts in the eastern spectra (see Fig. 5) will reduce the intensity of the red peaks and enhance the blue peaks when we make an average of the observational data over the same radial distance. Meanwhile, these averaged spectra exhibit increasing broadening towards the outer envelope while the broadening of our numerical results remain invariant. This phenomenon might be possibly attributed to the incremental influence of turbulence with increasing radius, which is typical in star-forming clouds (e.g. McKee & Tan, 2003). In order to reduce the number of free parameters and to better understand the dynamics of L1517B, we simply assumed that the intrinsic broadenings (which are characterized by the Doppler bparameter in the RATRAN code) take the same value everywhere in the numerical model fitting calculations.

All these averaged spectral profiles present evident red profile characteristics, which are indicative of cloud core expansions. Our fitting results reveal that L1517B is most likely undergoing a core collapse and expansion process as described by the polytropic EECC shock model.

Figure 8.— Model fittings of spatially resolved spectral line HCO observations for L1517B. Histograms are Delingha observational data produced by the average of line spectra with the same distance from the core center. A telescope efficiency of was used to convert the antenna temperature into the main beam brightness temperature . Solid curves present LRT fitting simulations averaged over a beam area of based on the same underlying polytropic EECC shock hydrodynamic model and molecular abundance pattern. Spectra correspond to the impact parameter b (distance of LOS from the core center) being , , , , , and , respectively. The same intrinsic broadening with Doppler bparameter km s and receding velocity of the cloud core km s as being used previously are adopted in our data fitting procedure.

3.5.3 Several Model Predictions for Cloud Core L1517B

Figure 9.— Red-skewed molecular spectral profiles with spatial resolutions for L1517B. Emissions from HCO, HCO(), HCO(), CS() and CS() are averaged over a radius of to simulate the observation of IRAM 30 m telescope. We present emission spectral profiles with the chosen impact parameter being , , , and respectively to reveal the variation of the relative difference in strengths between red and blue peaks as well as the overall magnitudes.

As predictions for future observations of L1517B, we compute spatially resolved relevant molecular line profiles shown in Fig. 9 for molecular line transitions with red asymmetries in their central spectra.

Our model simulation results in Fig. 9 show a gradual disappearance of red asymmetries in spectral line profiles and an explicit decrease in overall magnitudes as the LOS departs away from the cloud core center. This variation trend is attributed to the diminution in both cloud density and temperature, as well as optical depth (Gao & Lou, 2010). Therefore, it would help to assess the influence of optical depths on molecular line emissions by making comparisons between model results and observations for molecular transitions with spatial resolution.

4. Millimeter and Sub-millimeter Continua of Star-Forming Cloud Core L1517B

While model fittings to observed molecular spectral line profiles provide an important means to infer physical properties of cloud cores, molecular line emissions only present under certain conditions and some molecules may be frozen onto dust grain surfaces in certain regions (e.g. Walmsley et al., 2004). Meanwhile, the radial profiles of flux intensities at (sub)millimeter wavelengths from dust grain emissions are generally independent of chemical processes, thus can serve as an independent observational diagnostics to examine and constrain density and temperature distributions in star-forming clouds (see e.g. Bacmann et al., 2000; Shirley et al., 2000; André et al., 2004). Moreover, as the dust density is almost independent of depletion holes of molecular abundances, millimeter and (sub)millimeter continua would offer an additional probe to distinguish different models. For optically thin dust emissions, the integral form of the specific intensity from a spherically symmetric cloud globule along a LOS with an impact parameter is given by


(e.g. Adams, 1991) where is the outer radius of a cloud globule, and , , and are the dust temperature, cloud mass density, specific dust opacity and the Planck function, respectively.

4.1. Radio Continuum Emissions at 1.2 mm

The 1.2 mm wavelength radio continuum mapping of L1517B is centrally concentrated and appears grossly spherical (see figure 1 of Tafalla et al., 2004), so the spherical symmetry is a reasonable first-order approximation for L1517B. We use the same EECC shock hydrodynamic model described in Section 2 to fit the radial profile of 1.2 mm radio continuum emission data (see figure 2 in Tafalla et al., 2004). Tafalla et al. (2004) acquired the 1.2 mm radio continuum mapping from the IRAM 30 m telescope with a beam size of . Therefore, we integrate the intensities from model calculations within this beam size. Parameters of the EECC shock model are shown in Tables 1 and 2, and our theoretical model fitting to the observed radial profile of radio continuum is displayed in Fig. 10. In the model computations, the dust temperature is assumed the same as the gas temperature adopted for previous molecular spectral line profile calculations shown in the Section 4. We utilize the dust opacity model proposed by Ossenkopf & Henning (1994) with emissions and absorptions included but with scatters ignored. A typical gas to dust mass ratio of 100 to 1 is assumed in our cloud model calculations as usual. Numerical calculations are carried out using the RATRAN code with the spherical cloud divided into 256 spherical uneven shells. We have also doubled the number of shells (i.e. 512) in our calculations (see the dashed curve in the top panel of Fig. 10) and find that, except for a tiny increment at small radii, the result generally coincides with that produced by 256 shells. Therefore, our adoption of 256 shells provides sufficient computational accuracy for modeling the continuum emissions. Though the density and temperature of the cloud increase rapidly near the center, they do not lead to a noticeable central peak or a prompt increase in the dust continuum since our calculation results have been convolved with a telescope beam size of which smooths the singularities of density and temperature profiles at very small radii.

4.2. Continuum Emissions for 850 m and 450 m

High-quality sub-millimeter continuum observations at wavelengths of 850 m and 450 m reveal more structural complexities of pre-protostellar cores (e.g. Holland et al., 1999; Shirley et al., 2000; Dye et al., 2008). As a consistent check, we compute radial profiles for 850 m and 450 m flux intensities using the same EECC model and compare them with the observational data of SCUBA by Kirk et al. (2005). The same opacity model for (Ossenkopf & Henning, 1994) and the gas-to-dust mass ratio of 100 to 1 are adopted in these model calculations.

In order to be compatible with the JCMT beam FWHM of arcsec at wavelengths 850 m and arcsec at 450 m, our numerical calculations are convolved with a telescope beam size of for 850 m and of for 450 m, respectively. The fitting result of 850 m normalized flux density is shown in Fig. 10 and the model calculation of the 450 m normalized flux density radial profile is shown in Fig. 11. We find a reasonable agreement of the model results with observations of 850 m continuum inside AU; outside this radius, the quality of data becomes poor with considerable error bars. The 450 m flux density profile is grossly consistent in intensity with the observational data (see figure 3 of Kirk et al. 2005) at corresponding radius, with several substructures being ignored. As we have not included the influence of the central point source in our calculations, peak flux intensities of both 850 m ( 155 mJy beam) and 450 m ( mJy beam) are somewhat lower than those observed, namely mJy beam and mJy beam (see table 1 of Kirk et al. 2005). Since the radio observations of 1.2 mm continuum and 850 m data were carried out by two different telescopes (i.e. the IRAM for 1.2 mm and the SCUBA for 850 m) and our model computational result agrees well with the 1.2 mm continuum of dust emissions (see Fig. 10), the systematic difference of our fitting for the 850 m continuum profile may also be related to differences in calibrations between the two telescopes, which could be checked by future direct radio observations of 850 m continuum using the IRAM 30m telescope.

4.3. Column Number Densities and Dust Extinctions

The near-infrared dust extinction provides an independent diagnostics to trace the column density in molecular clouds without involving the dust temperature distributions (Alves, Lada & Lada, 2001). It serves as an important constraint to the density radial distribution of our general polytropic EECC shock hydrodynamic model, as a complementary yet independent diagnostics to various molecular line profiles and dust radio continuum emission measurements (Lou & Gao, 2011)121212Because the temperature profile is not isothermal (as often assumed in the literature) in our EECC model, the column density radial distribution cannot be directly mapped out by dust (sub)millimeter continuum emissions. . Moreover, column densities in the low-density outer region of molecular cloud, where (sub)millimeter dust continuum emissions may not be easy to detect, can still be traced by near-infrared dust extinction measurements (Kandori et al., 2005). As an essential and necessary check, we thus compute the column density radial profile of H molecules by simply integrating volume densities from the same EECC shock hydrodynamic model. Since the column density radial profile is sensitive to the outer radius of a cloud globule, our model fitting to the observed radial profile would identify a proper outer radius of the cloud core at  AU. We then invoke the empirical conversion relation cmmag (e.g. Kandori et al., 2005; Kirk, Ward-Thompson & André, 2005), in which is the H molecule column number density and is the visual extinction, to predict the column density and visual extinction for cloud core L1517B using our EECC shock model which can reasonably fit other currently available observational data.

Figure 10.— Radial profiles of the 1.2 mm (top) and 850 m (bottom) radio continuum emissions from L1517B. For the 1.2 mm continuum, the intensity scale is in mJy per beam size and the abscissa is the angular distance from the cloud core center in arcsec. Solid squares are data from figure 2 of Tafalla et al. (2004). The solid and dashed curves present the simulation fitting results based on the density and temperature profiles of our EECC shock hydrodynamic model by dividing the cloud into 256 and 512 shells, respectively. For 850 m, we show the normalized radial profile for the flux intensity (per beam size) with the peak flux density 155 mJy beam. Data points from SCUBA observations (Kirk et al. 2005) are shown here in a logarithmic scale at half beam spacings with 1 error bars ( mJy beam). The abscissa stands for radius from the center in unit of AU on a logarithmic scale. The scatter of data points, in both the 1.2 mm and 850 m continua, at the most outer part of the core (, where the inconsistency between our simulations and observations emerge) is highly likely due to the asymmetry of the core (see figure 1 of Tafalla et al., 2004). No conspicuous shock discontinuities are seen for either 1.2 mm or 850 continua, since this shock locates at the outer envelope () where both the temperature and density change slightly due to the shock wave.

In the column density calculation, we convolve numerical integration results with a typical beam resolution of . As a comparison and prediction, results of both the H molecule column density and the visual extinction derived from the EECC shock model are shown in the lower panel of Fig. 11. For example, the value of central column density predicted by our EECC shock model ( cm) is close to yet a bit lower than the model result of Kirk et al. (2005) [i.e. cm with a typical error range of percent (see table 4 of Kirk et al. 2005). Their central “column density”, however, instead of being directly obtained from infrared dust extinction measurements, is actually derived from the 850 m radio continuum profile modelled by a static isothermal ( K ) Bonnor-Ebert sphere and a constant dust opacity cm g. However, with those estimated parameters, such an isothermal Bonnor-Ebert sphere for L1517B turns out to be dynamically unstable. Besides, their central “column density” is model dependent without being constrained by other available observations of L1517B, so we propose here to test and examine their prediction as well as ours for the column density radial profile of L1517B by dust extinction observations.

Figure 11.— Radial profile for the 450 m radio continuum emissions (top); and column densities for hydrogen molecules and corresponding magnitudes of visual extinction (bottom) predicted for observations of L1517B using the EECC shock model. For the 450 m continuum, here we show the normalized radial profile for the flux density (per beam size) with the peak flux intensity mJy beam. Results are shown in a log-log scale with the abscissa being the radial distance in AU. For the column density and visual extinction profiles, the ordinate on the left axis marks the column number density distribution in unit of cm and the ordinate on the right axis indicates the magnitudes of visual extinction . The abscissa is the radial angular distance from the cloud core center in arcsec. These results are convolved with a typical observational resolution of .

5. Conclusions and Discussion

In this paper, we have systematically investigated the physical and chemical (abundance) properties of the starless cloud core L1517B from three complementary observational aspects, viz., molecular spectral line profiles (some of them are spatially resolved), (sub)millimeter continuum emissions and the column density radial profile through near-infrared measurements of dust extinction. We invoke a self-similar polytropic EECC shock model with six independent model parameters involved, viz. the upstream sound parameter , time scale , polytropic index , coefficients of asymptotic solution conditions and , and the shock location , as summarized in Table 1 to present reasonable data fittings to various asymmetric ‘red profiles’ together with optically thin single peak spectral profiles, as well as to 1.2 mm and 850 m continuum radial profiles from observations (Kirk et al. 2005; Tafalla et al., 2004, 2006). In addition, we have recently completed a (with a 50 spacing) spatially resolved spectral profile observation of the HCO line transition using the PMO 13.7 m millimeter-wave radio telescope at Delingha and made a comprehensive comparison with our model results. Furthermore, other predicted molecular spectral line profiles with higher spatial resolutions are expected to be tested by future high spatial resolution spectral observations of L1517B.

The underlying general polytropic self-similar EECC shock hydrodynamic model (Wang & Lou, 2008) provides radial velocity, mass density and temperature structures of the cloud globule consistently. Based on the inferred EECC model with a set of fitting parameters, we provide physical parameters for L1517B. (1) The L1517B molecular cloud core is likely evolving in an early phase of pre-protostar formation with a small central point mass of (incapable of thermal nuclear burnings) and a current central mass accretion rate of yr. (2) The total mass of the cloud core is estimated as revealing that L1517B belongs to a low-mass cloud core (e.g. Tafalla et al., 2002; McKee & Ostriker, 2007). (3) The cloud core L1517B has an expanding envelope with a typical outgoing speed of km s and a core collapsing at km s.131313We show the velocity profile from our model in Fig. 1; these typical values only imply the characteristic magnitudes of expanding and collapsing flow velocities, rather than implying the outgoing/infalling region has an uniform velocity of this typical value. (4) There is also an expanding shock at radius AU with an outgoing radial speed of km s. This outgoing shock is initiated when an expanding outflow runs into a slowly infalling gas in our dynamic EECC model scenario. More importantly, we demonstrate that the polytropic EECC shock dynamic phase of L1517B with a set of sensible parameters can give rise to the manifestation of various observed molecular spectral line profiles.

In contrast to the results of Tafalla et al. (2004, 2006) and Kirk et al. (2005) (see section 2.1 for details), our general polytropic EECC shock hydrodynamic model simply depends on physical properties of L1517B inside the observed outer boundary of and self-consistently produces thermal as well as dynamic properties of the cloud core in the framework of self-similar hydrodynamics to fit the molecular spectral line profiles, (sub)millimeter continuum radial profiles and to predict the dust extinction property and the column density radial profile in the plane of sky.

The uniqueness of our model fitting is generally supported by the theoretical work done by Gao & Lou (2010) and Lou & Gao (2011), which discussed the relationship between different sets of model parameters and the corresponding molecular line spectral profiles. Instead of other possibilities that might cause red-skewed molecular spectral line profiles as noted in the introduction of Gao & Lou (2010), we conclude that the ‘red’ asymmetries present in several molecular central spectra of L1517B are most plausibly produced by the coexistence of core collapse and envelope expansion motions in the cloud. Other scenarios, such as rotation and pure contraction without expansion were suggested by Tafalla et al. (2004). However, as discussed in Section 3.5.1, spatially-resolved molecular spectra around the core center reveal that this asymmetry is most likely an intrinsic property of the cloud dynamics, rather than being caused by a rotation. Besides, as the temperature profile in our model increases towards the center (a typical characteristics for most molecular clouds), it is impossible to produce molecular line profiles of red asymmetry from a pure contraction (see eq. B4 of Gao & Lou, 2010).

As indicated by Broderick et al. (2002), there exist radial breathing modes generated from nonlinear evolution of acoustic radial pulsations on large spatial and temporal scales, with long crossing times in protostar-forming molecular cloud cores. We advanced the following scenario for the possible emergence of the EECC phase of L1517B. The molecular cloud is undergoing an envelope expanding phase of the breathing mode when the central contracting region collapses due to nonlinear instabilities. This central infalling dynamics is promoted by self-gravity and the collapsed region expands towards the envelope, similar to the EWCS in Shu (1977). Eventually, all gas materials would continue to fall towards the center to form a proto-stellar core. This scenario is consistent with both theoretical (e.g. Stahler & Yen, 2010) and observational (e.g. Lee & Myers, 2011) pictures concerning the evolution of internal motions in starless cores.

We find in our model analysis that a comparison between optically thin and thick molecular spectral line profiles (or between different energy level transitions of the same molecule) serves as an effective and sensible method to study the influence of optical depth variations. In addition, analyses of the same molecular line transition along LOS with distinct impact parameter from the center may further provide an effective approach to examine effects of both optical depths and radial variations of cloud physical properties (see Section 3.5). We have prescribed the abundance patterns of different molecules with central depletion holes, which has been widely invoked for molecular line profiles from star-forming clouds as radiative diagnosis (e.g. Tafalla et al., 2002, 2004; Walmsley et al., 2004). Intensities of ‘red’ peaks appear to be much more sensitive to the size of a depletion hole, as compared to ‘blue’ peaks. This may allow us to assess molecular depletions by comparing the relative differences between intensities of ‘red’ and ‘blue’ peaks in their molecular spectral line profiles.

We thank the anonymous referee for suggestions to improve the quality of the manuscript. We would like to thank the hospitality of PMO Delingha Observatory and the 13.7 m telescope staff for their support during the observation. D.R. Lu and Y. Sun are acknowledged for assistance in the data processing. We thank J. Yang of PMO for advice and suggestions. This research was supported in part by the National Undergraduate Innovation Research Program 091000344 for two consecutive years at Tsinghua University from the Ministry of Education.

Appendix A    General Polytropic Self-Similar Hydrodynamic Model

In spherical polar coordinates , general polytropic nonlinear hydrodynamic partial differential equations (PDEs) for a spherically symmetric molecular cloud under the self-gravity and gas pressure force are given by


where the mass density , the radial bulk flow velocity , the thermal gas pressure and the enclosed mass within radius at time depend on and in general; dyne cm g is the gravitational constant and is the polytropic index. PDEs (A1) and (A2) are the two complementary forms of the mass conservation. PDE (A3) is the radial momentum conservation of a molecular cloud under the gas pressure force and the self-gravity but in the absence of random magnetic fields (Wang & Lou, 2007, 2008). PDE (A4) requires the specific entropy conservation along streamlines corresponding to a general polytropic EoS with formally related to the specific entropy that varies in both time and radius in general.

To derive an important subset of nonlinear self-similar solutions, these hydrodynamic PDEs can be cast into a set of coupled nonlinear ordinary differential equations (ODEs) with the following transformation (Wang & Lou, 2008),


where is the independent self-similar variable combining and in a proper manner, and , , and are the dimensionless reduced radial flow speed, mass density, thermal pressure and enclosed mass, respectively. Two constants and are the sound parameter and scaling index which consistently make , , , and dimensionless and control the spatial and temporal scalings of physical variables in a hydrodynamic cloud. The thermal gas temperature is given by the ideal gas law


where , and are Boltzmann’s constant, the hydrogen atomic mass and the mean molecular weight, respectively. For a finite as in the central free-fall solution (e.g. Wang & Lou, 2008), the last expression in equation (A6) gives the central mass accretion rate as


with being the reduced central enclosed point mass. For smaller or larger than 1, this central mass accretion rate decreases or increases with increasing , respectively. For , the central mass accretion rate is constant.

Substituting self-similar transformation (A5)(A6) into PDEs (A1)(A4), we obtain two coupled nonlinear ODEs


with for , where parameter . These ODEs are related to the mass conservation, the radial momentum conservation and the specific entropy conservation along streamlines, respectively. In the limit of , we have the asymptotic self-similar solution to the leading orders as


(see Lou & Shi 2011 in preparation for more specific comments on the term) where and are two integration constants, referred to as the mass and velocity parameters, respectively.

In the other limit of and to the leading order, we have the asymptotic central free-fall solution


where the constant is the reduced enclosed point mass, representing the dimensionless proto-stellar mass.

By solving coupled nonlinear ODEs (A9) and (A10) with analytic asymptotic solutions (A11) and (A12) as “boundary conditions”, and by taking proper care of the sonic critical curve, we can derive the radial profiles of velocity, density and thermal temperature simultaneously and self-consistently [see Wang & Lou (2008) for more details].

Appendix B Hydrodynamic shock conditions in the self-similar form

Self-similar shocks may appear in dynamic molecular clouds and can be constructed within the framework of a general polytropic hydrodynamic model. For a hydrodynamic shock, we apply the three conservation laws in the shock comoving framework, namely conservations of mass, radial momentum and energy, across the shock front


where , , and represent the flow velocity, shock front speed, gas mass density and thermal gas pressure, respectively. We use a pair of square brackets outside each expression embraced to denote the difference between the upstream (marked by subscript ‘1’) and downstream (marked by subscript ‘2’) quantities, as has been done conventionally for shock analyses (Landau & Lifshitz, 1959).

As the parameter in the self-similarity transformation equation (A5) is related to the sound speed which are generally different in the upstream and downstream sides of a shock, there are two parallel sets of self-similar transformation in the two flow regions across a shock. We set with the dimensionless ratio representing this difference. The relation is then required for consistency. With this similarity scaling relation, hydrodynamic shock jump conditions (B1) can be readily cast into the following self-similar form (Wang & Lou, 2008), namely


These three self-similar shock conditions (B2)(B4) can be solved explicitly and more relevant details can be found in Appendix D of Wang & Lou (2008) in the absence of random magnetic fields.


  • Adams (1991) Adams, F. C. 1991, ApJ, 382, 544
  • Aguti et al. (2007) Aguti, E. D., Lada, C. J., Bergin, E. A., Alves, J. F., & Birkinshaw, M. 2007, ApJ, 665, 457
  • Alves, Lada & Lada (2001) Alves, J., Lada, C. J., & Lada, E. A. 2001, Nature, 409, 159
  • André et al. (2004) André, P., Bouwman, J., Belloche, A., & Hennebelle, P. 2004, ApSS, 292, 325
  • Bacmann et al. (2000) Bacmann, A., André, P., Puget, J. L., Abergel, A., Bontemps, S., & Ward-Thompson, D. 2000, A&A, 361, 555
  • Benson & Myers (1989) Benson, P. J., & Myers, P. C. 1989, ApJS, 71, 89
  • Bodenheimer & Sweigart (1968) Bodenheimer, P., & Sweigart, A. 1968, ApJ, 152, 515
  • Bonnor (1956) Bonnor, W. 1956, MNRAS, 116, 351
  • Broderick et al. (2002) Broderick, A. E., Keto, E., Lada, C. J., & Narayan, R. 2007, ApJ, 671, 1832
  • Caselli et al. (2002) Caselli, P., Benson, P. J., Myers, P. C., & Tafalla, M. 2002, ApJ, 572, 238
  • Dye et al. (2008) Dye, S., et al. 2008, MNRAS, 386, 1107
  • Dyson & Williams (1997) Dyson, J. E., & Williams, D. A. 1997, The Physics of the Interstellar Medium. 2nd ed., IOP Publishing Ltd., Bristol and Philadelphia
  • Ebert (1955) Ebert, R. 1955, Z. f Astr., 37, 217
  • Elias (1978) Elias, J. H. 1978, ApJ, 224, 857
  • Evans et al. (2009) Evans, N. J. II, et al. 2009, ApJS, 181, 321
  • Fuller, Williams & Sridharan (2005) Fuller, G. A., Williams, S. J., & Sridharan T. K. 2005, A&A, 424, 949
  • Galli et al. (2002) Galli, D., Walmsley, M., & Concalves, J. 2002, A&A, 394, 275
  • Gao, Lou & Wu (2009) Gao, Y., Lou, Y.-Q., & Wu, K. 2009, MNRAS, 400, 887
  • Gao & Lou (2010) Gao, Y., & Lou, Y.-Q. 2010, MNRAS, 403, 1919
  • Goldsmith (2001) Goldsmith, P. F. 2001, ApJ, 557, 736
  • Goldsmith & Langer (1978) Goldsmith, P. F., & Langer, W. D. 1978, ApJ, 222, 881
  • Gregersen et al. (2000) Gregersen, E. M., Evans, N. J. II, Mardones, D., & Myers, P. C. 2000, ApJ, 533, 440
  • Gregersen et al. (1997) Gregersen, E. M., Evans, N. J. II, Zhou, S., & Choi, M. 1997, ApJ, 484, 256
  • Harvey, Wilner & Myers (2003) Harvey, D. W. A., Wilner, D. J., & Myers, P. C. 2003, ApJ, 583, 809
  • Herbst & Klemperer (1973) Herbst, E., & Klemperer, W. 1973, ApJ, 185, 505
  • Ho & Townes (1983) Ho, P. T. P., & Townes, C. H. 1983, ARA&A, 21, 239
  • Hogerheijde & van der Tak (2000) Hogerheijde, M. R., & van der Tak, F. F. S. 2000, A&A, 362, 697
  • Holland et al. (1999) Holland, W. S., et al. 1999, MNRAS, 303, 659
  • Hunter (1977) Hunter, C. 1977, ApJ, 218, 834
  • Jijina, Myers & Adams (1999) Jijina, J., Myers, P. C., & Adams, F. C. 1999, ApJS, 125, 161
  • Jørgensen et al. (2004) Jørgensen, J. K., et al. 2004, A&A, 415, 1021
  • Kandori et al. (2005) Kandori, R., et al. 2005, AJ, 130, 2166
  • Kenyon & Hartmann (1995) Kenyon, S. J., & Hartmann, L. W. 1995, ApJS, 101, 117
  • Keto et al. (2006) Keto, E., Broderick, A. E., Lada, C. J., & Narayan, R. 2006, ApJ, 652, 1366
  • Kirk, Ward-Thompson & André (2005) Kirk, J. M., Ward-Thompson, D., & André, P. 2005, MNRAS, 360, 1506
  • Kirk, Ward-Thompson & Crutcher (2006) Kirk, J. M., Ward-Thompson, D., & Crutcher, R. M. 2006, MNRAS, 369, 1445
  • Lada et al. (2003) Lada, C. J., Bergin, E. A., Alves, J. F., & Huard, T. L. 2003, ApJ, 586, 286
  • Landau & Lifshitz (1959) Landau, L. D., & Lifshitz, E. M. 1959, Fluid Mechanics, Springer
  • Larson (1969) Larson, R. B. 1969, MNRAS, 145, 271
  • Larson (1981) Larson, R. B. 1981, MNRAS, 194, 809
  • Lee, Myers & Tafalla (1999) Lee, C. F., Myers, P. C., & Tafalla, M. 1999, ApJ, 526, 788
  • Lee & Myers (2011) Lee, C. W., & Myers, P. C. ApJ, in press, 2011arXiv1104.2950L
  • Lou & Shen (2004) Lou, Y.-Q., & Shen, Y. 2004, MNRAS, 348, 717
  • Lou & Gao (2011) Lou, Y.-Q., & Gao, Y. 2011, MNRAS, 412, 1755