CORRELATING INFALL WITH DEUTERIUM FRACTIONATION IN DENSE CORES

Starless cores are dense (n ∼ 10⁵ cm⁻³), cold (T ∼ 10 K), and compact (r ∼ 0.05 pc) objects embedded within molecular clouds (Di Francesco et al. 2007). Starless cores are the transition phase between the comparatively diffuse molecular cloud material and a new generation of protostars. The cores are often identified from their (sub)millimeter dust continuum emission and classified as starless from the lack of infrared emission indicative of an embedded protostar (e.g., Nutter & Ward-Thompson 2007; Hatchell et al. 2007; Enoch et al. 2008). Recent (sub)millimeter interferometric observations have shown that low-luminosity protostars have been overlooked by Spitzer in roughly 10% of supposedly starless cores (Schnee et al. 2012).

In the cold and dense interiors of starless cores, chemical reactions lead to deuterium enrichment of some molecules that remain in the gas phase, such as N₂H⁺(see, e.g., Caselli et al. 2008). In a survey of low-mass cores in nearby molecular clouds, Crapsi et al. (2005) measured N(N₂D⁺)/N(N₂H⁺) values ranging from 0.01 to 0.44, much larger than the [D]/[H] elemental abundance of ≃ 1.5 × 10⁻⁵ (Oliveira et al. 2003) found in the interstellar medium. Crapsi et al. (2005) concluded that the N(N₂D⁺)/N(N₂H⁺) ratio is an indicator of core evolution such that a high ratio implies that a core is close to transitioning from being starless to forming a protostar. Friesen et al. (2013) showed that the deuterium fraction in cores in the Perseus molecular cloud tends to increase with H₂ column density and with increasing core concentration, in agreement with the Crapsi results. Emprechtinger et al. (2009) found that the N₂D⁺/N₂H⁺ ratio is also an indicator of evolutionary state for Class 0 protostars, with the youngest protostars having the highest deuterium fractionation and clearest signs of ongoing collapse. In the Perseus sample, the most highly deuterated protostellar cores are still associated with significant envelope mass (Friesen et al. 2013). Therefore, in the evolution of a collapsing core, the N(N₂D⁺)/N(N₂H⁺) ratio is expected to peak just as the protostar is formed. In a sample of starless cores, one might expect the N₂D⁺/N₂H⁺ ratio to be roughly indicative of the likelihood of collapse, with scatter due to differences in core properties such as mass, temperature, non-thermal motions, and magnetic field strength.

Infall and outflow motions in dense cores are revealed by asymmetries in their spectral line profiles. In a simple model of collapse, density, optical depth, and excitation temperature increase toward the center of the core. Inward motions of the side of the core nearer to the observer induce a redshift relative to the average core velocity, and, as a result, one expects to see redshifted self-absorption of optically thick lines. Myers et al. (1996) and De Vries & Myers (2005) present sets of equations that can be used to derive infall velocity in idealized cases of core velocity and density profiles, and these models match observations with reasonable accuracy. In a real core, the inferred infall velocity depends on the point at which the observed line becomes optically thick and the variation of the infall velocity with radius, complicating the interpretation of the observed spectra. Nevertheless, surveys of infall motions in dense cores have been successfully carried out using various molecular tracers, such as HCO⁺ (Gregersen & Evans 2000), HCN (Sohn et al. 2007), CS (Williams & Myers 1999; Lee et al. 1999, 2001), and N₂H⁺ (Williams & Myers 1999; Lee et al. 1999, 2001; Crapsi et al. 2005).

Infall in a starless core will result in a localized density increase due to mass accumulation and a decrease in the central core temperature due to the additional shielding from the interstellar radiation field. These conditions promote increased deuteration, so one might expect to see large infall velocities in those cores with the highest degree of deuteration. To test the hypothesis of a correlation between chemical and dynamical evolution, we observed the HCO⁺ (3–2) and H¹³CO⁺ (3–2) lines toward every core in the Crapsi et al. (2005) sample for which the N(N₂D⁺)/N(N₂H⁺) ratio was measured. If those cores with the highest N(N₂D⁺)/N(N₂H⁺) ratio are really on the verge of forming a protostar, then we expect them to also show infall signatures. Those cores with smaller observed deuterium enrichment would have slower, or non-existent, infall motions.

