OBSERVATIONAL DETERMINATION OF THE TURBULENT AMBIPOLAR DIFFUSION SCALE AND MAGNETIC FIELD STRENGTH IN MOLECULAR CLOUDS

Interstellar magnetic fields and turbulence have been the most debated mechanisms of support against gravitational collapse in molecular clouds. Recent studies have evolved toward acknowledging both mechanisms by studying turbulence in a magnetized gas (e.g., Shebalin et al. 1983; Goldreich & Sridhar 1995; Ostriker et al. 2001; Biskamp 2003; Basu et al. 2009; Falceta-Gonçalves et al. 2010). It is, however, a challenge to combine the existing models for turbulence with observational data for the measurement of the interstellar magnetic fields to fully understand the dynamics in molecular clouds. There is yet no closed theory for magnetized turbulence due to its complicated and intermittent nature. On the other hand, interstellar magnetic fields are weak (∼few hundred μG to a few mG) and difficult to measure with the very few observational techniques currently at hand (e.g., Crutcher et al. 1999; Houde 2004; Falgarone et al. 2008; Hildebrand et al. 2009; Houde et al. 2009). It is therefore essential to innovate with more observing techniques to provide reliable data and test existing models.

Turbulence is basically a dissipative process. In an incompressible turbulent fluid, energy is transferred from larger eddies to smaller ones, until it is dissipated by viscosity at small enough scales (Kolmogorov 1941). The dissipation of turbulent energy has important effects on the interstellar medium (ISM), such as initiating collapse in supercritical cores (Basu et al. 2009), modifying the chemical evolution of the gas molecules (Xie et al. 1995), and causing the decoupling of gas and grains in a cloud (Falgarone & Puget 1995).

Li & Houde (2008) established a new observational technique, expanding on the original work of Houde et al. (2000a, 2000b) on the effect of magnetic fields on the spectral line widths of ion and neutral molecular species, to calculate the ambipolar diffusion length scale and the strength of the plane of the sky component of the ambient magnetic field (B_pos) in molecular clouds. They studied the turbulent velocity dispersion spectra of a pair of coexistent ion and neutral molecular species (HCN and HCO⁺) in M17, with the assumption that the turbulent energy dissipation process is associated with ambipolar diffusion.

In this paper, we explain the method of Li & Houde (2008) and our purpose for testing their turbulent energy dissipation model using optically thin coexistent ion and neutral molecular species in Section 2, and present our observations in Section 3. The data analysis is explained in Section 4 and we end with a summary in Section 5.

Ostriker et al. (2001) introduced a numerical method to study the kinetic structure in a single-fluid turbulent cloud by measuring the distribution of the velocity dispersion against the size and mass of the simulated regions over which it was averaged. They showed how the lower envelope of the calculated velocity dispersions depicts the actual turbulent velocity spectrum of the gas in their simulations. Li & Houde (2008) used this technique to study the velocity dispersions of HCN (J = 4 → 3) and HCO⁺ (J = 4 → 3) spectral line profiles in M17 as a function of length scale. They noticed that the lower envelope of the HCO⁺ data was consistently downshifted from that of HCN by a constant amount, independent of scale, and interpreted this result as the signature of turbulent ambipolar diffusion acting at smaller spatial scales than probed by their observations. Based on these results, Li & Houde (2008) suggested a scenario where ambipolar diffusion sets in at small enough length scales where neutrals decouple from the flux-frozen ions and the magnetic field, and the friction between the ions and the drifting neutrals causes the dissipation of turbulent energy. Tilley & Balsara (2010) tested this model by simulating a two-fluid (ion and neutral) turbulence with a range of different ionization fractions. Their results showed that at length scales comparable to the ambipolar diffusion scale, certain MHD waves, including Alfvén waves, are strongly damped, and this damping causes the velocity spectrum of ions to drop below that of the neutrals. They simulated ion and neutral spectral line profiles for ionization fractions ranging from 10⁻² to 10⁻⁶ and plotted their corresponding velocity dispersions against different beam sizes, producing a result very similar to Li & Houde (2008). Thus they were successful in numerically duplicating the observational results of Li & Houde (2008), including their determination of the turbulent ambipolar diffusion scale.

