OPACITY BROADENING OF ¹³CO LINEWIDTHS AND ITS EFFECT ON THE VARIANCE–SONIC MACH NUMBER RELATION

Supersonic magnetized turbulence is observed in multiple tracers across several different interstellar medium phases. This includes the neutral medium, traced by H i (see Chepurnov et al. 2010; Burkhart et al. 2010; Peek et al. 2011); the warm ionized medium, traced by Hα (see Hill et al. 2008); electron density fluctuations (Armstrong et al. 1995; Chepurnov & Lazarian 2010); synchrotron polarization (Gaensler et al. 2011; Burkhart et al. 2012); and the molecular medium (see Heyer et al. 1998; Goldsmith et al. 2008), which includes a variety of molecular traces including the often used carbon monoxide (CO) line.

Turbulence and magnetic fields in these environments are interrelated with a variety of physical processes, including cosmic ray transport (Yan & Lazarian 2004; Beresnyak et al. 2011), magnetic reconnection (Lazarian & Vishniac 1999; Kowal et al. 2009), and star formation (McKee & Ostriker 2007 and references therein). Furthermore, it is clear that, in order to understand MHD turbulence and related physical mechanisms, one needs to be able to measure basic plasma parameters such as the sonic and Alfvénic Mach numbers (M_s ≡ V_turb/c_s, where c_s is the sound speed, and M_A ≡ V_turb/V_Alfven, respectively).

It is clear from simulations, observations, and theoretical works that compressible turbulence, which generates shocks, is important for creating filaments and local regions of high density contrast (Kowal & Lazarian 2007). Shocks broaden the gas/dust density and column density probability distribution function (PDF) as well as increase the peaks and drag out the distribution's tails toward higher gas densities (Burkhart et al. 2009). Based on this observation, several authors (e.g., Vazquez-Semadeni 1994; Padoan et al. 1997; Federrath et al. 2010; Burkhart & Lazarian 2012, to name a few) have developed relationships between the PDF moments, such as the variance or standard deviation of density (σ_ρ/〈ρ〉), and the sonic Mach number, for example,

$$\begin{equation} \sigma ^2_{\rho /\langle \rho \rangle }=b^2M_s^2, \end{equation} \tag{ 1 }$$

where b is a parameter that depends on the type of turbulence forcing with b = 1/3 for pure solenoid forcing and b = 1 for pure compressive forcing (Federrath et al. 2008). Once gas becomes dense enough, collapse can occur via self-gravity and the PDF forms power-law tails toward higher density regions (Klessen et al. 2000; Collins et al. 2012). Additional variants of Equation (1) have also been developed, e.g., including plasma β (see Molina et al. 2012).

Equation (1) is hard to constrain observationally as volume densities and the three-dimensional (3D) velocity structure are not available from observations. Important information on turbulent supersonic motions in molecular clouds comes from the observed non-thermal broadening of the linewidths of different spectral lines, e.g., from carbon monoxide (CO) emission. However, CO, in particular ¹³CO, is often partially or fully optically thick and traces only a limited dynamic range of column densities (see Goodman et al. 2009; Burkhart et al. 2013a, 2013b).

Fortunately, dust extinction column density maps of infrared dark clouds (IRDCs), including mid infrared (MIR) and near infrared (NIR) wavelengths, can be used to trace a much larger dynamic range of densities in order to probe the PDF (A_v = 1–25 for NIR and A_v = 10–100 for MIR, see Lombardi & Alves 2001; Kainulainen et al. 2011; Kainulainen & Tan 2013, henceforth KT13). However, extinction maps carry no dynamical information regarding velocities and for this, molecular line profiles are needed to measure the dynamics of clouds, including the sonic Mach number for a given cloud temperature.

In this Letter, we investigate the robustness of measuring the sonic Mach number from maps of ¹³CO using synthetic observations derived from 3D MHD turbulence simulations. We further create synthetic dust maps in order to compare the PDF standard deviation with the measured sonic Mach number, following the approach in the observational work of KT13.