Similarly, those cores more massive than their thermal Jeans mass might also be expected to exhibit infall motions. We also compare the observed motions with the Jeans stability of each core to test the hypothesis that those cores with M/M_J > 1 will have inward motions while cores with M/M_J < 1 will appear to be static. More details on this survey are reported below.

The 26 dense cores in this survey are those for which the N(N₂D⁺)/N(N₂H⁺) ratios were measured by Crapsi et al. (2005). The Crapsi et al. (2005) sample is meant to be representative of evolved (i.e., dense) starless cores, with their high densities derived from N₂H⁺ (1–0) or 1.2 mm dust continuum emission. Although thought to be starless by Crapsi et al. (2005), the cores L1521F (Bourke et al. 2006) and L328 (Lee et al. 2009) were later found to be protostellar. The Crapsi et al. (2005) cores are nearby, i.e., within 250 pc, making them comparatively easy to study. We chose to observe HCO⁺ (3–2) to trace inward and outward motions in these dense cores. HCO⁺ (3–2) was shown by Gregersen & Evans (2000) to be a good tracer of line asymmetries in starless cores. The effective critical density of HCO⁺ (3–2), 6.3 × 10⁴ cm⁻³ (Evans 1999), is well-matched to the densities of the cores in the Crapsi et al. (2005) sample, 8 × 10⁴ cm⁻³ to 1.4 × 10⁶ cm⁻³. The observed spectra are likely to be affected by optical depth and depletion such that HCO⁺ (3–2) may not be tracing exactly the same portion of every core in our sample. Although this source of uncertainty is unavoidable in a survey such as ours, we outline a possible path around this uncertainty in Section 5.

Our observations were centered on the position of the peak N(N₂D⁺)/N(N₂H⁺) given in Crapsi et al. (2005). The [N₂D⁺]/[N₂H⁺] ratio is expected to peak in the region with the highest density and coldest temperature. High densities and low temperatures (and therefore low thermal support) are conditions that are likely to lead to gravitational collapse. The coordinates of our pointed observations and the degree of deuteration are given in Tables 1 and 2. The noise values of the spectra are given in Table 1. The observed HCO⁺ (3–2) spectra are shown in Figures 1(a)–(d).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The data were obtained at the James Clerk Maxwell Telescope⁷ (JCMT) using the 1 mm Receiver RxA3 front-end and the ACSIS backend throughout semesters 10A and 10B. Each core was observed in position-switching mode with an off-position located (−300,−300) arcsec (in the J2000 coordinate frame) from the core position (see Table 1). ACSIS was configured to cover two 250 MHz wide windows, one centered on HCO⁺ (3–2) at 267.558 GHz and another centered on H¹³CO⁺ (3–2) at 260.255 GHz. Each window was divided into 8192 channels for a resolution of ∼0.05 km s⁻¹. Each core was observed for about 45 minutes of integration and the beam size of the observations was ∼20'' FWHM. The JCMT beam size is fairly well matched to the beam sizes of the N₂H⁺ (1–0) and N₂D⁺ (2–1) observations presented in Crapsi et al. (2005), i.e., 26'' and 16'' FWHM, respectively.

After acquisition, the data were reduced using the Starlink package. Each integration was first visually checked for baseline ripples and extremely large spikes, and such data were removed from the ensemble, if necessary. Spectral baselines were subtracted, frequency axes converted to velocities, and spectra edges trimmed using Starlink scripts kindly provided by T. van Kempen. The integrations were visually inspected and written out as FITS format files using other standard Starlink routines. The spectra for each core were co-added outside of Starlink, using the IDL programming language. Spectra were converted to main beam brightness temperatures using an assumed efficiency of 0.75.⁸ HCO⁺(3–2) emission was detected toward each core except for L1495A-S, but H¹³CO⁺ (3–2) emission was only detected toward a few cores (L1521F, Oph D, L1689B, and L328; see Figures 1(a)–(d)). Hence, the H¹³CO⁺ data will not be discussed further in this paper.

