DEUTERIUM FRACTIONATION AS AN EVOLUTIONARY PROBE IN MASSIVE PROTOSTELLAR/CLUSTER CORES

In some cores, the emission of N₂D⁺ is undesirably weak and it becomes more challenging to constrain the spectral fit parameters. Instead of performing an independent spectral fit for N₂D⁺, we improve the model described in Chen et al. (2010) to perform a self-consistent fit for the N₂H⁺ and N₂D⁺ spectra using D_frac as a scaling factor. Because of the numerous hyperfine components, we fit each spectrum of every core with a synthetic spectrum comprised of 38 hyperfine components with updated line frequencies and spontaneous emission rates (Pagani et al. 2009). For each individual source, all the hyperfine components of every J-level are assumed to be in thermal equilibrium at a single excitation temperature, T_g, adopted from the ammonia gas temperature in Sakai et al. (2008). The synthetic spectra are described by three more parameters: the total column density, N(N₂H⁺), the systemic velocity, υ_LSR, and the FWHM as line width, Δυ. Model spectra are optimized with the minimization of the reduced χ² value, $\overline{\chi ^2}$, and the results are listed in Table 2. Note that the uncertainties of the adopted T_g (Sakai et al. 2008) are also incorporated into the uncertainties of all the derived spectral parameters.

Table 2. Parameters of Spectral Fits for N₂H⁺ and N₂D⁺

Core	$\sigma _\mathrm{N_2H^+}$	$\sigma _\mathrm{N_2D^+}$	N(N2H+)	υLSR	Δυ	Dfrac	$\overline{\chi ^2}$	τmax
	(mK)	(mK)	(1012 cm−2)	(km s−1)	(km s−1)
G34.43+0.24
MM1	44	9.5	17 ± 2	57.550 ± 0.008	4.81 ± 0.06	(39 ± 6) × 10−4	3.7	0.77
MM2	20	7.7	13 ± 1	57.579 ± 0.005	4.54 ± 0.05	(80 ± 7) × 10−4	9.8	0.61
MM3	20	7.8	6.0 ± 0.6	59.66 ± 0.01	4.20 ± 0.03	0.007 ± 0.002	3.1	0.32
MM4	21	8.8	5.8 ± 0.5	57.56 ± 0.01	4.76 ± 0.03	0.022 ± 0.002	4.2	0.27
MM5	19	7.1	1.8 ± 0.3	57.89 ± 0.03	3.45 ± 0.06	0.033 ± 0.005	1.6	0.12
MM6	23	10.5	0.59 ± 0.09	58.27 ± 0.06	2.3 ± 0.1	0.11 ± 0.02	0.8	0.05
MM8	24	8.5	0.33 ± 0.04	57.2 ± 0.1	3.4 ± 0.3	0.00 ± 0.02	1.2	0.02
MM9	17	8.9	1.6 ± 0.2	58.89 ± 0.02	3.04 ± 0.06	0.058 ± 0.007	2.5	0.11
IRAS 18151−1208
MM1	25	10.6	3.9 ± 0.4	33.203 ± 0.008	2.94 ± 0.03	0.010 ± 0.002	6.0	0.27
MM2	22	7.9	6.6 ± 0.8	29.706 ± 0.007	4.01 ± 0.04	0.019 ± 0.001	3.9	0.36
MM3	20	8.3	1.1 ± 0.3	30.67 ± 0.02	2.40 ± 0.06	0.064 ± 0.007	2.4	0.10
IRAS 18223−1243
MM1	30	8.6	3.5 ± 0.4	45.41 ± 0.02	3.69 ± 0.05	0.021 ± 0.002	3.4	0.20
MM2	38	7.9	1.9 ± 0.3	45.20 ± 0.05	3.7 ± 0.1	0.033 ± 0.005	1.1	0.11
MM3	40	8.6	4.2 ± 0.4	45.54 ± 0.03	4.42 ± 0.07	0.021 ± 0.002	1.5	0.21
MM4	43	8.4	0.9 ± 0.2	45.93 ± 0.07	2.5 ± 0.2	0.015 ± 0.009	1.2	0.08