Assuming the gas to be well coupled to the field lines at large scales and that ion–neutral decoupling occurs at dissipation length scales, Li & Houde (2008) obtained an expression for B_pos by setting the effective magnetic Reynolds number R_m ∼ V_nL/β to 1 at the decoupling length scale. Here β = B²/4πn_iμν_i is the effective magnetic diffusivity with n_i the ion density, ν_i the collision rate of an ion with the neutrals, and μ the mean reduced mass characterizing such collisions. The aforementioned expression for B_pos is

$$\begin{eqnarray} B_{\mathrm{pos}} &=& \left(\frac{L^{\prime }}{0.5\:\mathrm{mpc}}\right)^{1/2}\left(\frac{V_{n}^{\prime }}{1\:\mathrm{km\,s^{-1}}}\right)^{1/2}\nonumber\\ &&\times \left(\frac{n_{n}}{10^{6}\:\mathrm{cm^{-3}}}\right)\left(\frac{\chi _{e}}{10^{-7}}\right)^{1/2}\:\mathrm{mG}, \end{eqnarray} \tag{ 1 }$$

where L' is the ion–neutral ambipolar diffusion decoupling length scale, V'_n is the neutral velocity dispersion at L', n_n is the neutral volume density, and χ_e is the fractional ionization. It should be noted that Equation (1) provides only an approximate value for B_pos as it assumes that ambipolar diffusion sets in exactly when R_m = 1, which is most likely too simple of an approximation. The decoupling scale and the velocity dispersion V'_n are determined by fitting the lower envelopes of the square of the ion and neutral velocity dispersion data (σ²_i and σ²_n, respectively) to Kolmogorov-type power laws. If the coexisting ions and neutral species are well coupled at large scales, their turbulent power spectra will have the same spectral index over the inertial range (i.e., at scales larger than the dissipation scale and smaller than the injection scale). Moreover, the cutoff range for the ion and neutral energy dissipation do not coincide, and assuming the ion energy spectrum to drop significantly at the ambipolar diffusion scale, the difference between σ²_i and σ²_n is the difference of the integral of the ion and neutral energy spectra between the ambipolar diffusion scale and neutral viscosity dissipation scale (this difference, a (see Equation (2) below), is the shaded area in Figure 3 of Li & Houde 2008). With all the above assumptions, the σ²_i and σ²_n data are thus fitted to the following expressions:

$$\begin{eqnarray} \sigma _{i}^{2}(L) & = & bL^{n}+a\\ \sigma _{n}^{2}(L) & = & bL^{n}.\nonumber \end{eqnarray} \tag{ 2 }$$

The observational determination of velocity dispersions from spectral line profiles is a measurement of a velocity field that is integrated along the line of sight within a column of gas of volume l² × L, where l is the beam size and L is the depth through the cloud. However, a given velocity field is associated with turbulent substructures that span volumes different than l² × L. It should therefore be examined whether the observational value of σ² at a given beam size l truly represents the actual velocity dispersion at that length scale. Falceta-Gonçalves et al. (2010) investigated this issue with spectral line profile synthesis and three-dimensional MHD numerical simulations. They showed that although the minimum observed velocity dispersions scale as the actual σ², the coincidence of individual values at a given scale depends on l and the sonic Mach number of the system. The observed dispersion is generally a good approximation of the actual dispersion, although the latter is slightly underestimated for supersonic turbulence and overestimated for a subsonic regime. They also showed that the power-law functions fitted to the lower envelopes and actual three-dimensional velocity dispersions from the simulations correspond well. The dispersions obtained using the lower envelope therefore appear to be reliable estimates for the actual velocity dispersions, as long as a large enough number of lines of sight are considered.