The HCO⁺ spectra were fit with a Gaussian model as well as two simple infall/outflow models, the "two-layer" model from Myers et al. (1996) and the "HILL5" model from De Vries & Myers (2005). The two-layer model describes the case of two regions with differing excitation temperatures (lower for the front, higher for the rear) and identical velocity dispersions and total optical depth. The HILL5 model describes the case of a core whose excitation temperature peaks at the center and falls off linearly with radius. The HILL5 core also has its front and rear halves moving toward or away from each other. Both models reproduce spectra with one peak and a shoulder reasonably well, though the HILL5 model is better than the two-layer model at fitting spectra with two peaks (De Vries & Myers 2005). The five free parameters in both infall/outflow models are inward/outward velocity, kinetic temperature, peak optical depth, velocity dispersion, and velocity with respect to the local standard of rest. The dynamic models assume that the two peaks in the prototypical infall/outflow profile arise from self-absorption and not from having two unrelated sources with slightly different velocities along the same line of sight. This assumption can be tested through observations of an optically thin tracer that should peak at the same velocity as the self-absorption dip in the optically thick spectrum if no radial motion is present. The H¹³CO⁺ (3–2) observations described in Section 2 were intended to be this optically thin tracer, but detections were rare. The N₂H⁺ and N₂D⁺ observations presented in Crapsi et al. (2005) are sufficient to show that there are no confounding superpositions in our sample. In the relatively few cases where H¹³CO⁺ (3–2) is detected, the line peaks at the same velocity as the N₂H⁺ and N₂D⁺ spectra. The Gaussian and dynamic models were fit with the MPFIT suite of functions (Markwardt 2009).

For each core, the Gaussian and dynamic model fits were compared with F-tests, taking into account the χ² values for the fits and the degrees of freedom in the fits. For those cores where the F-test prefers a dynamic model and the inward/outward motions are at least three times the uncertainty of the motion, we treat the core as having inward or outward motions for the rest of this paper. We describe as "static core" those cores for which the F-test prefers the Gaussian model or the uncertainty in the radial motions is greater than three times the magnitude of the motion. The signal-to-noise ratio of the L1495A-S observations is too low to attempt a fit to either model. The peak intensities, integrated intensities, and velocity dispersions (standard deviations, or σ for a Gaussian velocity distribution) measured from the best fits to the HCO⁺ spectra are given in Table 1, with the intensities given on T_mb scale. Note that the cores L1521F and L328 are actually protostellar, not starless, and L1521F has the largest infall velocity in our survey. The infall velocity of L1521F was measured to be 0.2–0.3 km s⁻¹ by Onishi et al. (1999). The presence of an embedded protostar likely reduces the accuracy of the HILL5 and two-layer models (note the relatively poor fit to the blue peak of the spectrum), so we trust the detailed model of Onishi et al. (1999) to provide a more accurate estimate of the infall velocity than we derive here. The kinematic state and preferred model of each core as determined by the model fits are given in Table 3, along with previous kinematic classifications from the literature (see Section 4). In Table 3, we also give the velocity with respect to the local standard of rest from the best-fit model to our HCO⁺ spectra, along with the LSR velocity fit to the N₂H⁺ (1–0) spectra from Crapsi et al. (2005). We find that the median of the absolute value of the velocity difference between the LSR velocity from HCO⁺ and N₂H⁺ is 0.07 km s⁻¹ with a standard deviation of 0.1 km s⁻¹. This velocity difference likely comes from velocity gradients along the line of sight due to inward/outward motions, since HCO⁺ and N₂H⁺ will sometimes probe different layers within a core. The uncertainties given in Table 3 are the 1σ standard deviations.

Masses and radii for the cores in this survey were taken from the literature and calculated from (sub)millimeter continuum surveys. Most core masses are derived from the observed flux at 850 μm or 1.2 mm using the formula

$$\begin{eqnarray} M &=& 0.12 [ {\rm e}^{\left(\frac{1}{1.439} \frac{\lambda }{\rm mm} \frac{T}{10\,{\rm K}} \right)^{-1}}-1 ] \left(\frac{\kappa _\nu }{0.01\,{\rm cm}^2\,{\rm g}^{-1}}\right)^{-1} \left(\frac{S_\nu }{{\rm Jy}}\right) \nonumber\\ && \times \left(\frac{D}{100\,{\rm pc}}\right)^2 \left(\frac{\lambda }{{\rm mm}}\right)^3 {{M}}_\odot. \end{eqnarray} \tag{ 1 }$$

