The Multi-phase Turbulence Density Power Spectra in the Perseus Molecular Cloud

TL;DR

This study analyzes the spatial power spectra of various interstellar medium tracers in the Perseus molecular cloud, revealing differences in turbulence characteristics and the impact of opacity and chemistry on observational tracers.

Contribution

It provides the first detailed comparison of turbulence power spectra across multiple tracers in Perseus, highlighting the effects of opacity and chemistry on turbulence measurements.

Findings

HI traces a non-gravitating, transonic medium

Dust extinction indicates a self-gravitating, supersonic medium

CO lines are affected by opacity, complicating turbulence analysis

Abstract

We derive two-dimensional spatial power spectra of four distinct interstellar medium tracers, HI, $^{12}$CO($J$=1--0), $^{13}$CO($J$=1--0), and dust, in the Perseus molecular cloud, covering linear scales ranging from $\sim$0.1 pc to $\sim$90 pc. Among the four tracers, we find the steepest slopes of $-3.23\pm0.05$ and $-3.22\pm0.05$ for the uncorrected and opacity-corrected HI column density images. This result suggests that the HI in and around Perseus traces a non-gravitating, transonic medium on average, with a negligible effect from opacity. On the other hand, we measure the shallowest slope of $-2.72\pm0.12$ for the 2MASS dust extinction data and interpret this as the signature of a self-gravitating, supersonic medium. Possible variations in the dust-to-gas ratio likely do not change our conclusion. Finally, we derive slopes of $-3.08\pm0.08$ and $-2.88\pm0.07$ for the ${}^{12}$CO(1--0) and ${}^{13}$CO(1--0) integrated intensity images. Based on theoretical predictions for an optically thick medium, we interpret these slopes of roughly $-3$ as implying that both CO lines are susceptible to the opacity effect. While simple tests for the impact of CO formation and depletion indicate that the measured slopes of ${}^{12}$CO(1--0) and ${}^{13}$CO(1--0) are not likely affected by these chemical effects, our results generally suggest that chemically more complex and/or fully optically thick media may not be a reliable observational tracer for characterizing turbulence.

Figures

Tables

Table 1. Table 1 Characteristics of the data used in our study

	$N$(Hi)	Corrected $N$(Hi)	$A_V$	$I$(12CO(1–0))	$I$(13CO(1–0))
Image Sizea	${14.8}^{\circ }\times {9.0}^{\circ }$	${14.8}^{\circ }\times {9.0}^{\circ }$	${15.1}^{\circ }\times {15.5}^{\circ }$	${6.6}^{\circ }\times {3.5}^{\circ }$	${6.6}^{\circ }\times {3.5}^{\circ }$
	(77.5 pc $\times$ 47.1 pc)	(77.5 pc $\times$ 47.1 pc)	(79.1 pc $\times$ 81.2 pc)	(34.6 pc $\times$ 18.3 pc)	(34.6 pc $\times$ 18.3 pc)
Pixel Sizea	4.3′	4.3′	5.0′	46.0′′	46.0′′
	(0.4 pc)	(0.4 pc)	(0.4 pc)	(0.07 pc)	(0.07 pc)
Angular Resolution	4.3′	4.3′	5.0′	46.0′′	46.0′′
Median 1$\sigma$	6 $\times$ 1019 cm-2	8 $\times$ 1019 cm-2	0.2 mag	0.5 K km s-1	0.2 K km s-1

Table 2. Table 2 Summary of Moduli Properties

	Hi	$A_V$	12CO(1–0)	13CO(1–0)
Grid Size [pixels]	206$\times$126	200$\times$200	520$\times$273	520$\times$273
Angular Extent [deg]	14.8	16.7	6.6	6.6
Modulus Pixel Resolution [$\lambda$]	3.9	3.4	8.6	8.6
Longest Spatial Frequency [$\lambda$]	397.6	341.0	2241.0	2241.0
Smallest Effective Baseline [$\lambda$]	10.1	8.9	27.5	27.5
Maximum Linear Scale [pc]	77.7	87.9	34.8	34.8
Minimum Linear Scape [pc]	0.8	0.9	0.1	0.1
Maximum Effective Linear Scale [pc]	29.6	33.7	10.9	10.9

Table 3. Table 3 Power spectral slopes of turbulence for various environments

Environment	3D Density Spectruma	Reference
${ℳ}_s$ $\lesssim$ 1	$\approx k^{-11/3}$	Kolmogorov (1941); Goldreich & Sridhar (1995); Cho & Lazarian (2002, 2003)
${ℳ}_s$ $>$ 1	shallower than $k^{-11/3}$	Beresnyak et al. (2005); Kowal et al. (2007)
${ℳ}_s$ $>$ 1 and self-gravitating	shallower than pure compressible (positive values possible)	Fleck (1996); Federrath & Klessen (2013); Collins et al. (2012); Burkhart et al. (2015a)
optically thick	$\approx k^{-3}$	Lazarian & Pogosyan (2004); Burkhart et al. (2013);

Table 4. Table 4 Properties of the dark and star-forming regions a in Perseus

Property	B5	IC 348	B1	NGC 1333	Westend
b$\tau$(12CO)	1.0	0.5	0.5	0.3	0.9
c$\tau$(13CO)	0.3	0.4	0.3	0.3	0.4

Table 5. Table 5 Summary of Previous SPS Studies

Object	Tracer	Optical Depth	Scale	Slope	Reference
Polaris Flare	Dust	Thin	0.01$\sim$8 pc	$-2.7\pm 0.1$	(1)
High Latitude Cirrus	Dust	Thin	0.01$\sim$50 pc	$-2.9\pm 0.1$	(2)
CNM	Hi Absorption	Thin	0.01$\sim$3 pc	$-2.75\pm 0.25$	(3)
WNM	Hi Emission	Thin	–	$-4$	(4)
Ursa Major	Hi Emission	Thin	0.1$\sim$25 pc	$-3.6\pm 0.2$	(5)
MBM16	Hi Emission	Thin	0.1$\sim$20 pc	$-3.7\pm 0.2$	(6)
North Ecliptic Pole	Hi Emission	Thin	–	from $-2.86\pm 0.04$	(7)
				to $-2.59\pm 0.07$