We generate a database of 12 3D numerical simulations of isothermal compressible (MHD) turbulence with resolution 512³ (all models are listed in Table 1). For models with M_s < 20 we use the MHD code detailed in Cho & Lazarian (2003) with large scale solenoidally driven turbulence. Our two models with M_s > 20 are from the AMUN code (Kowal et al. 2011). The magnetic field consists of the uniform background field and a time-dependent fluctuating field: B = B_ext + b(t); b(0) = 0. For more details on the simulations scheme see Burkhart & Lazarian (2012, henceforth BL12) and Burkhart et al. (2013b).

Table 1. Description of the 12 Simulations

	Super-Alfvénic (MAa ≃ 7.0)			Sub-Alfvénic (MAa ≃ 0.7)
Density	$\sigma ^{{\rm 1D}}_{^{13}{\rm CO}}$	Msb	τ	$\sigma ^{{\rm 1D}}_{^{13}{\rm CO}}$	Msb	τ
	(km s−1)			(km s−1)
	Msc ≃ 26.2			Msc ≃ 25.3
9n	6.75	34.4	0.0017	4.21	21.5	0.0022
275n	6.29	32.1	0.03	3.85	19.6	0.092
8250n	6.25	31.9	1.2	4.77	24.4	2.5
82, 500n	6.42	32.8	2.0	5.64	28.8	29.6
	Msc ≃ 9.0			Msc ≃ 7.9
9n	1.83	9.3	0.0066	1.26	6.4	0.0003
275n	1.92	9.8	0.3	1.41	7.2	0.3
8250n	2.20	11.2	3.0	1.76	9.0	3.2
82, 500n	2.32	11.8	54.7	1.94	9.9	74.4
	Msc ≃ 7.1			Msc ≃ 6.8
9n	1.34	6.8	0.008	1.00	5.1	0.0037
275n	1.47	7.5	0.1	1.12	5.7	0.4
8250n	1.71	8.7	3.0	1.54	7.9	10.0
82, 500n	1.81	9.2	68.9	1.68	8.6	96.1
	Msc ≃ 4.3			Msc ≃ 4.5
9n	0.84	4.3	0.0067	0.68	3.5	0.015
275n	0.84	4.3	0.3	0.85	4.4	0.4
8250n	1.00	5.1	13.0	1.09	5.6	13.0
82, 500n	1.06	5.4	56.6	1.18	6.0	128.0
	Msc ≃ 3.1			Msc ≃ 3.2
9n	0.49	2.5	0.016	0.57	2.9	0.0064
275n	0.57	2.9	0.5	0.61	3.1	0.1
8250n	0.73	3.7	20.0	0.80	4.1	14.0
82, 500n	0.78	4.0	140.0	0.86	4.4	20.3
	Msc ≃ 0.7			Msc ≃ 0.7
9n	0.16	0.8	0.082	0.16	0.8	0.094
275n	0.18	0.9	3.8	0.17	0.9	4.0
8250n	0.24	1.2	71.0	0.21	1.1	81.0
82, 500n	0.27	1.4	480.0	0.23	1.2	806.0

Notes. Each simulation has four different density values (given in Column 1) in order to probe varying optical depth. Columns 2–4 give cloud parameters for super-Alfvénic simulations while Columns 5–7 give parameters for sub-Alfvénic simulations. Those parameters are the one-dimensional (1D) velocity dispersion, estimated sonic Mach number (see Section 2) and the average optical depth of the cloud given by SimLine3D (see Ossenkopf 2002). The actual M_s as measured from the original MHD simulations with no radiative transfer is shown above each group of density scaling factors. ^aInitial Alfvén Mach Number. ^bEstimated M_s from the measured linewidths. ^cActual M_s from original MHD simulations.

Download table as:  ASCIITypeset image

Our models are divided into two groups corresponding to sub-Alfvénic (B_ext = 1.0) and super-Alfvénic (B_ext = 0.1) turbulence. For each group we compute five models with different gas pressures falling into subsonic and supersonic regimes with M_s ranging from 0.7 to 26 (see Table 1).