We assume that all cores are isothermal at 10 K, as assumed in the core mass calculations of Crapsi et al. (2005) and consistent with the average temperatures derived from NH₃ observations of starless cores (Tafalla et al. 2002; Schnee et al. 2009). Real cores have temperature gradients ranging from ∼13 K at large radii to ∼6 K at the center (e.g., Crapsi et al. 2007; Pagani et al. 2007; Schnee et al. 2007b), but these variations in temperature do not affect the conclusions of this paper. We assume a value of κ_ν at 850 μm of 0.01 cm² g⁻¹ and κ_ν at 1.2 mm of 0.005 cm² g⁻¹, as assumed in Sadavoy et al. (2010) and Crapsi et al. (2005), respectively. These values for κ_ν are consistent with each other for an emissivity spectral index of 2, in agreement with observations of the starless cores L1498 (Shirley et al. 2005) and TMC-1C (Schnee et al. 2010). Distances (D) to the cores and references used to determine core masses are given in Table 2. We were unable to find previously published values for the masses and radii of the cores L1495 and L1400A.

Following Sadavoy et al. (2010), we calculate the Jeans mass for the cores in this survey using the formula

$$\begin{equation} M_J = 1.9 \left(\frac{T}{10\,{\rm K}}\right) \left(\frac{R}{0.07\,{\rm pc}}\right) {M}_\odot. \end{equation} \tag{ 2 }$$

As in the calculation of core masses, when calculating the Jeans mass, we assume that all cores are isothermal at T = 10 K and use the radii (R) given in Table 2.

A plot of infall velocity (see Table 3) versus degree of deuteration (see Table 2) is shown in Figure 2. Nine cores show no evidence for infall or outflow, 5 cores have outward motions, and 11 have inward motions. Below an N(N₂D⁺)/N(N₂H⁺) ratio of 0.1, there is no correlation with infall velocity. For the eight cores with N(N₂D⁺)/N(N₂H⁺) ⩾ 0.1, six cores have infall velocities and two cores appear to be static. We therefore find a general trend, though not a quantitative correlation, that inward motions are more likely to be found at higher levels of deuteration, whereas at lower levels of deuteration, inward and outward motions are roughly equally likely to be found. The HCO⁺ spectra therefore support the idea that the most chemically evolved starless cores are also dynamically evolved.

Zoom In Zoom Out Reset image size

Simpson et al. (2011) found that the Jeans stability of a pre-stellar core provides a reasonable predictor of its dynamical state. As shown in Figure 2 (see also Tables 2 and 3), we find that three Jeans-stable cores (those with M/M_Jeans < 1) show inward motions, two cores appear to be static, and four show outflow motions. Those cores with M/M_Jeans > 1 are much more likely to have line asymmetries associated with infall (eight cores) than outward motions (no cores), and six cores with M/M_Jeans > 1 are best fit by Gaussian profiles. Our results are in agreement with those of Simpson et al. (2011), in that the Jeans stability of dense cores is a useful indicator of likely infall candidates, and the identification of collapsing cores becomes more likely with increasing M/M_Jeans. As shown in Figure 2, those cores with N(N₂D⁺)/N(N₂H⁺) > 0.1 also have M/M_Jeans > 1.

Many cores with lower degrees of deuteration, however, also have high M/M_Jeans (L1495A-S, TMC-1C, CB 23, L1517B, L1689B, L492, and L328) and these cores sometimes exhibit infall motions that are a consequence of their Jeans instability. There are no cores with M/M_Jeans < 1 and N(N₂D⁺)/N(N₂H⁺) > 0.1. Those cores with M/M_Jeans > 1, but no indication of inward motions, may be supported by mechanisms other than thermal pressure (e.g., turbulent or magnetic pressure). Alternatively, infall may indeed be on-going, but the signature of this motion could be missing from the HCO⁺ (3–2) spectra in some cores due to freeze-out (Roberts et al. 2010) or rotation (Redman et al. 2004). We suggest that those cores with super-Jeans masses but no indication for infall are a good sample for future observations.