Download table as:  ASCIITypeset image

Taking the cosmic background temperature, T_bg = 2.7 K, into account, we derive the optical depth of an observed spectrum with

$$\begin{equation} \tau (\upsilon) = {-}\!\ln \left[ 1 - \frac{T_\mathrm{mb}(\upsilon)}{J(T_g) - J(T_\mathrm{bg})} \right], \end{equation} \tag{ 2 }$$

where T_mb(υ) is the main beam temperature of the spectra and J(T) ≡ (hν/k)/(e^hν/kT − 1). In general, the spectra have fairly small optical depths if a beam filling factor of unity is assumed. The optimized spectral model also delivers an estimate for the N₂H⁺ optical depth by integrating optical depths of all the hyperfine components. The maximum optical depth, τ_max, is about 0.8 in G34.43–MM1 (Table 2). Emissions of N₂H⁺ and N₂D⁺ in all the observed cores remain optically thin. Although sub-structures within our observing beam sizes cannot be completely ruled out, the emission in these early phases of core evolution is mostly attributed to large scales. In the example of G28.34+0.06, the integrated flux observed with the Submillimeter Array is less than 10% of the total flux observed with the single-dish telescope SMT (Chen et al. 2010).

The CO depletion factor, f_D, is defined as the ratio of the canonical CO abundance, x(CO)_can, to the observed CO abundance, x(CO)_obs,

$$\begin{equation} f_D \equiv \frac{x(\mathrm{CO})_\mathrm{can}}{x(\mathrm{CO})_\mathrm{obs}} = \frac{x(\mathrm{C^{18}O})_\mathrm{can}}{x(\mathrm{C^{18}O})_\mathrm{obs}}. \end{equation} \tag{ 3 }$$

Using the C¹⁸O abundance of 1.7 × 10⁻⁷ in the solar neighborhood (Frerking et al. 1982) and the abundance gradients of Δlog [C/H]/ΔR = −0.066 dex kpc⁻¹ and Δlog [O/H]/ΔR = −0.065 dex kpc⁻¹ in the Galactic disk (Wilson & Matteucci 1992), the canonical abundance of C¹⁸O is estimated to be

$$\begin{equation} x(\mathrm{C^{18}O})_\mathrm{can} = 1.7 \times 10^{-7} \times 10^{-0.131 \, (D_\mathrm{GC}-D_\odot)}, \end{equation} \tag{ 4 }$$

where D_GC is the Galactocentric distance of the core and D_☉ = 8.5 kpc is the distance of the Sun to the Galactic center. Given the location of our IRDCs, x(C¹⁸O)_can = 3.80 × 10⁻⁷, 3.92 × 10⁻⁷, and 4.68 × 10⁻⁷ for G34.43, I18151, and I18223, respectively.

Overall, the observed brightness temperature, T_mb, is smaller than the kinetic temperature, T_g, and renders fairly small optical depths. The maximum optical depth estimated with Equation (2) in individual IRDC is 0.5 for G34.43–MM2, 0.5 for I18151–MM3, and 0.8 for I18223–MM1. Assuming that all rotational levels are thermalized, we determine the column density of C¹⁸O with the method based on Caselli et al. (2002) that accommodates the effect of background emission at T_bg,

$$\begin{equation} N(\mathrm{C^{18}O}) = \frac{3 h}{8 \pi ^3} \frac{1}{\mu ^2 S_{21}} \, W(\mathrm{C^{18}O}) \, \frac{Q(T_g)}{T_g - T_\mathrm{bg}} \, \frac{e^{E_\mathrm{up}/k T_g}}{(e^{h\nu _0/k T_g} -1)}, \end{equation} \tag{ 5 }$$

where Q(T_g) is the partition function, μ² S₂₁ = 0.02440 Debye² for the J = 2–1 transition, ν₀ is the transition rest frequency, E_up = 15.8 K is the upper level energy, and W(C¹⁸O) is the integrated brightness temperature in velocity (Table 1). Since some of the spectra do not resemble a Gaussian profile, direct integration in the channels with significant emission is performed to obtain W(C¹⁸O) without fitting a Gaussian profile, and the derived N(C¹⁸O) is listed in Table 3.