We post-process the simulations to include radiative transfer effects from the ¹³CO J = 2–1 line. For more information on this code and the assumptions involved see Ossenkopf (2002) and Burkhart et al. (2013a, 2013b). Our original simulation set the sonic Mach number using the sound speed. We rescale our simulation's velocity field in such a way that each has the same value of temperature and retain the original sonic Mach number. When applying radiative transfer post-processing we vary the density scaling factor by increasing and decreasing it by a factor of 30 from the standard numerical density of 275 cm⁻³(cm⁻³ = n henceforth), which represents a typical value for a giant molecular cloud. We also include a very high density case with 82, 500n. Each line of sight (LoS) has a spectral emission of ¹³CO, with resolution of 0.5 kms⁻¹ and is taken perpendicular to the mean magnetic field.⁶ The cube size of our clouds is 5pc, the gas temperature is 10K, and the CO abundance [¹³CO/H] = 1.5 × 10⁻⁶.

2.1. Analysis of the Synthetic Observations

Once we generate the synthetic ¹³CO line profile maps we measure the dispersion of the velocity profile using a Gaussian fit. Using a gas temperature of T = 10 K we can calculate the sonic Mach number as

$$\begin{equation} M_s=\frac{\sigma ^{{\rm 3D}}_{{\rm CO}}}{c_{s}}=\frac{\sqrt{3}\sigma ^{{\rm 1D}}_{{\rm CO}}}{c_{s}}, \end{equation} \tag{ 2 }$$

where $\sigma ^{{\rm 1D}}_{{\rm CO}}$ is the velocity dispersion in one dimension and c_s is the sound speed:

$$\begin{equation} c_{s}=\sqrt{\frac{K_{b}T}{\mu {\rm m}_{{\rm H}}}}. \end{equation} \tag{ 3 }$$

In addition to the ¹³CO simulations, we create synthetic dust maps using the column density maps from our MHD simulations with the LoS taken perpendicular to the mean magnetic field. We scale our column density mean value of unity to 2 × 10²² cm⁻² and then take a simple scaling law (see Bohlin et al. 1978) from column density to extinction as

$$\begin{equation} N_{\rm H}=1.9\times 10^{21}{\rm \,cm}^{-2}\frac{A_v}{{\rm mag}}. \end{equation} \tag{ 4 }$$

We then disregarded the column densities with A_v < 7 mag, following the procedure of KT13. In the left panel Figure 1, we show an example plot of the ¹³CO line profiles for same supersonic super-Alfvénic simulation with different density values, and hence different optical depths. The opacity broadening effect slightly widens the observed profile. In the right panel of Figure 1, we show an example of the synthetic dust PDFs (with A_v cut off of 7) of simulations with subsonic (supersonic) Mach number. An increase in sonic Mach number broadens the dust column density PDF distribution.

Zoom In Zoom Out Reset image size

In the following sections, we calculate the observed sonic Mach numbers from ¹³CO line maps and compare these with the actual sonic Mach numbers as measured in the simulations. We also calculate the PDFs of the synthetic dust extinction maps and examine the standard deviation–sonic Mach number relationship as if it were calculated from the observations.

Figure 2 plots M_s as measured from the ¹³CO line profiles (measured M_s) versus actual M_s of the simulations. We plot four different values of density, which effectively change the optical depth (see Table 1). We also overplot linear fits to the sub-Alfvénic (blue lines) and super-Alfvénic (black lines) simulations separately with the y-intercept held to pass through the origin. We fit slopes to each optical depth regime and list the values in Table 2.

Zoom In Zoom Out Reset image size

We find that the optically thin cases (triangle and diamond symbols) reproduce well the actual sonic Mach number from the measured one. However, for our highest optical depth cases the ratio between the observationally measured M_s and the actual M_s increases to between ≈1.16–1.3 (fitted slopes are reported in Table 2). Additional line broadening due to opacity enlarges the discrepancy between the measured and actual M_s as the optical depth increases. This occurs because the center of line emission becomes saturated and photons escape the cloud at the wings of the emission line. Indeed, CO opacity broadening can be explained by the curve of growth relationship as well as by numerous past studies (Lepine 1991; Leung & Brown 1977; Leung & Liszt 1976, see Table 2).