Li et al. (2010) studied the differences in the velocity dispersion for the spectral lines of H¹³CN (J = 1 → 0) and H¹³CO⁺ (J = 1 → 0) in DR21(OH) using the data of Lai et al. (2003b). They performed a similar analysis for HCN (J = 3 → 2) and HCO⁺ (J = 3 → 2) data which they obtained with the Submillimeter Array in NGC 2024, and also HCN (J = 4 → 3) and HCO⁺ (J = 4 → 3) lines previously observed by Houde et al. (2002) in M17 and used in Li & Houde (2008). In each source, the velocity dispersion comparisons were performed for two locations with measured B_pos line densities. They found that the ion–neutral velocity difference is almost constant in M17, which has a uniform density of magnetic field lines, while in NGC 2024 and DR21(OH) the ion–neutral velocity dispersion difference correlates well with the magnetic field line density. This behavior is predicted and consistent with the model of Li & Houde (2008). They also considered other possible parameters that could cause differences between the velocity dispersion of the two tracers, including optical depth, outflow enhancement, and the presence of hyperfine structures in the spectral lines, but only turbulent ambipolar diffusion could consistently explain their observations.

The spectra of an HCN and HCO⁺ pair, however, usually appear to be optically thick and are subject to saturation. The aim of this work is to test the turbulent energy dissipation model of Li & Houde (2008) with optically thin isotopologues, i.e., H¹³CN and H¹³CO⁺, to investigate whether the constant down shift of the velocity dispersion spectrum of the ion from that of the neutral is consistently reproducible with different species. Accordingly, the lower optical depths of the aforementioned molecules will allow us to potentially probe magnetic fields at greater depths in star-forming regions, where they may be stronger (Li & Houde 2008; Basu 2000; Crutcher 1999). Also, just like their main isotopologues they have very similar molecular masses and the critical densities associated with their J = 4 → 3 transition are approximately 10⁷ cm^-3 and 10⁶cm^-3, respectively. In particular, their line profiles will not be saturated or self-absorbed, which is likely to be the case with HCN and HCO⁺, and therefore the calculated line widths will be better representatives of the true velocity dispersions characterizing the medium under study.

We obtained position-switch grid maps of H¹³CN(J = 4 → 3) at 345.336 GHz and H¹³CO⁺(J = 4 → 3) at 346.998 GHz in weather grade 3 (0.08 < τ₂₂₅ < 0.12) in DR21(OH) during the months of 2008 August to November at the James Clerk Maxwell Telescope (JCMT) on the summit of Mauna Kea, Hawaii. The observations were performed using the HARP-B (Buckle et al. 2007) single sideband receiver and the ACSIS correlator configured for a bandwidth of 500 MHz and channel resolution of 61 kHz. HARP-B is a 4 × 4 element heterodyne focal plane array that uses superconductor–insulator–superconductor detectors. The 16 detectors have receiver temperatures of 94–165 K and are separated by 30'' with a foot print of 2'. To fully sample the 15'' JCMT beam, we mapped with sample offsets of 15'' around the reference position (α = 20^h39^m01^s and δ = 42°22'377, J2000) of the source. Since the difference between the frequencies of H¹³CN(J = 4 → 3) and H¹³CO⁺(J = 4 → 3) is 1.66 GHz, we made use of the special configuration capability of ACSIS to observe the two spectra simultaneously within the 1.8 GHz passband of the receiver. We used the previously published main beam efficiency of ≃62% determined with observations of Saturn at 345 GHz by Buckle et al. (2007).