Table 3. Parameters of C¹⁸O Spectra and Derived Quantities

Core	$\sigma _\mathrm{C^{18}O}$	W(C18O)	N(C18O)	Tda	S1.2 mm	fD	log xe
	(mK)	(K km s−1)	(1015 cm−2)	(K)	(Jy/(34'')2)
G34.43+0.24
MM1	51	10.1 ± 0.1	4.71 ± 0.09	34	6.17	19 ± 4	−4.98
MM2	66	23.8 ± 0.1	11.2 ± 0.2	>34b	4.97	<6 ± 1	−5.30
MM3	63	5.1 ± 0.1	2.32 ± 0.06	32	1.35	9 ± 2	−5.27
MM4	77	13.0 ± 0.1	6.0 ± 0.1	32	2.74	7 ± 1	−5.77
MM5	57	6.64 ± 0.08	3.0 ± 0.1	...	0.92	14 ± 3	−6.01
MM6	52	5.37 ± 0.09	2.44 ± 0.09	...	0.67	13 ± 3	−6.62
MM8	55	5.2 ± 0.1	2.42 ± 0.07	...	0.72	10 ± 2	...
MM9	54	4.68 ± 0.08	2.12 ± 0.07	...	0.75	16 ± 3	−6.28
IRAS 18151−1208
MM1	54	16.13 ± 0.09	7.8 ± 0.2	27.3	2.41	6 ± 1	−5.37
MM2	59	12.0 ± 0.1	5.6 ± 0.1	19.4	1.85	10 ± 2	−5.71
MM3	49	8.33 ± 0.08	3.8 ± 0.2	...	0.81	8 ± 2	−6.33
IRAS 18223−1243
MM1	69	20.3 ± 0.1	9.4 ± 0.2	31	1.74	3.7 ± 0.7	−5.80
MM2	60	10.10 ± 0.09	4.6 ± 0.1	...	0.54	6 ± 1	−6.03
MM3	72	8.9 ± 0.1	4.1 ± 0.1	18	0.71	7 ± 1	−5.77
MM4	73	4.6 ± 0.1	2.09 ± 0.09	...	0.24	6 ± 1	−5.65

Notes. ^aDust temperature based on SED fits in literature (Rathborne et al. 2005; Marseille et al. 2008; Beuther et al. 2010). ^bDust temperature adopted from G34.43–MM1.

Download table as:  ASCIITypeset image

The observed C¹⁸O fractional abundance, x(C¹⁸O)_obs, depends on the column density of H₂, $N_\mathrm{H_2}$, and $x(\mathrm{C^{18}O})_\mathrm{obs} = N(\mathrm{C^{18}O})_\mathrm{obs}/N_\mathrm{H_2}$. In the attempt to reduce the uncertainty in deriving x(C¹⁸O)_obs, we first match the angular resolutions between the C¹⁸O and dust continuum observations by convolving the 1.2 mm continuum maps in the literature (beam size = 11''; Beuther et al. 2002a; Rathborne et al. 2006) with the 34'' beam of our C¹⁸O observations. The column density of H₂ is estimated with the 1.2 mm peak flux density, S_1.2 mm, arising from warm dust

$$\begin{equation} N_\mathrm{H_2} = \frac{S_\mathrm{1.2 \, mm}}{\kappa _\mathrm{1.2 \, mm} \, B_\nu (T_d) \, \Omega _b} \frac{1}{\mu m_\mathrm{H_2}}, \end{equation} \tag{ 6 }$$

where κ_1.2 mm = 0.005 cm² g⁻¹ is the dust opacity assuming a gas-to-dust mass ratio of 100 (Shepherd & Watson 2002), B_ν(T_d) is the Planck function at dust temperature, T_d, Ω_b is the solid angle subtended by the convolved beam size of 34'', μ = 1.36 is the mean molecular weight, and $m_\mathrm{H_2}$ is the mass of H₂ molecule.