Additionally, we observe a magnetic field dependency: super-Alfvénic turbulence shows increased line broadening compared with sub-Alfvénic turbulence for all values of optical depth for supersonic turbulence. This effect may be due to the velocity dispersion of super-Alfvénic turbulence having a greater dependence on density fluctuations for supersonic turbulence (see Burkhart et al. 2009, Figures 9 and 10), which causes an additional broadening. This occurs regardless of the optical depth of the line (see Table 2). We investigated the effect of opacity broadening on a LoS parallel to the mean magnetic field and found results similar to those reported in Figure 2: the slopes of the measured versus actual M_s covering the range of optical depth values were 1.06 for sub-Alfvénic simulations and 1.28 for super-Alfvénic simulations. This indicates that the anisotropy present in MHD turbulence is not responsible for the line broadening and instead the strength of the magnetic field is more important.

We utilize both our synthetic dust maps and ¹³CO line dispersion data to estimate the standard deviation–sonic Mach number relation as it would be calculated from the observations and compare this with data from KT13. We calculate the standard deviation using two different methods. One is direct calculation of the column density standard deviation (σ_direct) divided by the mean value as

$$\begin{equation} \sigma _{N/\langle N\rangle }=\sqrt{\frac{1}{n}\sum ^n_{i=1}\left(\frac{N_i}{\langle N\rangle }-\left<\frac{N_i}{\langle N\rangle }\right>\right)^2}. \end{equation} \tag{ 5 }$$

The other method assumes a log-normal distribution and calculates the variance ($\sigma ^2_{{\rm ln}N}$) and mean (μ) values of the logarithm of the extinction maps by fitting a Gaussian distribution, illustrated in Figure 1, right panel. This method does not assume that we know a priori the actual mean value of the extinction distribution, which is indeed the case in the observations; however, the assumption of log normality is not always appropriate for all molecular clouds as gravity creates deviations from the log-normal distribution (see Collins et al. 2012). From the fit we are able to determine the parameters μ and $\sigma ^2_{{\rm ln}N}$, which are then used to calculate the mean (〈N〉) and standard deviation (σ_N) of the corresponding log-normal distribution as

$$\begin{equation} \langle N\rangle =e^{\mu +\sigma ^2_{{\rm ln}N}/2} \end{equation} \tag{ 6 }$$

and

$$\begin{equation} \sigma _N=\left[\left(e^{\sigma ^2_{{\rm ln}N}}-1\right)e^{2\mu +\sigma ^2_{{\rm ln}N}}\right]^{1/2}. \end{equation} \tag{ 7 }$$

For both methods, we calculate the standard deviation for an extinction distribution with an A_v cut off of ≈7 for compatibility with the method used in KT13 (hence we denote these as σ_{dir, cut} and σ_{lnN, cut}). We also investigate the standard deviation calculations for the complete A_v distribution. This is not available from observations; however, it is interesting to see how large of a difference in the standard deviation is observed between this idealized case and a more realistic distribution with a low A_v cut off.

We find that there is not a substantial difference between the values of the standard deviation calculated either with Equation (5) or Equations (6) and (7). This is expected for our data since we are dealing with log-normal distributions in the case of pure isothermal MHD turbulence. The differences between σ_lnN and σ_direct extend up to values of 0.15. We additionally find that there is not a large difference between the standard deviation of the full column density distribution and the column density distribution which employs a cut off value of A_v = 7. The differences between σ_lnN and σ_direct with a cut off value of A_v = 7 extend up to values of 0.13. However, for a subset of real clouds, it may be difficult to discern if one is dealing with a true log-normal distribution or the beginnings of a high density power-law tail and thus the difference in using direct calculation versus a log-normal fit to calculate the standard deviation of the column density in observational clouds might be higher.