Equations

$$P(\nu )=\mathfrak{F}\left(\nu \right)\times {\mathfrak{F}}^{\ast }\left(\nu \right)$$

$$L[pc]\sim \frac{d[pc]}{k[\lambda ]},$$

$$MAD=median\left(|X_i-median\left(X\right)|\right).$$

Full text

1. Introduction

Turbulence has been considered as one of the key physical processes for the evolution of the interstellar medium (ISM) (e.g., Elmegreen & Scalo 2004; McKee & Ostriker 2007; Lazarian 2009). Yet many questions still remain open, in particular with regard to the role of turbulence in the formation and evolution of molecular clouds (MCs) and subsequent star formation. Some of the important questions include: What are the injection sources of turbulent energy and on what scales? How is the transition from atomic (Hi) to molecular hydrogen (H2; Bialy et al. 2017), likely a prerequisite of star formation, influenced by turbulence? How does turbulence manifest itself as a function of ISM phase (e.g., cold and warm neutral medium (CNM and WNM) and molecular gas)? How does opacity affect an observer’s ability to accurately characterize the properties of turbulence?

As for the origin of turbulence, it is believed that turbulence is driven on a variety of spatial scales and cascades down to smaller scales, as evidenced by the the fractal structure of the ISM (e.g., Stutzki et al. 1998; Stanimirović et al. 1999; Elmegreen et al. 2001a). The accretion of the circumgalactic medium and galactic-scale gravitational instabilities likely trigger turbulence on large scales (e.g., Klessen & Hennebelle 2010; Krumholz & Burkhart 2016, while stellar feedback such as outflows and supernova explosions inject energy on smaller scales (e.g., Krumholz et al. 2006; Zamora-Avilés et al. 2012; Padoan et al. 2016).

Depending on the specific driver, the characteristics of turbulence will then be imprinted within the ISM mainly as three-dimensional density and velocity fluctuations, and these fluctuations have been traditionally studied via correlation functions such as the spatial power spectrum (SPS) (e.g., Crovisier & Dickey 1983), $\Delta$-variance (e.g. Stutzki et al. 1998, and structure function (e.g., Padoan et al. 2002; Burkhart et al. 2015b). In particular, the SPS approach has been applied to observations of various Galactic and extragalactic environments (e.g., Plume et al. 2000; Dickey et al. 2001; Elmegreen et al. 2001a; Burkhart et al. 2010; Combes et al. 2012; Zhang et al. 2012; Pingel et al. 2013), showing power spectral slopes $\beta$ roughly ranging from $-2.7$ to $-3.7$ depending on the used tracers (e.g., Hi, carbon monoxide (CO), and dust). These slopes essentially provide information on the relative amount of structures as a function of spatial scale and can be compared with theoretical models of turbulence (mainly numerical simulations) to characterize turbulence cascade (e.g., Burkhart et al. 2010), to determine the influence of shocks (e.g., Beresnyak et al. 2005), to reveal the injection and dissipation scales of turbulent energy (e.g., Kowal & Lazarian 2007; Federrath & Klessen 2013; Chen et al. 2015), and to trace the evolution of MCs (e.g., Burkhart et al. 2015a). The proximity and a wealth of multi-wavelength observations make MCs in the solar neighborhood an ideal laboratory for probing the impact of turbulence on the formation and evolution of MCs. In this paper, we focus on the Perseus molecular cloud, which is a nearby ($\sim$300 pc; e.g. Herbig & Jones 1983; Černis 1990), low-mass ($\sim$2 $\times$ 104 M*⊙*; e.g. Sancisi et al. 1974; Lada et al. 2010) cloud. Its star formation activities, as well as atomic and molecular gas contents, have been extensively examined over the past decade (e.g., Ridge et al. 2006; Jørgensen et al. 2007; Pineda et al. 2008; Lee et al. 2012, 2014, 2015; Mercimek et al. 2017), revealing that the cloud consists of several individual dark and star-forming regions (e.g., B5, B1, B1E, IC348, and NGC1333) and is actively forming low- to intermediate-mass stars (see Bally et al. 2008 for a review).

The properties of turbulence in Perseus has been studied in the past, mainly analyzing the observations of CO isotopologues (12CO($J$=1–0, 3–2) and 13CO($J$=1–0, 2–1)111In this paper, we quote the 12CO($J$=1–0) and 13CO($J$=1–0) lines as 12CO and 13CO, while specifying other CO transitions.; e.g., Padoan et al. 1999; Bensch et al. 2001; Sun et al. 2006) and dust (e.g., Schneider et al. 2011; Burkhart et al. 2015c) with the SPS, probability distribution function, and $\Delta$-variance methods.

In this paper, we extend these previous studies, with an aim of probing the characteristics of turbulence in the multi-phase ISM more comprehensively. To this end, we perform a SPS analysis of multi-wavelength data, including recent high resolution ($\sim$4.3*′* or $\sim$0.4 pc scales) Hi observations of Perseus (Lee et al. 2012; Stanimirović et al. 2014; Lee et al. 2015), and systematically compare the derived power spectra. In addition, we examine several biases that could affect SPS analyses, such as the impact of opacity and CO chemistry, and discuss the optimal tracer(s) for studying the turbulent environments of MCs.

Our paper is structured as follows. In Section 2, we describe the observations of each tracer used in this study. In Section 3, we then summarize the SPS technique and derivation of uncertainties. Section 4 presents the derived power spectra, and Section 5 discusses how our results relate to the theoretical predictions of the SPS, as well as previous observational studies. Finally, in Section 6, we summarize our conclusions.

2. Data

In this section, we describe the gas and dust data used in our analysis of the SPS. Several specifics of the data, e.g., image/pixel sizes, angular resolutions, and median 1$\sigma$ uncertainties, are summarized in Table 1.

2.1. HI

We use the Hi column density image of Perseus from Lee et al. (2012). To derive the Hi column density on 4.3*′* scales ($\sim$0.4 pc at the distance of Perseus), Lee et al. (2012) used Hi cubes from the GALFA-Hi survey222https://purcell.ssl.berkeley.edu/ (Peek et al., 2011) and integrated the Hi emission from $V_{\rm LSR}$ = $-$5 to $+$15 km s*-1* under the optically thin assumption. This velocity range was determined based on the correlation between the derived $N$(Hi) and the 2MASS $A_{V}$ from the COMPLETE survey333https://www.cfa.harvard.edu/COMPLETE/ (Ridge et al., 2006). The final Hi column density image is centered at (R.A., decl.) = (03h29m52s,$+$30*∘34′01′′) with a size of $\sim$14.7∘* $\times$ 9.0*∘*.

In addition, we make use of the Hi column density image from Lee et al. (2015), which was corrected for optical depth effects, with an aim of probing the impact of optically thick Hi on the slope of a power spectrum. In essence, Lee et al. (2015) derived the empirical correction factor for the optically thick Hi in Perseus based on Arecibo Hi absorption measurements toward 26 background continuum sources (Stanimirović et al., 2014) and applied the correction to the Hi column density image from Lee et al. (2012) on a pixel-by-pixel basis. The amount of the cold Hi was found not substantial (median cold-to-total Hi ratio of $\sim$0.3) and the opacity correction was hence small (only up to $\sim$1.2). The final opacity-corrected Hi column density image is presented in Figure 1, and we refer to Lee et al. (2012) for details on the derivation of the Hi column density images.

2.2. $A_{V}$

We use the $A_{V}$ image from the COMPLETE survey, which was produced by applying the NICER technique (Lombardi & Alves, 2001) to 2MASS data444http://www.ipac.caltech.edu/2mass/releases/allsky/doc/explsup.html at 5*′* resolution. For our analysis, we regrid the $A_{V}$ image to have independent pixels with a size of 5*′* (Figure 1). Note that Lee et al. (2012) also derived the $A_{V}$ image of Perseus based on 60 $\mu$m and 100 $\mu$m data from the IRIS survey (Miville-Deschênes & Lagache, 2005). We, however, do not use this image due to the presence of masked pixels, which will cause the “Gibbs ringing” in our power spectrum analysis (see Section 3 for details). In Lee et al. (2012), $\sim$20% of the total pixels were blanked to exclude the lines of sight where dust emission modeling was unreliable due to possible contaminations from foreground stars, a background Hii region, and the nearby Taurus molecular cloud. A pixel-by-pixel comparison shows that the Lee et al. (2012) $A_{V}$ and the 2MASS $A_{V}$ are in agreement: they have overall a 1:1 relation with a small median difference of $\sim$0.2 mag (which is comparable to the median 1$\sigma$ uncertainty of the 2MASS $A_{V}$ image; Table 1). This agreement primarily arises from the fact that Lee et al. (2012) estimated $A_{V}$ by calibrating the derived dust optical depth at 100 $\mu$m with the 2MASS $A_{V}$. The scatter around the 1:1 relation does exist, however, and this can be mostly attributed to the variations in the dust temperature along a line of sight (Goodman et al., 2009).

2.3. CO

Finally, we use the 12CO($J$=1–0) and 13CO($J$=1–0) cubes from the COMPLETE survey (Ridge et al., 2006). These cubes were obtained with the 14-m FCRAO telescope at 46*′′* resolution ($\sim$0.07 pc at the distance of Perseus), covering an area of $\sim$6.6*∘* $\times$ 3.5*∘* with a center of (R.A.,decl.) = (03h36m03s,$+$31*∘44′10′′). To derive the 12CO and 13CO integrated intensity images (Figure 1), we first regrid the cubes to have a pixel size of 46′′* and integrate the CO emission from $V_{\rm LSR}$ = $-$5 to $+$15 km s*-1*.

3. SPS Analysis

3.1. Derivation of the SPS

The SPS, i.e. the Fourier transform of the two-point auto-correlation function, provides information on the properties of turbulence, including the injection and dissipation scales, as well as inertial range scaling. The two-dimensional SPS is defined as:

[TABLE]

where $\nu$ is the wavenumber ($\nu=\frac{2\pi}{L}$ where $L=$ length scale) and $\mathfrak{F}\left(\nu\right)$ is the Fourier transform of the field under study, which for our purpose is the Hi/dust column density and 12CO/13CO integrated intensity images of Perseus.

Our derivation of the SPS is similar to previous studies (e.g., Crovisier & Dickey 1983; Green 1993; Stanimirović et al. 1999; Elmegreen et al. 2001b; Muller et al. 2004; Pingel et al. 2013).

In practice, we take the 2D Fourier transform of each column density or integrated intensity image and square the modulus of the transform — <$$\Re^{2}+\Im^{2}$$>, where $\Re$ and $\Im$ are the real and imaginary components of the 2D Fourier transform, respectively — in order to derive what we define as the modulus image with which we perform our statistical analysis.

We relate each baseline to a linear length scale by adopting the distance to Perseus to be 300 pc. Specifically, the linear length scale, is related to the spatial frequency ($k$; measured in units of wavelength $\lambda$) by

[TABLE]

where $d$ is the distance to Perseus and $k$ represents the spatial frequency. For the Hi data, the maximum angular scale of 14.8*∘* corresponds to a linear length scale of 77.5 pc, which sets the pixel resolution of 3.9$\lambda$. The origin is set at the center of the Fourier plane, and the spatial frequencies range from 0 to 397.6$\lambda$. The grid size, pixel resolution, and spatial frequency range for the Hi and other tracers are summarized in Tables 2.

To derive the SPS, we first place twelve annuli (uniformly spaced in log space) on each modulus image such that the innermost annulus corresponds to the largest length scale, while the outermost annulus corresponds to the smallest length scale. In this case, the linear length scale ranges from 77.7 pc down to 0.8 pc for the Hi, 87.9 pc down to 0.9 pc for the $A_{V}$, and 34.8 pc down to 0.1 pc for both CO data sets. As an example, we present the modulus image of the Hi with the annuli boundaries overlaid in Figure 2. Next, assuming azimuthal symmetry in the modulus image, we calculate the median value of the pixel distribution contained within each annulus and plot the median values as a function of decreasing linear scale. The derived SPS for each tracer is presented in Figure 3.

Our use of the median values is motivated by the significant number of bright pixels along the axes of each modulus image caused by the Gibbs phenomenon. The Gibbs phenomenon causes $sin(x)/x$ (where $x$ is the pixel coordinate) ringing along the axes of the modulus image due to the Fourier transform of discontinuous image edges. The resulting “cross” in the modulus image introduces a high-end tail in the pixel distributions within individual annuli. Traditionally, a power spectrum analysis has been done using mean values (e.g.,  Stanimirović et al. 1999; Dickey et al. 2001; Muller et al. 2004). We experimented with padding the images with zeros, as well as applying Gaussian tapers similar to the methods by Muller et al. (2004) to smoothly reduce the intensity to zero at the image edges, but found that the pixel distributions, especially in the Hi image, still possess a high-end tail of pixel values.

To demonstrate the impact of the high-end tail, we present the probability distribution function (PDF) of the pixel distribution within the largest annulus for each tracer in Figure 4. On average, we find that each pixel distribution can be characterized reasonably well by a Gaussian with some high-end tail. We thus fit a Gaussian function (red dashed line) to each PDF and overlay the arithmetic mean (black line), the fitted mean (red solid line), and the median (green line). For each tracer, the median of the PDF is consistently within $<$ 1% of the fitted mean, suggesting that the median is a good representative value of each pixel distribution. The arithmetic mean, however, is consistently higher than what is expected from a Gaussian distribution, mainly due to the high power values created by the Gibbs phenomenon. The Hi PDF clearly demonstrates the impact of these artificially high power values on the computed mean: only 2% of the pixels have power values greater than 3$\sigma$ as determined by the Gaussian fit, but the arithmetic mean of the distribution is an order of magnitude higher than the median and mean returned by the fit. On the other hand, the mean values of the pixel distributions for other tracers are not as susceptible to the effects of the Gibbs phenomenon as Hi. The relatively well-behaved distributions for these tracers are due to the gradual decrease of integrated emission towards the edge of their respective images. The general uniformity in the Hi integrated emission (the pixel values vary by only a factor of two over a large area of $\sim$15${}^{\circ}\times 10^{\circ}$), and to some extent the dust tracer, enhances the edge effects of the Fourier transform, which in turn increases the discrepancy between the median and the arithmetic mean. Considering how well the median characterizes a pixel distribution even in the presence of outlying values, using the median as opposed to the mean of the pixel distributions of the modulus image is justified to derive power spectral slopes. We examine the PDFs of the pixel distributions in all individual annuli for each tracer and indeed confirm that the median is a robust representation of the distributions.

Our final concern about the measured median power within each annulus is the contributions from noise and beams. We first find that the power spectra of Hi, 12CO(1–0), and 13CO(1–0) derived from emission-free channels only are flat (which is expected for white noise) and three or four orders of magnitude lower than those in Figure 3. We therefore conclude that the contribution from noise to our derived power spectra is negligible.

As for the effect of the beam, we expect no significant impact since the pixels in each of our tracer maps are independent with each other with the size equal to the respective beamsize. We confirm this by deriving the SPS of an idealized Gaussian beam with a FWHM set to the respective beamsize for each tracer and finding that the beam has a negligible contribution to the power at scales larger than the pixel resolution.

3.2. Derivation of SPS Uncertainties

Following Pingel et al. (2013), we initially estimated the errors of the median values by (1) creating 1000 simulated integrated intensity or column density images (each pixel value of a simulated map was determined by drawing from a Gaussian distribution centered on the original map pixel value and width characterized by the associated 1$\sigma$ uncertainty), (2) deriving the SPS for each image, and (3) taking the final uncertainty for each annulus as the standard deviation of the measured 1000 median power values. For each tracer, we then found that this Monte Carlo approach results in an order of magnitude lower uncertainty on the largest length scale as compared to smaller length scales.

This discrepancy in uncertainties between different scales mainly arises from a broad difference in the number of pixels available for each annulus, i.e., the small number of pixels for the largest length scale tends to cause a small spread in individual simulated values, which drives the small standard deviation between the 1000 realizations. The lower uncertainties on large length scales in turn leads to artificially shallow slopes, since the corresponding data points are weighted more in our linear regression fit.

We therefore require a method for error characterization that is robust against the difference in sample sizes, as well as high-end tail pixel values. To this end, we employ the median absolute deviation (MAD) method. The MAD of some distribution $X$, which depends on a univariate data set made up of individual samples ($X_{1}$,$X_{2}$,…$X_{n}$), is defined as

[TABLE]

In practice, we calculate the absolute deviations of the pixel values from the computed median for each annulus and use the median of these deviations as a final error. This MAD is more resilient against outlying values than standard deviation. For example, in the case of standard deviation, the large differences between pixel values and the artificially high mean due to the Gibbs ringing are squared causing these large discrepancies to be weighted more. In the MAD calculation, on the other hand, a small number of large absolute deviations becomes less relevant due to the use of the median.

4. Results

The derived power spectra in Figure 3 show that the four tracers have a range of slopes, suggesting that turbulence manifests itself in different ways depending on the medium: the 2MASS $A_{V}$ ($-2.72\pm 0.12$) shows the shallowest slope, while the 12CO(1–0) ($-3.08\pm 0.08$), 13CO(1–0) ($-2.88\pm 0.07$), and Hi ($-3.23\pm 0.05$ and $-3.22\pm 0.05$ for the uncorrected and corrected Hi respectively) show the intermediate and steepest slopes. Among these four tracers, 12CO(1–0), 13CO(1–0), and dust have been previously examined in the context of the turbulent environment of Perseus. One such study was Padoan et al. (2006b). In their SPS analysis, the same 13CO(1–0) image from Ridge et al. (2006) was used to derive the spectral slope ($\beta_{I}$) of $-1.99\pm 0.05$ by measuring the total power, rather than the average or median value, within each wavenumber shell. Since we take a similar derivation approach for the power spectral index as Lazarian & Pogosyan (2000), our slope ($\gamma$) is related to the estimate by Padoan et al. (2006b) in the following way: $\beta_{I}$ $=$ $1-\gamma$. In order to properly compare with our result, a value of one must then be subtracted from the Padoan et al. (2006b) slope. Considering other differences between the two studies (e.g., angular resolution of 46*′′* and 92*′′* for our study and Padoan et al. (2006b) respectively), we conclude that our 13CO(1–0) power spectral slope is consistent with Padoan et al. (2006b)

Sun et al. (2006) is another study which utilized similar 12CO(1–0), 13CO(1–0), and 2MASS $A_{V}$ datasets to investigate the properties of turbulence in Perseus. To measure power spectral indices, they employed the $\Delta$-variance technique and found the slopes of $-3.08\pm 0.04$, $-3.09\pm 0.09$, and $-2.55\pm 0.02$ for 12CO(1–0), 13CO(1–0), and dust respectively. Considering slight differences in the datasets, as well as the fitting ranges, we conclude that our results are in relatively good agreement with Sun et al. (2006).

Interestingly, for all tracers, we do not find a break over two orders of magnitude in length scales ($\sim$0.1 pc to $\sim$10 pc). This implies that turbulence is most likely driven on the scales larger than the size of the entire molecular cloud and dissipated down to the scales smaller than $\sim$0.1 pc (see Section 5.5 for more discussions). That said, that if turbulent energy is injected at multiple scales, the density slope is shown to be particularly influenced by large-scale drivers, which make it difficult to trace the injection scale turnover (Yoo & Cho, 2014). In addition, our power spectra provide critical information on the statistical properties of density fluctuations in three dimensions. The power spectrum of a column density or integrated emission image can be used to infer the power spectrum of the 3D density field when the tracer is in the “thick slice” limit. The “thick slice” limit is formally defined by Lazarian & Pogosyan (2000) as the case where the velocity dispersion of the tracer is smaller than the channel width over which emission is integrated. In this case, velocity fluctuations are averaged out, and the measured 2D power spectrum can be directly related to the 3D density power spectrum. For the Hi, 12CO(1–0), and 13CO(1–0) in and around Perseus, this is indeed the case: the measured velocity dispersions ($\sim$1–2 km s*-1* for 12CO(1–0) and 13CO(1–0) and $\sim$5–7 km s*-1* for Hi) are smaller than the velocity range ($\delta$v = 20 km s*-1*) over which each emission is summed. The same argument cannot be made for the 2MASS $A_{V}$ though, since it lacks velocity information. However, considering that the dust column density traces the total gas along a line of sight (atomic + molecular gas), we expect that the 2MASS $A_{V}$ is also most likely in the “thick slice” limit.

Whether or not velocity fluctuations are sufficiently averaged out can be more thoroughly examined by the Velocity Channel Analysis (VCA) developed by Lazarian & Pogosyan (2000). For this analysis, the velocity resolution of a data cube is varied from its instrumental resolution to coarser and coarser channel widths until the final integrated intensity image is created. At each velocity resolution, 2D power spectra are calculated and averaged together to estimate a representative slope. The velocity slices are then considered “thick” if the average slope saturates over several values of the velocity thickness. As a preliminary analysis, we applied the VCA to the Hi, 12CO(1–0), and 13CO(1–0) cubes of Perseus, finding that the power spectrum slope saturates well before the final channel width of 20 km s*-1*. This suggests that we are indeed in the “thick slice” limit, probing density fluctuations in our 2D power spectrum analysis. While the VCA enables us to probe 3D velocity fluctuations which are potentially important for disentangling various drivers of turbulence (e.g., Offner & Arce 2015), it is out of the scope of this paper, since we mainly focus on comparing the statistical properties of density fluctuations revealed by several tracers. We plan to discuss the impact of velocity fluctuations through the VCA method in a forthcoming paper.

5. Discussion

In this section, we discuss the results of the individual tracers and compare them with previous studies. Several implications of our results will be also discussed.

5.1. Theoretical Predictions for the SPS of Turbulence

We begin our discussion by first summarizing theoretical expectations for the SPS of the multi-phase ISM (Table 3). In the presence of supersonic turbulence without self-gravity, which is likely the case for the CNM in and around molecular clouds, the density spectral slope is expected to be shallower than the Kolmogorov index of $-11/3$ due to shock compression in three dimensions. Shocks can create small-scale density enhancements (e.g., Beresnyak et al. 2005; Kowal & Lazarian 2007; Burkhart et al. 2010), which in turn induce more power on small scales as compared to incompressible turbulence. In a weakly magnetized incompressible medium or a medium with no shock, the spectrum is very close to the Kolmogorov index of $-11/3$. Increasing the strength of the magnetic field in this limit increases the magnetic pressure and can steepen the power spectrum to $-13/3$ (e.g., Kowal et al. 2007).

If self-gravity plays a role in shaping density structures, which is expected to be the case for the dense molecular medium, the 3D density and column density spectral slopes can become increasingly shallow as compared to non-gravitating supersonic turbulence (e.g., Ossenkopf et al. 2001; Collins et al. 2012; Federrath & Klessen 2013; Burkhart et al. 2015a). This is due to gravity enhancing over-densities in a supersonic flow. For example, Burkhart et al. (2015a) examined the observational diagnostics of self-gravity by employing MHD simulations and found that self-gravitating supersonic turbulence can produce density structures that drive the spectral slope even up to positive values. This is in contrast to non-gravitating turbulence where the power is dominated by large-scale structures and decreases toward smaller scales.

One of the important caveats for the above predictions is that the optical depth of the medium can significantly alter the column density/integrated intensity power spectral slope. For example, Lazarian & Pogosyan (2004) and Burkhart et al. (2013) showed that the integrated intensity images of the optically thick medium will have a power spectrum whose slope saturates to a universal value of $-3$ regardless of the presence of shocks, gravity, or magnetic fields. This would likely apply to 12CO(1–0), which is typically optically thick in molecular clouds, and possibly also to 13CO(1–0) and the cold Hi.

5.2. Hi

Theoretical models of heating and cooling in the ISM predict the presence of two main atomic phases, the CNM and WNM (e.g., Field et al. 1969; McKee & Ostriker 1977). With distinct physical properties (density $n$ and kinetic temperature $T_{\rm k}$; $n$ $\sim$ 5–120 cm*-3* and $T_{\rm k}$ $\sim$ 40–200 K for the CNM; $n$ $\sim$ 0.03–1.3 cm*-3* and $T_{\rm k}$ $\sim$ 4100–8800 K for the WNM; e.g., Wolfire et al. 2003), the CNM and WNM would then contribute to the overall turbulent properties of Hi in different ways. Recently, Burkhart et al. (2015c) attempted to disentangle the turbulent properties of the CNM and WNM in Perseus by estimating $\mathcal{M}_{s}$ using two different constraints: (1) the Hi spin temperatures derived by Stanimirović et al. (2014) based on the Arecibo Hi emission and absorption observations and (2) the width of the Hi column density PDF. While the $\mathcal{M}_{s}$ values based on the Hi absorption data range from $\sim$1 to $\sim$40 with a median of $\sim$4, the average $\mathcal{M}_{s}$ deduced from the Hi PDF width is $\sim$1. The authors then suggested that this discrepancy most likely arises from the fact that the Hi absorption data mainly trace the CNM, while the Hi PDF comes from a mix of the CNM and WNM. They also advocated that the Hi-to-H2 transition could also contribute in narrowing the PDF (see also Bialy et al. (2017)). Considering that our Hi column density ($N$(${\mbox{\sc Hi}})$) traces both the CNM and WNM and the CNM fraction is relatively low in Perseus ($\sim$30%; Stanimirović et al. 2014, we conclude that our derived power spectral slope of $-3.23\pm 0.05$ is consistent with a non-gravitating transonic medium on average.

Given that the CNM fraction in and around Perseus is only $\sim$30%, it is not surprising to find that the power spectral slope derived using the opacity-corrected data ($-3.22\pm 0.05$) agrees well with the uncorrected case. However, at the same time, the opacity may impact the measurement of the power spectral slope for molecular clouds with higher CNM fractions, e.g., W43, where the opacity correction increases the Hi mass by a factor of two (Motte et al., 2014; Bihr et al., 2015).

5.3. Dust

Among the four tracers under study, the 2MASS $A_{V}$ shows the shallowest slope of $-2.72\pm 0.12$, much shallower than the 3D Kolmogorov density slope of $-11/3$. Based on various numerical models of interstellar turbulence (Table  3), this implies that the medium traced by the 2MASS $A_{V}$ is on average self-gravitating and supersonic. For example, Burkhart et al. (2015a) examined how the slope of a column density power spectrum changes as a function of time by running simulations of non-gravitating and self-gravitating MHD turbulence with varying sonic and Alfvénic Mach numbers ($\mathcal{M}_{\rm s}$ and $\mathcal{M}_{\rm A}$) and found that the slope drastically changes as the model cloud evolves. The slope is initially $-11/3$ and becomes increasingly shallow for supersonic turbulence. Once gravity turns on, the slope increases beyond the purely supersonic case and eventually becomes positive. In their simulations, the slope values between $-2.7$ and $-2.5$ indeed appear during the self-gravitating supersonic phase (e.g., Figure 10 of Burkhart et al. 2015a). Several independent studies of $A_{V}$ PDF also support our conclusion: a mix of log-normal and power-law shapes has been found for Perseus, which has been mainly interpreted as the presence of both supersonic turbulence and gravity (e.g., Kainulainen et al. 2011; Burkhart et al. 2015a; Stanchev et al. 2015).

In our analysis, we essentially use the 2MASS $A_{V}$ as a tracer of total hydrogen abundance based on the assumption that interstellar dust and gas are well mixed. This mixing of dust and gas has been inferred from the strong correlation between the color excess and the total hydrogen column density, $E(B-V)$/($N$(Hi) + 2$N$(H2)) = 1.7 $\times$ 10*-22* mag cm2 (e.g., Bohlin et al. 1978; Rachford et al. 2009). With the total-to-selective extinction ratio $R_{V}$ = 3.1, the typical value for the diffuse ISM (e.g., Mathis 1990), the correlation becomes $A_{V}$/($N$(Hi) + 2$N$(H2)) = 5.3 $\times$ 10*-22* mag cm2. This dust-to-gas ratio, however, changes in different environments, e.g., $R_{V}$ is higher with $\sim$4–6 in dense molecular clouds, resulting in higher dust-to-gas values (e.g., Cardelli et al. 1989; Fitzpatrick 1999). To assess the impact of the dust-to-gas ratio variation on our power spectral slope, we then perform a simple test by taking any pixel with $A_{V}$ $\geq$ 5 mag and replacing it with the half of the original value, e.g., the pixels with $A_{V}$ = 6 mag are replaced with 3 mag. This test is to simulate the case where the dust-to-gas ratio increases in dense regions by a factor of two, and we find that the power spectral slope of the simulated $A_{V}$ image is almost the same as before, $-2.72\pm 0.15$, suggesting that our interpretation of the self-gravitating and supersonic medium traced by the 2MASS $A_{V}$ is likely robust.

Another interesting thing to note is that the 2MASS $A_{V}$, which traces both the atomic and molecular media, shows the slope that is distinctly different from the Hi slope. In terms of mass, the neutral medium of Perseus is dominated by Hi (Hi and H2 mass of $\sim$2 $\times$ 104 M*⊙* and $\sim$6 $\times$ 103 M*⊙* respectively; calculated over the Hi coverage in Figure 1), and a naive expectation would then be that the 2MASS $A_{V}$ slope is similar to the Hi slope. Our finding of the shallower slope for $A_{V}$ hence suggests a strong influence of gravity, which has in fact been noted in several numerical studies of interstellar turbulence. For example, Burkhart et al. (2015a) showed that self-gravitating high-density cores on small scales have a significant impact on the global turbulence statistics of the model cloud (e.g., column density PDF and the SPS) and found that the introduction of a delta-function like density profile on small scales can mimic self-gravity in the turbulence power spectrum.

5.4. 12CO(1–0) and 13CO(1–0)

As for 12CO, we find a relatively steep slope of $-3.08\pm 0.08$, which is consistent with $-3$, implying that the 12CO(1–0) emission in Perseus is likely subject to the optical depth effect. The presence of self-absorbing medium makes a substantial impact on the slope of an intensity power spectrum, as shown by several authors including Lazarian & Pogosyan (2004) and Burkhart et al. (2013). In particular, Lazarian & Pogosyan (2004) predicted that absorption induces a universal slope of $-3$ for intensity fluctuations, and Burkhart et al. (2013) numerically confirmed this prediction by analyzing 13CO(2–1) power spectra of MHD simulations with varying sonic and Alfvénic Mach numbers. In the case of optically thick emission with $\tau\gg 1$, Burkhart et al. (2013) found that the integrated intensity spectral slope saturates to $-3$ regardless of $\mathcal{M}_{\rm s}$ and $\mathcal{M}_{\rm A}$ values. On the other hand, for the mildly optically thick medium with $\tau\sim 1$, the behavior was mixed: super-Alfvénic turbulence shows an increasingly shallow slope as $\mathcal{M}_{\rm s}$ increases (from $-3$ to $-2$ for $\mathcal{M}_{\rm s}$ $\sim$ 0.4–8), while the slope is always around $-3$ for sub-Alfvénic turbulence. For Perseus, we find that 12CO is somewhat optically thick with $\tau$(12CO) $\sim$ 0.3–1 across the cloud, based on the curve-of-growth analysis of 13CO by Pineda et al. (2008) (Table 4). These estimates are, however, for individual dark and star-forming regions on $\sim$10 pc scales: $\tau$(12CO) will be most likely $\gg$ 1 for smaller and denser cores on $\sim$0.07 pc scales (12CO pixel size). Indeed, Pineda et al. (2008) found that $I$(12CO) becomes saturated at $A_{V}$ $\gtrsim$ 4 mag in Perseus on $\sim$0.4 pc scales (2MASS $A_{V}$ pixel size), indicating that the emission becomes optically thick. The same argument applies for 13CO. While the measured slope of $-2.88\pm 0.07$ is not as consistent with $-3$ as it is for 12CO, $\tau$(13CO) is in fact not drastically different from $\tau$(12CO) (Table 4): only two regions, B5 and Westend, show $\tau$(12CO) more than a factor of two higher than $\tau$(13CO). In addition, the 13CO(1–0) emission was found to saturate at $A_{V}$ $\gtrsim$ 5 mag (Pineda et al., 2008). Finally, we note that 12CO(3–2) and 13CO(2–1) are also likely subject to the optical depth effect, based on Sun et al. (2006) who measured the slopes of $-3.15\pm 0.04$ and $-3.03\pm 0.14$ for the two transitions.

Other than the opacity, the 12CO(1–0) and 13CO(1–0) emission could also be subject to the effects of CO chemistry, e.g., formation (CO abundance sharply drops in diffuse regions with $A_{V}$ $\lesssim$ 1–2 mag due to insufficient dust shielding; e.g., Wolfire et al. 2010) and depletion (CO is frozen onto dust grains in dense regions with $A_{V}$ $\gtrsim$ 5–10 mag; e.g., Bergin et al. 1995). To evaluate the impact of CO chemistry on our power spectrum analysis, we then perform a similar test on the 2MASS $A_{V}$ image as we do for the dust-to-gas ratio variation in Section 5.3. Specifically speaking, we replace any pixels with $A_{V}$ $<$ 2 mag with zero values and re-measure the power spectral slope to mimic the effect of CO formation. Similarly, all pixels with $A_{V}$ $>$ 5 mag are fixed to 5 mag to simulate CO depletion in dense regions. In the CO depletion experiment, the power spectral slope of $-2.77\pm 0.13$ is comparable to the original value of $-2.72\pm 0.12$. In addition, we found that the new slope of $-2.39\pm 0.07$ for the CO formation test is reasonably close to the original value, implying that CO formation and depletion likely do not significantly affect our SPS results. While this conclusion is consistent with Padoan et al. (2006a), who performed essentially the same tests on Taurus, another low-mass star-forming region, we note that more molecular clouds with diverse properties should be further examined to confirm the conclusion. All in all, our results suggest that particular attention is needed to use 12CO(1–0) and 13CO(1–0) as a probe of the turbulent properties of molecular gas, since the emission may suffer from the effects of opacity and chemistry.

5.5. Comparison to Previous Studies

In this section, we make a comparison with previous SPS studies to place our results into a larger context. For the comparison, we consider density fluctuations in various Galactic environments (e.g., non-star-forming clouds and CNM-dominated regions) and focus on Hi and dust power spectra due to potential opacity and chemical effects with the molecular gas tracers. The details on the previous studies, including the measured SPS slopes and spatial scales, are presented in Table 5.

First of all, we find that our results are generally consistent with the previous studies for the Milky Way: the Hi SPS shows a Kolmogorov-like slope (from $-4$ to $-3.6$), while the dust SPS has a much shallower slope (from $-2.9$ to $-2.7$). This finding suggests that on average, Hi emission traces a non-gravitating transonic/subsonic medium, while dust probes a self-gravitating supersonic medium. The only exceptions are Deshpande et al. (2000) and Martin et al. (2015), where the Hi SPS slopes were found to be shallower than $-3$. In the case of Deshpande et al. (2000), their data strictly trace the CNM seen in absorption against Cassiopeia A and Cygnus A, and the measured slope of $-2.75\pm 0.25$ in fact agrees with the expectation for the supersonic CNM. On the other hand, the analysis by Martin et al. (2015) for intermediate Galactic latitudes is based on Hi emission observations (tracing both the CNM and WNM), and whether their shallow slopes of $\gtrsim$ $-3$ result from a substantial amount of the CNM remains to be further examined.

Interestingly, no break has yet been observed in the SPS analyses probing the spatial scales from $\sim$0.01 pc to a few tens pc. One possible explanation is that interstellar turbulence is driven on the scales larger than individual molecular clouds. Possible turbulence drivers then include galactic-scale thermal, gravitational, and magnetorotational instabilities, accretion of circumgalatic material, as well as spiral shock waves (e.g., Dobbs & Bonnell 2007; Brunt et al. 2009; Krumholz & Burkhart 2016). Small-scale stellar feedback such as outflows and bubbles is likely not the dominant source for turbulence (e.g., Nestingen-Palm et al. 2017), considering comparable SPS slopes between non-star-forming (e.g., Polaris flare, Ursa Major cirrus clouds, and MBM16) and star-forming regions over a wide range of length scales (e.g., Perseus and the SMC).

In particular, for Perseus, Padoan et al. (2009) investigated both density and velocity fluctuations in NGC1333, an active star-forming and outflow-driving region in Perseus. If turbulence is primarily driven by outflows on the scale of individual stellar clusters, the power spectral slope should flatten beyond the energy injection scale (e.g., Matzner 2007), which Padoan et al. (2009) estimated to be 0.3 pc for NGC1333. Our measured power spectrum for 13CO(1–0) does not show any characteristic length scale. Finally, the clear power-laws down to $\sim$0.01 pc scales suggests that turbulence dissipation, likely via viscous dissipation and ambipolar diffusion (e.g., Hennebelle & Falgarone 2012), occurs on the scales smaller than 0.01 pc.

Lastly, we consider how the SPS performs overall as a diagnostic tool for characterizing turbulence over other statistical techniques. Recently, Boyden et al. (2016) evaluated several astrostatistics (e.g., Principal Component Analysis (PCA), spectral correlation function, SPS, etc.) as a probe of stellar wind feedback by applying these statistics to synthetic 12CO(1–0) cubes generated from MHD simulations. In their experiments, the SPS showed a strong response to time evolution, but ambiguous behavior to wind activities and very little response to magnetic fields (although due to opacity effects, 12CO(1–0) may not be an optimal tracer for such a study). In addition, by having two different scales of energy injection in MHD simulations, Yoo & Cho (2014) demonstrated that the small-scale driving could have only a limited influence on the density fluctuations. These results suggests that the SPS maybe not a good diagnostic tool to study small-scale stellar feedback such as stellar winds, and thus utilizing a suit of statistical tools whose individual strengths probe different observable characteristics would be necessary to fully characterize turbulence in the multi-phase ISM.

6. Conclusion and future work

In this study we presented a comparison of the two dimensional spatial power spectrum (SPS) derived from four distinct ISM tracers for the Perseus molecular cloud. In our derivation of the SPS, we take median value of the pixel distribution in each annulus as opposed to the standard practice of taking the mean to mitigate the effects of the Gibbs phenomenon. We find the Hi data produce the steepest slopes with values of $-3.23\pm 0.05$ and $-3.22\pm 0.05$ for the uncorrected and corrected Hi, respectively. The dust gives the shallowest slope values of $-2.72\pm 0.12$. In the case of 12CO(1–0) and 13CO(1–0) we measure slopes of $-3.08\pm 0.08$ and $-2.88\pm 0.07$, respectively. The comparison between the relative slopes of the tracers reveals several important characteristics about the turbulent environment in Perseus:

• •

Our derived slope value for the Hi in Perseus suggests that the Hi is largely non-gravitating and, in general, transonic. The consistency of the slope between the opacity corrected and uncorrected column density images is not surprising since the CNM fraction in Perseus is only about 30%. This may not be the case for molecular clouds with significant CNM fractions.

• •

We assume the dust dependably traces the gas in Perseus. It is therefore interesting to see the dust have such a shallow slope as compared to the Hi (since Hi dominates the mass of Perseus). Simple tests which increases the dust-to-gas ratio in dense regions of the dust image reveal a similar slope value, which according to numerical studies, indicates the dust in Perseus is influenced by self-gravity due to the relatively shallow slope.

• •

While the dust is rather robust to small-scale dust-to-gas variations, we find the 12CO(1-0) and 13CO(1–0) power spectra are particularly susceptible to opacity effects — confirming theoretical prediction of a slope of \sim$$-3 when the SPS is applied to an optically thick medium. Furthermore, we show through simple models of CO formation/depletion that such chemical processes do not significant affect our power spectrum analysis for Perseus specifically. This conclusion may not be applicable for other molecular clouds with more complicated chemistry. We conclude careful consideration is needed when utilizing 12CO(1-0) and 13CO(1-0) intensity maps as a probe for turbulent environments due to possible opacity and chemical effects.

• •

Previous results of dust power spectra in diffuse, high-latitude clouds agree well with our dust slope. The similar slope values measured between regions of active star formation (i.e., Perseus) and non-star forming (Polaris Flare and Ursa Major cirrus clouds) hint that small-scale stellar feedback is relatively unimportant in regulating the turbulent environment traced by interstellar dust. A comparison of previous Hi results for the SPS derived over a diverse set of environments (e.g. star forming vs. non-star forming and CNM dominated regions) shows our Hi slope is in good general agreement, and is indicative that Hi traces a non-gravitating transonic to supersonic medium in Galactic environments.

• •

The lack of a break in any of the power spectra signifies we detect neither the injection nor the dissipation scale in Perseus for any of our tracer. Additionally, the absence of a break in any previous Galactic SPS analysis suggests large-scale turbulence drivers, such as galactic-scale instabilities driven by gravitational, magnetorotational or supernovae dominate over small-scale drivers like stellar outflows.

B.B. is supported by the NASA Einstein Postdoctoral Fellowship and the ITC Postdoctoral Fellowship at the Harvard-Smithsonian Center for Astrophysics.