For warmer cores in our sample, the averaged dust temperature, T_d, has been found from studies of the spectral energy distributions (SEDs; Table 3; Rathborne et al. 2005; Marseille et al. 2008; Beuther et al. 2010). However, most cores in our sample remain as extinction features in the near- or mid-IR, and it is challenging to derive their dust temperatures. Sensitive mid- and far-IR observations, such as Herschel, will offer better opportunities to constrain dust emission properties, including temperature, in cold clumps of IRDCs (e.g., Beuther et al. 2010; Stutz et al. 2010). Alternatively, we adopt the gas temperature, T_g, for the conversion between S_1.2 mm to $N_\mathrm{H_2}$ when dust temperature is unavailable. In case the thermal coupling between dust and gas is not ideal, this assumption may underestimate T_d and hence overestimate N(H₂) by a factor of <2 in our sample. In the case of G34.43–MM2, no SED study is found to constrain its dust temperature. Since the core is associated with a UC H ii region, we expect its averaged dust temperature to be warmer with respect to MM1 and sets a lower limit of T_d ⩾ 34 K for MM2.

Once the observed C¹⁸O abundance, x(C¹⁸O)_obs, is determined, the CO depletion factor, f_D, can be computed with Equation (3) accordingly, and the results are summarized in Table 3.

Overall, a monotonically decreasing trend in D_frac with both increasing gas temperature, T_g, and fitted line width, Δv, is discerned (Figures 4(a) and (b)). While examining cores in individual clouds, there seems to be a slightly better correlation, implying a possible cloud-to-cloud variation due to the influence of environments as previously suggested in the studies of low-mass cores (Crapsi et al. 2005; Emprechtinger et al. 2009). Miettinen et al. (2011) also suggest relatively large environmental variations in the cosmic-ray ionization rates in massive IRDC cores. In particular, quiescent cores with no outflow activities, i.e., G34.43–MM6, G34.43–MM9, I18151–MM3, and I18223–MM2, all have the lowest temperature (Table 1) and the largest D_frac (Table 2) as well. Since warm and turbulent gas is more likely to be associated with evolved cores, the observed D_frac suggests a decreasing dependence with evolutionary stage. For a better determination of their evolutionary stage, one may desire to compare with theoretical evolutionary models, which depend on several physical parameters, such as bolometric temperature, bolometric luminosity, and envelope mass (Froebrich 2005). However, most of our selected cores, especially those with the largest D_frac, do not have observations in the mid- and far-IR wavebands to better constrain their dust temperature and bolometric luminosity. Alternatively, we compare D_frac with gas temperature and line width instead of bolometric temperature and luminosity.

A few previous studies have reported an anti-correlation between deuterium fractionation and evolutionary stage of massive protostellar cores. Chen et al. (2010) found D_frac = 0.017–0.052 in three cores at different evolutionary stages within the IRDC G28.34+0.06, including a massive starless core (MM9), a core with fragmentation and outflow activities (MM4; Wang et al. 2011), and a UC H ii region (MM1), and suggested a decreasing trend in D_frac with evolutionary stage. A subsequent study by Fontani et al. (2011) significantly improved the statistics with a sample of 27 cores and also revealed this decreasing trend with values of D_frac in the range of 0.012–0.7, 0.017 − ⩽0.4, and 0.017–0.08 for their observed HMSCs, HMPOs, and UC H ii regions, respectively. In particular, they found an anti-correlation between D_frac and T_g. Miettinen et al. (2011) further reported a decreasing trend in the ratio of N(DCO⁺)/N(HCO⁺) = 0.0002–0.014 with gas temperature in their sample of seven IRDC cores. They also reported D_frac = 0.002–0.028 for the four cores with higher gas temperature. Our values of D_frac are in the range of 0.0039–0.11, which are comparable to the values obtained by Fontani et al. (2011) and Miettinen et al. (2011). The anti-correlation between D_frac and T_g is also seen in our results (Figure 4(a)).