The origin of outward motions in five starless cores in our survey, as indicated by negative infall velocities in Table 3 and Figure 2, is not fully understood. Such outward motions have been seen in previous surveys (e.g., Sohn et al. 2007; Lee & Myers 2011) and may indicate oscillatory behavior (Lada et al. 2003; Broderick & Keto 2010). Indeed, oscillatory behavior may also be responsible for some of the inward motions observed in our HCO⁺ survey. Smoothed particle hydrodynamics simulations of nominally super-Jeans cores modeled as Bonnor–Ebert spheres may oscillate rather than collapse if the exterior temperatures of the cores are high enough (Anathpindika & Di Francesco 2013). As the sometimes-different classifications of cores as "expanding," "oscillating," or "contracting" in Table 3 make clear, observing cores at slightly different positions with different molecular tracers can lead to very different conclusions, an effect which surely adds considerable scatter to plots like Figure 2.

Only 3 of the 16 cores included in 4 surveys of inward and outward motions (Crapsi et al. 2005; Sohn et al. 2007; Lee & Myers 2011, this work) showed the same behavior in all cases (L1544, L694-2, and L1197, which all show inward motions). Infall in the cores L1544 and L694-2 has been particularly well-studied, with evidence for extended infall at speeds ⩽0.1 km s⁻¹ seen in N₂H⁺ observations (Tafalla et al. 1998; Williams et al. 2006; Keto & Caselli 2010). The lack of a 1:1 correlation between infall velocity and high N(N₂D⁺)/N(N₂H⁺) or M/M_Jeans is likely to be at least partly a result of the uncertainties in measuring the dynamics of cores with a single molecular tracer at a single position. In addition, because starless cores will not all reach exactly the same level of N(N₂D⁺)/N(N₂H⁺) just before collapse, one would not expect a perfect correlation with infall even in a comprehensive set of observations.

It may seem surprising that some cores in this survey have significant infall velocities without having similarly large N(N₂D⁺)/N(N₂H⁺) ratios. As mentioned above, in some cases, the inward motions may be due to oscillations in a stable core rather than being caused by collapse, in which case the low N(N₂D⁺)/N(N₂H⁺) ratio is a true indicator of the stability of a core. Alternatively, HCO⁺ may be preferentially tracing accretion motions (i.e., motions of material onto the core) instead of contraction or expansion motions of the core itself, as HCO⁺ is also very abundant at low densities (i.e., in the core envelope and the surrounding molecular cloud), unlike N₂H⁺. For example, molecular abundances in the contracting core TMC-1C are affected by accretion from its environment, from which it is gaining material rich in N₂H⁺ but likely not in N₂D⁺ (Schnee & Goodman 2005; Schnee et al. 2007a). It is also possible that cores with significant infall velocities but low levels of deuterium fractionation are too young to have significant deuterium enrichment despite their high central densities. For instance, the cores L1521E (Tafalla & Santiago 2004) and L1689B (Lee et al. 2003) have been identified as dense cores with little chemical evolution, likely due to their relative youth (⩽1.5 × 10⁵ yr). In these chemically young cores, where CO is not highly depleted, HCO⁺ may be a better tracer of infall in the central high density region, compared to chemically evolved cores. Therefore, a mismatch between the chemical and dynamical states of a starless core inferred by this survey could be due to (1) a single HCO⁺ (3–2) spectrum being an incomplete tracer of the dynamics of an oscillating dense core, (2) the N(N₂D⁺)/N(N₂H⁺) ratio within a core being lowered by accretion from the surrounding medium, or (3) the timescale of chemical evolution being significantly longer than the dynamical timescale.

The magnitude of the inward and outward velocities are low, similar to the sound speed (∼0.2 km s⁻¹ for a 10 K core). This result suggests that the initial conditions from which the motions are derived were near hydrostatic equilibrium. That is, the infall motions are not dominated by highly supersonic flows attributed to cloud turbulence but more likely derive from gravitational collapse of a core that has lost pressure support. The outward motions observed are consistent with either oscillations of a near gravitational equilibrium object (Broderick & Keto 2010; Anathpindika & Di Francesco 2013) or thermal expansion of an over-pressured but gravitationally unbound core, as would be expected for a core with M/M_Jeans ≲ 1.