Figure 3 shows the column density dispersion for four different methods of calculating the standard deviation: σ_direct, σ_{dir, cut}, σ_lnN, and σ_{lnN, cut} (direct calculation, direct calculation with an A_v = 7 cut off, log-normal fit, and log-normal fit with an A_v = 7 cut off, respectively) versus estimated M_s for sub-Alfvénic (top panel) and super-Alfvénic cases (bottom panel). The values of M_s plotted are the average M_s calculated from the linewidths of the four density cases with error bars representing the standard deviation between the values. We perform a linear fit to the (measured M_s, σ_N/〈N〉) data, as used in KT13,

$$\begin{equation} \sigma _{N/\langle N\rangle }=a_{1}\times M_{s}+a_{2}, \end{equation} \tag{ 8 }$$

where a₁ is the slope between the sonic Mach number and column density dispersion and a₂ is the intercept.

Zoom In Zoom Out Reset image size

Along with the synthetic observations presented in this Letter, Figure 3 overplots data from the IRDCs presented in KT13 (their Figure 7, left frame, M_s from a Gaussian fit, and reported in KT13's Table 1) which are shown with red diamonds. The KT13 values of σ_N/〈N〉 plotted here are calculated by direct calculation, and thus are most compatible with our method of direct calculation with an A_v = 7 cut off (blue triangles). Additionally, we include the predicted relation for the column density standard deviation–sonic Mach number relation from BL12, Equation (4) with b = 1/3 for solenoidal mixing and A = 0.11:

$$\begin{equation} \sigma _{N/\langle N\rangle }=\sqrt{\left(b^2M_{s}^2+1\right)^A-1}. \end{equation} \tag{ 9 }$$

The observational points of KT13 consistently show higher values of σ than the MHD simulations, regardless of the method of calculation for the standard deviation or the Alfvén Mach number. Upon further inspection, the KT13 points also have a spread in the values of σ for a given sonic Mach number, suggesting that other physics might come into play in the interpretation of the variance of the data. Processes that could create larger values of variance than expected from solenoidally driven turbulence include gravitational contraction (Collins et al. 2012) and compressive forcing (Federrath et al. 2010; Kainulainen et al. 2013).

Comparing the mean value of the slopes for the cases of direct calculation of the standard deviation with the A_v cut off, a₁ = 0.0103 ± 0.0047 (0.0136 ± 0.0050 for sub-Alfvénic; 0.0069 ± 0.0044 for super-Alfvénic), with KT13's Figure 7, left frame, which reports a₁ = 0.0095 ± 0.0066, we see that there is a good agreement. We see significantly less agreement of our a₁ values with the values of standard deviation derived from the fitted log-normals, a₁ = 0.0110 ± 0.0067 (0.0133 ± 0.0089 for sub-Alfvénic; 0.0086 ± 0.0046 for super-Alfvénic), compared to those of KT13 (their Figure 8, left frame), a₁ = 0.051 ± 0.018. The discrepancy may be due to PDF broadening by gravity in the observations (i.e., an unresolved power-law tail), thus distorting the standard deviation's relationship with the sonic Mach number.

The other possible explanation for a steeper relation could be that the driving in the KT13 data is compressive forcing; however, they derived an expression for b,

$$\begin{equation} b=\frac{\sigma _{\rho /\langle \rho \rangle }}{M_s}=\frac{\sigma _{N/\langle N\rangle }}{M_s}R^{-1/2}=a_1\times R^{-1/2}, \end{equation} \tag{ 10 }$$

where R is the 3D–2D variance ratio and was found to be between R = [0.03, 0.15] (Brunt 2010). KT13 estimated a b value with a 3σ uncertainty of $b=0.20^{+0.37}_{-0.22}$, which indicates that the driving is generally solenoidal to mixed (b = 1/3 for solenoidal and b = 1 for compressive driving; Federrath et al. 2008).

However, it is clear that applying a sonic Mach number opacity correction in the case of high optical depths will increase the measured value of b, since the slope a₁ will steepen by the correction factor. In the case of high optical depth this factor is as high as ≈1.3. Thus the slope, a₁, steepens by this factor and, assuming the observations are optically thick, the value of b calculated by KT13 will increase to $b=0.25^{+0.25}_{-0.15}$. These values still indicate solenoidal to mixed driving. Hence the most likely explanation for the higher variance values in KT13 as compared with simulations is the influence of gravity. This is not unexpected, as Collins et al. (2012, Table 1) shows that the variance measured from the PDFs in 3D density increases as gravity acts on the cloud.

The models of star formation must invoke stirring by turbulence, including collecting matter by compressible turbulent motions (see McKee & Ostriker 2007 and references therein) and the removal of magnetic flux from the collapsing region by the process of reconnection diffusion, which is much faster than the traditionally considered process of ambipolar diffusion (Lazarian et al. 2012). This motivates observational quantitative studies of turbulence and this Letter is a part of such studies.

This Letter studies the effect of self-absorption on the observational measurement of the sonic Mach number as measured from ¹³CO linewidths with a range of densities and optical depths. ¹³CO line broadening and extinction PDFs employed in this Letter to study the sonic Mach number are in no way the only methods available to researchers to study turbulence in molecular clouds. Sonic⁷ Mach numbers can be also successfully obtained with the kurtosis and skewness of the PDF distribution (Kowal et al. 2007; Burkhart et al. 2009, 2010). The power spectra of density and velocity can be obtained with velocity channel analysis and velocity coordinate spectrum techniques (Lazarian & Pogosyan 1999, 2004, 2006) which were successfully tested numerically (see Chepurnov et al. 2008) and applied to H i and CO data sets (see Padoan et al. 2009; Chepurnov et al. 2010; Lazarian 2009, for a review). In view of the complexity of astrophysical turbulence, simultaneous use of the combination of these techniques presents an unquestionable advantage.

Our empirical finding is that the evaluation of the Mach number from the linewidth is possible even for strongly self-absorbing species, but a correction factor should be applied for opacity broadening. For the range of absorption depths that we studied we found that this correction factor is at most 1.3. It also corresponds to the earlier one-dimensional studies of spectral line broadening in the presence of absorption in, e.g., Leung & Brown (1977) and Leung & Liszt (1976). An additional consideration is that IRDCs are cold (T = 10–40 K) and are likely to be close to isothermal, but if instead of 10 K the cloud had a different temperature there would be a factor of $1/\sqrt{T}$ change in M_s.

We create synthetic ¹³CO emission maps, with varying optical depth, and dust column density maps from a set of 3D MHD simulations. We derive a standard deviation–sonic Mach number relation, as found from the observations, and compare this with recent results from real clouds discussed in KT13. We find the following.

• 1.  
  Calculations of M_s from linewidths of ¹³CO are robust for optically thin ¹³CO but are overestimated by a factor of up to ∼1.3 for optically thick clouds. This is due to the well-known, but often overlooked, effect of opacity broadening.
• 2.  
  This overestimation of the sonic Mach number as derived from ¹³CO linewidths will cause the slope of the σ_N/〈N〉–M_s relation to become shallower and this will result in lower values of the measured b parameter.
• 3.  
  The σ_N/〈N〉values of clouds reported in KT13 are larger than values found from ideal solenoidally driven simulations of turbulence. This could be due to the fact that in real molecular clouds gravitational influences exist, which will increase the measured column density standard deviation.

We thank the anonymous referee for helpful comments. This work is supported by continuous grants from INEspaço/FAPERN/CNPq/MCT. C.C. acknowledges a graduate PDSE/CAPES grant Process n^o9392/13–0. B.B. acknowledges support from the Wisconsin Space grant. A.L. acknowledges the Vilas Associate Award and the NSF grant AST-1212096. A.L. and B.B. acknowledge the Center for Magnetic Self-Organization in Laboratory and Astrophysical Plasmas and the IIP/UFRN (Natal) for hospitality. V.O. and J.S. acknowledge support from the Deutsche Forschungsgemeinschaft (DFG) project n^o0s ∼ 177/2 − 1. V.O., J.S., and J.K. were supported by the central funds of the DFG-priority program 1573 (ISM-SPP).