Given the gas temperature range of T_g = 14–21 K, the corresponding thermal line width of N₂H⁺ is merely Δv_th = 0.15–0.18 km s⁻¹. The observed line width in the range of 2.3–4.8 km s⁻¹ is dominated by nonthermal motions, possibly arising from turbulent motions among clumps. In Figure 5, we compare the line width of our N₂H⁺ (3–2) spectra with line widths of N₂H⁺ (1–0) and NH₃ (1, 1), (2, 2), and (3, 3) spectra observed by Sakai et al. (2008). The N₂H⁺ (1–0) observations had a smaller beam of 18'' while the NH₃ observations had a much larger beam of 73'' with respect to our 27'' beam for N₂H⁺ (3–2). Between the two transitions of N₂H⁺, the J = 3–2 transition with higher E_up = 26.8 K shows a broader line width than the J = 1–0 with lower E_up = 4.5 K, as warmer gas is expected to be more turbulent (Figure 5(a)). In general, NH₃ lines have a much lower critical density of n_crit ∼ 2 × 10³ cm⁻³ (Evans 1999) and tend to trace the outer part of the cores with respect to N₂H⁺ lines with n_crit ∼ 10⁵ cm⁻³. Except in the two quiescent cores, G34.43–MM6 and I18151–MM3, line widths of the NH₃ (1, 1) and (2, 2) spectra are comparable to those of N₂H⁺ (1–0) but smaller than those of N₂H⁺ (3–2) (Figures 5(b) and (c)). The two quiescent cores are probably in a very early stage where turbulence dissipation may occur to produce smaller line width in the inner part of higher density (Goodman et al. 1998). In a number of more evolved cores, the NH₃ (3, 3) emissions with much higher E_up = 124.5 K are detected and show much larger line width, even up to 7.2 km s⁻¹ in the case of G34.43–MM3 (Figure 5(d)). Since all these evolved cores are associated with outflow activities, the NH₃ (3, 3) emission is tracing the hot gas which could be in the outflows (Zhang et al. 2002).

Zoom In Zoom Out Reset image size

Unlike NH₃ and NH₂D, which are affected by grain surface reactions (Gürtler et al. 2002; Bottinelli et al. 2010), N₂D⁺ and N₂H⁺ are pure gas-phase reactants and do not participate condensation and subsequent sublimation of ice mantles. Compared to other molecular species, the deuterium fractionation of N₂H⁺ better reflects the physical conditions at the present time without being confused by evaporation of mantles which had formed at earlier times with enhanced deuteration (Emprechtinger et al. 2009). In a sample of Taurus cores, a clear increasing trend in the deuterium fractionations of both NH₃ and N₂H⁺ was observed in prestellar cores (Crapsi et al. 2005; Hatchell 2003), whereas in protostellar cores, D_frac shows a faster decreasing trend with dust temperature than does the N(NH₂D)/N(NH₃) ratio (Emprechtinger et al. 2009). For high-mass protostellar cores within a similar gas temperature range, we note that the correlation between D_frac and T_g appears stronger in our cores than does the N(NH₂D)/N(NH₃) ratio (Pillai et al. 2007).

To further examine the relationship between deuterium fractionation and the CO depletion, we computed the C¹⁸O column density and the CO depletion factor, f_D, as described in Section 3.2. The results are listed in Table 3 and shown in Figure 4(c). When excluding G34.43–MM1, which has the largest value of f_D, we find an increasing trend in D_frac with f_D. This agrees with the general expectations from chemical models that CO is the major destroyer of H⁺₃, H₂D⁺, N₂H⁺, and N₂D⁺ (Caselli et al. 1998). As the envelope heats up, the CO abundance is expected to quickly rise based on the dramatical drop of the CO sublimation timescale from 10⁸ yr at T_d ≃ 12 K to 0.1 yr at T_d ≃ 20 K (Collings et al. 2003). The warmer temperature together with the return of gas-phase CO can alter the competition in the chemical networks and brings a drop in D_frac (Roueff et al. 2005; Aikawa et al. 2005; Aikawa 2008). When CO returns to the gas phase, it will quickly react with N₂H⁺ and N₂D⁺ to form HCO⁺ and DCO⁺, respectively (Lee et al. 2004). However, one should be cautious about possible chemical stratification once the CO sublimation starts in the central warm region. As a core warms up, the N₂H⁺ and N₂D⁺ emissions tend to trace the cold outer region whereas the CO emission is dominated in the central region. In the study of the Ophiuchus B2 core, Friesen et al. (2010) found that D_frac increases at greater projected distances from the embedded sources. When the observed emission arises from a partially filled volume, the beam-averaged abundance for each species will start to deviate from the actual abundance used in chemical models. Observations with high angular resolutions will help to image the spatial distributions and reduce the confusion in comparing abundances.