As outlined in this paper, there is a need for a uniform multidimensional study of the kinematics and chemistry of dense cores. The work presented in this paper is limited by the single-pointing observations. Spectral line maps would be better able to determine the location of peak infall or outflow, and to determine whether a core shows infall or outflow spectral signatures consistently across its extent. Part of the difference between the infall surveys presented in Table 3 can be attributed to the different pointing centers of the observations. To the extent that the velocity field changes across a core, these positional differences will confuse the conclusions of pointed surveys. Furthermore, the "center" of a core is not always a well-defined location. For example the center of L183 is offset by ∼20'' between the studies of Crapsi et al. (2005) and Pagani et al. (2007). A second conclusion that can be drawn from Table 3 is that there is no single best molecular tracer of inward and outward motions in dense cores. Although HCO⁺ is relatively bright and shows infall motion clearly in many cores, it may not trace the regions in starless cores with the lowest temperatures, highest densities, and most depletion. In this case, the infall velocity measured from HCO⁺ spectra is likely to be an underestimate of the maximum infall speed. To map precisely the inward and/or outward motions in a core, one would want to have several spectral maps with different molecules and rotational transitions. These observations could then be modeled with a radiative transfer code to constrain self-consistently the density, temperature, chemistry, and kinematic profiles of the cores. Such a survey would require a substantial amount of observing time. We note that the pointed observations of a single infall tracer in this paper required over 30 hr of telescope time (including overheads) over multiple semesters with the JCMT. The combination of high sensitivity with high spatial and spectral resolution provided by the Atacama Large Millimeter/submillimeter Array may make it possible to carry out a comprehensive survey of the connection between the chemistry, kinematics, and role of environment in dense cores.

Crapsi et al. (2005) predicted that dense cores with N(N₂D⁺)/N(N₂H⁺) > 0.1 will show evidence for infall. To test this hypothesis, we observed HCO⁺ (3–2) toward the Crapsi et al. (2005) sample of 26 dense cores (24 starless and 2 protostellar) and analyzed the spectra to look for inward and outward motions. We find that those cores with the largest values of N(N₂D⁺)/N(N₂H⁺) and M/M_Jeans are indeed more likely to show the signature of infall than the population as a whole, though some cores that have lower levels of deuterium fractionation and greater Jeans stability also show inward motions. Since asymmetries from inward and outward motions vary by position within cores (Lee et al. 2001; Lada et al. 2003) and vary by tracer observed (e.g., Lee & Myers 2011), it is not surprising that we do not see very tight correlations between the degree of deuterium fractionation or Jeans stability with infall velocity. The presence of outward motions, and perhaps some of the inward motions, in these presumably long-lived cores is possibly a result of oscillations (Lada et al. 2003; Broderick & Keto 2010; Anathpindika & Di Francesco 2013) and such outward motions have been seen before in previous surveys of dense cores (e.g., Sohn et al. 2007; Lee & Myers 2011). The magnitudes of the inward and outward motions are on the order of the sound speed, suggesting that the cores formed near hydrostatic equilibrium rather than out of supersonic flows.

We thank Malcolm Currie for his help reducing JCMT spectra. We thank our anonymous referee for suggestions that have significantly improved the clarity and content of this paper. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. J.D.F. and D.J. acknowledge support by the National Research Council of Canada and the Natural Sciences and Engineering Council of Canada (NSERC) via a Discovery Grant. A.P. is partially supported by the NSERC graduate scholarship program. R.F. is a Dunlap Fellow at the Dunlap Institute for Astronomy and Astrophysics, University of Toronto. The Dunlap Institute is funded through an endowment established by the David Dunlap family and the University of Toronto. The James Clerk Maxwell Telescope is operated by the Joint Astronomy Centre on behalf of the Science and Technology Facilities Council of the United Kingdom, the Netherlands Organisation for Scientific Research, and the National Research Council of Canada. The JCMT program ID numbers associated with this project are M10AC002, M10BC001, and M11BC003. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. P.C. acknowledges the financial support of the European Research Council (ERC; project PALs 320620).

Facility: JCMT - James Clerk Maxwell Telescope