Figure 1 shows the intensity maps of H¹³CN (J = 4 → 3) and H¹³CO⁺ (J = 4 → 3), integrated over a velocity range of −20 km s^-1 to 20 km s^-1 in DR21(OH). The high critical density of the J = 4 → 3 transition results in detection of radiation from high density regions only, as can be seen from the intensity contours congregating in an area of ≈40'' × 40'' around the reference position of the source. We performed the spectral line analysis with the Starlink⁴ and Gildas⁵ software packages. The spectra of neighboring regions were averaged across the maps to simulate beam sizes 15'',  30'',  45'', and 60'' at different lines of sight in the observed region. We then fitted multi-Gaussian profiles to each spectrum and the corresponding velocity dispersion σ was then calculated using these fits. The plot of the square of the velocity dispersions for which the dispersion is at least three times larger than its corresponding uncertainty is shown in the top panel of Figure 2 as a function of length scale, with the H¹³CN data plotted in black and H¹³CO⁺ in red. We did not include the error bars in this plot for a clearer representation of the data points. The H¹³CN data corresponding to the largest beam size are not shown, as they do not fit in the limits of the plot; we had low statistics at the 60'' scale and therefore were not able to produce quantitative data points for either of the species at this scale. The presence of data points at each length scale above the lower envelopes can be associated with the interception of the corresponding line of sight with several turbulent substructures of different velocities that brings an increase in the observed velocity dispersions. Therefore, the minimum values for σ² may be attributed to the line of sight that intercepts the fewest number of such turbulent sub-structures (Falceta-Gonçalves et al. 2010). All the minimum σ² values we obtained for both species arise from one location in DR21(OH), circled in red in Figure 1, thus further confirming the coexistence of H¹³CN and H¹³CO⁺ and showing the specific region along the line of sight that contains the least number of turbulent substructures. The spectra for the different beam sizes at this location are displayed in Figure 3, with the H¹³CO⁺ profiles plotted in red and H¹³CN in black and scaled to the temperature of H¹³CO⁺ for a clearer comparison of their line widths.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The lower envelopes of σ² for H¹³CN and H¹³CO⁺ are plotted with their corresponding uncertainties in the bottom panel of Figure 2, with Kolmogorov-type power-law fits given by Equation (2). As mentioned above, we did not have sufficient statistics to find the true lower envelope at the largest scale (i.e., 60'') for both spectra and thus used the σ² data at the three smallest beam sizes for the power-law fits. To enforce the same power-law fit for both species, we first fitted the difference between the square of the velocity dispersion data of H¹³CN and H¹³CO⁺ to a constant function to obtain a, and then fitted the sum of the two data sets to a power law of the form 2bLⁿ + a to obtain b and n. The fit results provide a = −0.45 ±  0.01 km² s^-2, b = 0.49 ± 0.03 $\mathrm{km^{2}\: s^{-2}}\:\mathrm{arcsecond^{-\mathit {n}}}$, and n = 0.36 ± 0.02. Although our fit is not as equally good for the two species at all length scales as it was for Li & Houde (2008), it is still very satisfactory and perfectly adequate for our study. Assuming a significant drop in the ion velocity spectrum at the dissipation scale, the velocity dispersion of ions will be approximately equal to V'_n at L' i.e., V'_n(L') ≃ σ²_i(L'). We can then write (Li & Houde 2008)

$$\begin{eqnarray} V_{n}^{{\prime }2}(L^{\prime }) & \simeq & 0.37bnL^{{\prime }n}\\ & \simeq & a+bL^{\prime }\nonumber \end{eqnarray} \tag{ 3 }$$

and therefore

$$\begin{equation} L^{{\prime }n}=-a/[b(1-0.37n)]. \end{equation} \tag{ 4 }$$

Taking the distance from DR21(OH) to be 3 kpc (Genzel & Downes 1977), we obtain L' ≃ 12 or ≃17 mpc and V'_n ≃ 0.26 km s^-1.

We used the RADEX software package (Van der Tak et al. 2007) to calculate the gas volume density and optical depths associated with the J = 4 → 3 transition of H¹³CN and H¹³CO⁺. Taking the values for the line widths and line temperatures of the spectral lines observed at the reference position of the source, and the column densities and excitation temperatures of the two species (≃10¹³cm^-2 and ≃13 K, respectively) from Hezareh et al. (2008), we get n_n ≃ 10⁶ cm^-3 and assuming the observed core to be as deep as it is wide, i.e., 40'', we obtain N(H₂) ≃ 2 × 10²⁴ cm^-3. Additionally, RADEX computes optical depths of 0.17 and 0.21 for H¹³CN and H¹³CO⁺ lines, respectively. This optical depth comparison agrees with the discussion by Li et al. (2010) mentioned earlier. We also use χ_e = 3.2 × 10⁻⁸ for the ionization fraction as obtained by Hezareh et al. (2008) for DR21(OH). Applying the parameters above in Equation (1), we determine a value of B_pos ≃ 1.7 mG. This value can be inaccurate by a factor of a few due to the uncertainties in n_n and χ_e, and the fact that Equation (1) is derived with the aforementioned assumption of determining the turbulent ambipolar diffusion length scale by setting R_m = 1.

Our results are in general consistent with the calculations of Li & Houde (2008) although the measurements were done on different sources and with different pairs of molecular species. An exception is perhaps the measured turbulent ambipolar dissipation scale for M17 reported by Li & Houde (2008) was 1.8 mpc, almost an order of magnitude smaller than what we measure for DR21(OH). This, however, may simply be the result of a stronger magnetic field for our source. Indeed, using the same fractional ionization of 3.2 × 10⁻⁸ for M17 yields B_pos ≃ 0.57 mG for that source, approximately three times weaker than for DR21(OH). A relatively weaker field strength for M17 than for DR21(OH) is also corroborated with existing CN (N = 1 → 0) Zeeman measurements (Crutcher et al. 1999; Falgarone et al. 2008). Furthermore, our result for B_pos is consistent with that of Lai et al. (2003a), who obtained a polarization map of dust emission in DR21(OH) and used the Chandrasekhar & Fermi (1953) method to calculate B_pos ≃ 1 mG.

We can now combine our value for B_pos with the Zeeman measurement of B_los ≃ 0.4 mG by Crutcher et al. (1999) in DR21(OH) and obtain a total strength of B ≃ 1.8 mG at a volume density n(H₂) ≃ 10⁶ cm^-3. We can also use this result to calculate a mass to flux ratio for the observed source. Taking M/Φ_B = 1.0 × 10⁻²⁰N(H₂)/|B|cm² μG (Mouschovias & Spitzer 1976) with N(H₂) ≃ 2 × 10²⁴cm^-2, we obtain M/Φ_B ≃ 10, which implies that DR21(OH) exhibits a highly supercritical core. Although this number could be overestimated by the column density considering our assumption for the core geometry, it is consistent with the fact that DR21(OH) has indeed revealed signs of high-mass star formation. Detection of OH masers (Norris et al. 1982), H₂O masers (Genzel & Downes 1977), and high-velocity outflows from its two dense cores, MM1 and MM2 (Lai et al. 2003a), are indications of ongoing star formation in this source.

We simultaneously observed and mapped H¹³CN (J = 4 → 3) and H¹³CO⁺ (J = 4 → 3) in DR21(OH) to test the turbulent energy dissipation model of Li & Houde (2008). The down shift of the ion velocity dispersion spectrum from that of the neutral is readily explained by the presence of ambipolar diffusion that causes a steeper turbulent power spectrum for the ions compared with that of the neutrals at dissipation length scales (Li & Houde 2008). We used our observational data to calculate a turbulent ambipolar diffusion scale of L' ≃ 17 mpc and a plane of the sky magnetic field strength of B_pos ≃ 1.7 mG for this source.

Our analysis shows that the difference between the velocity dispersion of the optically thin H¹³CN and H¹³CO⁺ spectral lines at different length scales exhibits a trend similar to the results originally obtained by Li & Houde (2008) with the optically thick main isotopologues. It therefore follows that the ion line narrowing effect first discussed by Houde et al. (2000a) and the ensuing conclusions of Li & Houde (2008) cannot be due to a relative optical depth effect between the ion and neutral species. This further corroborates prior results from Houde et al. (2000b), which first investigated this possibility, and more recently those of Li et al. (2010).

Although, as previously mentioned that optically thin spectral lines are preferable for this analysis, it should be noted that H¹³CN and H¹³CO⁺ are less abundant than the main species and are more difficult to map. Indeed, we observed the aforementioned pair in three additional sources, DR21(Main), NGC2264, and NGC2068, but DR21(OH) was the only source in which we obtained enough detections away from the core to perform our analysis. Moreover, observing HCN and HCO⁺ in regions where they are optically thick in the line core should not be too detrimental to the line width calculations, as the wings of such spectral lines are still optically thin and suitable for this study. More precisely, since the wings are more heavily weighted in line width calculations, the (multi-) Gaussian profiles fitted to such sections, away from the line peaks can be used for obtaining the velocity dispersions in both species (Li et al. 2010).

Finally, this technique has the advantage that it is not restricted to nearby sources, since the data analysis mainly involves line width calculations and does not require a spatial resolution high enough to resolve the ambipolar diffusion scale.

The authors thank the referee for providing helpful comments, and also S. Basu and W. Dapp for thoughtful discussions. M.H.'s research is funded through the NSERC Discovery Grant, Canada Research Chair, Canada Foundation for Innovation, Ontario Innovation Trust, and Western's Academic Development Fund programs.