In the case of G34.43–MM1, an unusually large value of $N_\mathrm{H_2}$ and hence f_D is obtained. Millimeter interferometric studies (Cortes et al. 2008; Rathborne et al. 2008) reveal a very strong (L_bol ∼ 2 × 10⁴ L_☉) and compact (2R ≃ 0.03 pc) source that exhibits signatures of an HMC, which is thought to have a typical temperature closer to 100 K. This source has developed a steep temperature gradient, from 34 to 100 K across core scales from 0.1 to 0.015 pc (Rathborne et al. 2008), translating to a single power-law dependence of r^−0.57. This steep temperature gradient suggests the presence of an inner region where optical depth is large to the IR photons that the photon diffusion should be considered (Kenyon et al. 1993; Osorio et al. 1999; Chen et al. 2006). In our approach to estimate f_D, the dust emission is likely dominated by the hot inner region while the CO emission arises from a cold and large envelope. The dust temperature, T_d, derived from the SED fit depends on the flux densities in the mid-IR wavebands that may suffer from significant optical depth and does not reflect the physical conditions of the hot central region where the peak of the optically thin 1.2 mm emission is produced. Such a steep temperature gradient may cause an underestimate in T_d and overestimates in $N_\mathrm{H_2}$ and f_D in our current calculation. To avoid potentially misleading interpretation, we exclude the result of G34.43–MM1 when discussing the electron abundance in the following section.

The enrichment of primary deuterated ions, e.g., H₂D⁺, CH₂D⁺, and C₂HD⁺, will give rise to the enrichment of subsequent deuterated species, such as N₂D⁺, but with lower [D]/[H] abundance ratios due to the statistical nature of the fractionation process. In simple steady-state models based on gas-phase ion–molecular chemistry, an upper limit to the electron abundance, x_e, can be found assuming that all the deuterium enrichment originates in H₂D⁺ and that the recombination on negatively charged grains is negligible (Wootten et al. 1979; Caselli et al. 1998; Caselli 2002). Following the method described in Caselli (2002), the deuterium fractionation, D_frac, may be expressed as a function of x_e and abundances of HD, x(HD), and important neutral species, x(m),

$$\begin{equation} D_\mathrm{frac} = \frac{1}{3} \frac{k_\mathrm{HD} \, x(\mathrm{HD})}{k_e x_e + \displaystyle \sum _m k_m \, x(m) }, \end{equation} \tag{ 7 }$$

where k_HD = 1.5 × 10⁻⁹ cm³ s⁻¹ is the rate coefficient for the reaction in Equation (1), k_e = 6 × 10⁻⁸ (T/300)^−0.65 cm³ s⁻¹ is the dissociated recombination rate of H₂D⁺ (Caselli et al. 1998), and k_m is the destruction rate for H₂D⁺ due to reactions with neutral species m, such as CO and O. The numerical factor of 1/3 accounts for the statistical branching ratio of 1/3 to transfer the deuteron in the reaction of H₂D⁺ with N₂. The HD abundance is taken from the interstellar value of x(HD) = 2[D]/[H] = 3 × 10⁻⁵ (Oliveira et al. 2003; Caselli et al. 1998). The electron abundance x_e is then given by

$$\begin{equation} x_e = \frac{k_\mathrm{HD} \, x_\mathrm{HD}}{3 \, k_e} \frac{1}{D_\mathrm{frac}} - \frac{1}{k_e} \sum _m k_m x(m). \end{equation} \tag{ 8 }$$

Since CO is the dominant neutral species that destroys H₂D⁺, we make an approximation by neglecting other ion–neutral reactions to get

$$\begin{equation} \sum _m k_m x(m) \gtrsim k_\mathrm{CO} \, x(\mathrm{CO}) = k_\mathrm{CO} \, \left(\frac{x(\mathrm{CO})_\mathrm{can}}{f_D} \right), \end{equation} \tag{ 9 }$$

where k_CO = 6 × 10⁻¹⁰ (T/300)^−0.5 cm³ s⁻¹ is the H₂D⁺ destruction rate due to reactions with CO (Caselli et al. 1998), and x(CO)_can = 1.5 × 10⁻⁴ is the canonical CO abundance at the locations of our cores. We may set an upper limit for x_e to be

$$\begin{eqnarray} x_e &\lesssim & \frac{k_\mathrm{HD} \, x_\mathrm{HD}}{3 \, k_e} \frac{1}{D_\mathrm{frac}} - \frac{k_\mathrm{CO} \, x(\mathrm{CO})_\mathrm{can}}{k_e} \frac{1}{f_D}\qquad\qquad \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} &=& \frac{3.8 \times 10^{-8}}{D_\mathrm{frac}} - \frac{9.7 \times 10^{-7}}{f_D}\;\;\; (\mathrm{for}\; T_g = 16 \; \mathrm{K}). \end{eqnarray} \tag{ 11 }$$

For easy comprehension of the dependence on D_frac and f_D, the numerical values are provided for the case of T_g = 16 K. With the derived D_frac and f_D, we obtain the ionization degree in the range of x_e = 2 × 10⁻⁷ to 5 × 10⁻⁶ for our selected cores (Table 3). These values lie at the high end of the ionization degrees reported in early studies of low-mass dense cores (Caselli et al. 1998; Williams et al. 1998) and massive cores (Bergin et al. 1999). Recently, Miettinen et al. (2011) derived the first estimates for ionization degrees in IRDC cores with upper limits of x_e = 2 × 10⁻⁶ to 2.9 × 10⁻⁴ and lower limits of x_e = 3 × 10⁻⁹ to 5.6 × 10⁻⁸. Our estimates give smaller values compared to their upper limits but are within the range bracketed by their upper and lower limits. The smaller x_e upper limits are mainly attributed to the larger deuterium fractionation derived from N₂H⁺ instead of HCO⁺. Furthermore, our x_e values show a moderate correlation with the evolutionary stage. Since more evolved cores show smaller values of D_frac and f_D, the plausible correlation between x_e and the evolutionary stage may lead to the decreasing trend of D_frac. Most of the cores in our sample have shown outflow signatures and likely have begun to form clusters with low- and intermediate-mass protostars. This increasing degree of ionization could be arising from accretion shocks and/or heating of the gas related to the central protostellar/cluster objects (Stahler et al. 1980; Hosokawa et al. 2010; Calvet et al. 2004). Given the porous nature of the surrounding medium (Indebetouw et al. 2006), it is possible for part of the energetic photons to reach the outer part of the core and increase the overall electron abundance (Kim & Koo 2001).

The four quiescent cores, G34.43–MM6, G34.43–MM9, I18151–MM3, I18223–MM2, also have the lowest ionization degrees, x_e = 2 × 10⁻⁷ to 9 × 10⁻⁷, among cores in the same IRDC. In these cores, the degree of ionization may play a role in regulating star formation efficiency through ambipolar diffusion, in which magnetic fields drift relative to a background of neutrals (Mestel & Spitzer 1956; Shu et al. 1987). In a partially ionized medium, charged particles are coupled to magnetic fields while neutrals are supported against their self-gravity through the frictional drag that they experience when drifting through ions. In addition to ambipolar diffusion, evolution of massive star-forming cores can be affected by various effects such as turbulence (McKee & Tan 2003), rotation, and magnetic fields. For a more complete picture, one needs to consider all the important supports against gravity (Myers & Goodman 1988; McKee & Ostriker 2007). Here we just make an attempt to estimate the characteristic scale for ambipolar diffusion, ℓ_AD, which gives the smallest scale for which the magnetic field is well coupled to the bulk of the gas. Following the picture described in McKee & Ostriker (2007) and approximating the ion abundance with $n_i = x_e \, n_\mathrm{H_2}$, we estimate the characteristic ambipolar diffusion scale to be

$$\begin{eqnarray} \ell _\mathrm{AD} &=& \frac{v_A}{n_i \, \alpha _{{\rm in}}} \approx 0.05 \; \mathrm{pc} \left(\frac{v_A}{3 \; \mathrm{km \, s^{-1}}} \right) \left(\frac{n_i}{10^{-3} \; \mathrm{cm^{-3}}} \right)^{-1} \nonumber \\ &=& 0.002 \; \mathrm{pc} \left(\frac{B}{10 \; \mathrm{\mu G}} \right) \left(\frac{x_e}{10^{-7}} \right)^{-1} \left(\frac{n_\mathrm{H_2}}{10^4 \; \mathrm{cm^{-3}}} \right)^{-3/2},\qquad \end{eqnarray} \tag{ 12 }$$

where $v_A = B/\sqrt{4 \pi \, \mu m_\mathrm{H_2} \, n_\mathrm{H_2} }$ is the Alfvén speed associated with the large-scale magnetic field B and α_in ≈ 2 × 10⁻⁹ cm³ s⁻¹ is the ion–neutral collision rate coefficient (Draine et al. 1983). For the unknown magnetic field strength in supersonic cores, we take the median value given by McKee & Ostriker (2007; originally from Crutcher 1999),

$$\begin{equation} B_\mathrm{med} = 30 \left(\frac{2 \, n_\mathrm{H_2}}{10^3 \; \mathrm{cm^{-3}}} \right)^{1/2} \left(\frac{\Delta v_\mathrm{nt}/\sqrt{8 \ln 2}}{1 \; \mathrm{km \, s^{-1}}} \right) \; \mathrm{\mu G}, \end{equation} \tag{ 13 }$$

where $\Delta v_{\rm nt} = \sqrt{\Delta v^2 - \Delta v_\mathrm{th}^2}$ is the line width for nonthermal motions. Applying the values for $n_\mathrm{H_2}$ and Δv in the four quiescent cores, we find the magnetic field strength to be in a range from 210 to 572 μG. The characteristic scale for ambipolar diffusion is in a range from ℓ_AD = 6 × 10⁻⁴ to 4 × 10⁻³ pc, much smaller than the typical size of the cores.

Since most star-forming cores have Alfvén speeds comparable to the nonthermal component of the velocity dispersions, the turbulence in these cores should have a substantial magnetohydrodynamic (MHD) component (Myers & Goodman 1988; Crutcher et al. 1999). Myers & Lazarian (1998) have proposed a pressure-driven cooling flow associated with local dissipation of turbulence due to wave damping by ion–neutral friction in an inner core region. From a sufficiently turbulent outer region, the MHD wave power transmission, g, into a spherical inner core is described in terms of the field-neutral coupling parameter, W, which is defined as the ratio of the core size, R, to the minimum cutoff wavelength, λ₀ = πℓ_AD, for the propagation of MHD waves (Myers & Khersonsky 1995). Using the upper limits of x_e in the four quiescent cores, we obtain W = 23–119 with a mean value of W = 75, roughly in the regime of marginal field-neutral coupling (g ∼ 0.4 for W ≲ 100; Myers & Lazarian 1998). The MHD waves excited in the outer region would have a fairly limited range of allowed wavelengths to transmit the turbulence power into the dissipating inner regions, and a pressure-driven turbulent cooling flow may occur in these cores. For comparison, Bergin et al. (1999) reported W = 20 in massive cores with ionization degree inferred from chemical models while Miettinen et al. (2011) obtained W ⩽ 18.5 using their lower limits of ionization degree. But one should be cautious about a direct comparison among these values because of different methods to estimate the ionization degrees and the magnetic field strength.