The Astrophysical Journal, 595:824–841, 2003 October 1 # 2003. The American Astronomical Society. All rights reserved. Printed in U.S.A.
INTRINSIC, OBSERVED, AND RETRIEVED PROPERTIES OF INTERSTELLAR TURBULENCE Christopher M. Brunt,1,2 Mark H. Heyer,3 Enrique Va´zquez-Semadeni,4 and Ba´rbara Pichardo5 Received 2003 May 8; accepted 2003 June 6
ABSTRACT We generate synthetic observations of three-dimensional, self-gravitating MHD simulations of the interstellar medium (ISM) to evaluate the ability of principal component analysis (PCA) to measure the scale dependence of turbulent velocity fluctuations in molecular clouds. Scaling exponents, , observationally obtained from the coupled characteristic scales for line profile variability in velocity, v, and in space, L, where v / L , are compared with the intrinsic scaling exponents of the MHD velocity fields. We determine the approximate structure function order at which PCA operates in order to then verify a previously established calibration of the PCA method. We also analyze the statistical properties of projected velocity line centroid fields, including effects of intermittent velocity fluctuations, density inhomogeneity, and opacity, and examine the relationship of the projected two-dimensional statistics to the intrinsic three-dimensional statistics. Using PCA, we infer steep three-dimensional energy spectra in the molecular ISM, generally steeper than can be accounted for by Kolmogorov turbulence or possibly even shock-dominated turbulence. Subject headings: ISM: clouds — ISM: kinematics and dynamics — methods: statistical — MHD — radio lines: ISM — turbulence injection acting on a range of scales, and so will in general be dependent on the energy-injection spectrum in addition to the transfer of energy between scales; in this case still provides useful information on the distribution of kinetic energy as a function of scale within the clouds. In the absence of internal driving sources, measurement of provides fundamental physical information on the nature of interstellar turbulence, which is expected to differ substantially from turbulence experienced in terrestrial situations. The exponent may be obtained from three-dimensional numerical simulations of turbulence under a wider range of conditions than is accessible analytically. Ossenkopf & Mac Low (2002) have measured 2 (or perhaps slightly larger) from simulations of driven magnetohydrodynamic turbulence in three dimensions. Vestuto, Ostriker, & Stone (2003) measured from three-dimensional simulations with a range of physical conditions, resulting in 2 but with large variations, including many cases of > 2 (significantly in excess of the Burgers index) for both shear and compressible modes. However, these measurements are only over a limited spatial dynamic range and must be considered preliminary estimates. Brunt & Heyer (2002a, hereafter Paper I) established a means of estimating from spectral line imaging observations of CO isotopes based on the PCA method of Heyer & Schloerb (1997), using a calibration obtained from simulated observations of a series of simple molecular cloud models with varying . Brunt & Heyer (2002b) reported a mean ¼ 2:17 0:31 (with a distribution preferentially skewed toward higher rather than lower) obtained from PCA of 23 molecular regions in the outer Galaxy. Brunt (2003b) has analyzed many more observations that confirm the prevalence of higher , according to the PCA calibration. The realm of applicability of PCA is currently limited by the molecular cloud models used to establish the method. Padoan, Goodman, & Juvela (2003) have shown, for example, that the stochastic models used to calibrate the PCA method in Paper I can be observationally ruled out (compared with more realistic models) using the spectral
1. INTRODUCTION
Recovery of information on the nature of turbulence in the interstellar medium (ISM) is both very important and very difficult. On physical grounds, turbulence is a significant agent that determines the structure and evolution of interstellar clouds (see, e.g., the review by Va´zquez-Semadeni 2002). On observational grounds, the complexity of the nonlinear transformation of the physical fields into the observable radiation fields is a direct result of the presence of turbulence and complicates recovery of information (Ballesteros-Paredes & Mac Low 2002). A number of explorations of the complexity of this transformation have been conducted (e.g., Falgarone et al. 1994; Padoan et al. 1998; Brunt 1999; Pichardo et al. 2000; Lazarian & Pogosyan 2000; Ballesteros-Paredes & Mac Low 2002; Ossenkopf 2002; Ossenkopf & Mac Low 2002). The majority of these studies identify the velocity field as playing a significant or dominant role in the creation of the observed emission structures on the position-position-velocity axes. An important quantitative measure of the state of turbulence in the ISM is the exponent describing the scale dependence of turbulent velocity dispersion, often quoted as the spectral slope, , of the energy spectrum [EðkÞ / k ]. For incompressible turbulence 5=3 (Kolmogorov 1941), while for shock-dominated turbulence 2 (Burgers 1974). Molecular clouds may have many sources of energy 1 National Research Council, Herzberg Institute of Astrophysics, Dominion Radio Astrophysical Observatory, P.O. Box 248, Penticton, BC V2A 6K3, Canada;
[email protected]. 2 Department of Physics and Astronomy, University of Calgary, 2500 University Drive Northwest, Calgary AB, T2N 1N4, Canada. 3 Five College Radio Astronomy Observatory and Department of Astronomy, Lederle Research Building, University of Massachusetts, Amherst, MA 01003. 4 Instituto de Astronomı´a, Universidad Nacional Auto ´ noma de Me´xico, Campus Morelia, Apartado Postal 72-3 (Xangari), Morelia, Michoaca´n 58089, Mexico. 5 Instituto de Astronomı´a, Universidad Nacional Auto ´ noma de Me´xico, Apartado Postal 70-264, Me´xico DF 04510, Mexico.
824
INTERSTELLAR TURBULENCE correlation function (Rosolowsky et al. 1999; Padoan, Rosolowsky, & Goodman 2001). There is concern therefore that the previously established calibration may be subject to some uncertainty. In this paper we extend the assessment of the method to include more physically realistic, observed models obtained from MHD simulations of turbulence. This also requires a more complex analysis of turbulent velocity field statistics. In x 2 the numerical MHD simulations are introduced and the statistical measures needed to characterize the MHD velocity fields are described. The transformation of the physical fields onto the observationally accessible radiation fields is described in x 3. In x 4 we examine the statistical properties of projected velocity line centroid fields. Section 5 provides a brief summary of the PCA method. In x 6 we apply PCA to synthetic observations of the numerical simulations and examine the ability of PCA to recover statistical information describing the intrinsic velocity fields.
2. MHD FIELDS AND THEIR STATISTICAL PROPERTIES
2.1. The Simulated Velocity and Density Fields The numerical simulation analyzed in this paper has been previously described by Pichardo et al. (2000; see also Passot, Va´zquez-Semadeni, & Pouqet 1995) and by Lazarian et al. (2001). The MHD equations were solved under conditions characteristic of the warm neutral medium within a cubic box of 300 pc on each side, at a resolution of 128 grid points per dimension. Physical inputs included selfgravity, magnetic fields, parameterized heating and cooling, and modeled star formation (localized heating over a limited time, triggered when the density field exceeds an imposed threshold within a converging flow) that maintains the turbulence. The energy-injection results in expanding shells of material around the sites of star formation (Pichardo et al. 2000). Expanding shells of neutral material around H ii regions are commonly observed in high-resolution H i data (Kothes & Kerton 2001). The power spectra of the density and velocity fields were not well-developed power laws, owing to the low spatial dynamic range of the simulation. To maintain well-defined input statistics with which to compare the observable quantities, the power spectra of the density and velocity fields were forced into obeying a power-law form, as described in Lazarian et al. (2001). The fields were transformed into Fourier space, and each amplitude was modified so that a power law is obeyed while retaining the phase structure. This modifies the details of the density-velocity coupling only slightly. Lazarian et al. (2001) used the simulations to test the emissivity spectral index method of Lazarian & Pogosyan (2000). The spectral index, , of the energy spectrum [EðkÞ; see x 2.2] was forced onto values of 1.0, 1.5, 2.0, 2.5, and 3.0. The power spectrum of the density field was forced onto the following form: PðkÞ / k4 . In this paper we have taken advantage of this existing set of simulated fields to investigate the PCA method more thoroughly. We are mostly interested in analyzing fields for which there exists a more realistic coupling between the density field and velocity field (see Lazarian et al. 2001 for details) than was present in the fractional Brownian motion (fBm) fields used previously (Paper I). As the basic physical ingredients of the simulations are not critically dependent
825
on the overall scaling of the models and it is more important to compare the MHD fields with the fBm fields under identical observational circumstances than to maintain physical consistency with the initial model parameters, we have rescaled the original physical parameterization of the simulated fields into conditions more typical of molecular clouds (as detailed in x 3). The simulations used here thus represent a significant improvement over the previously used fBm fields. The most important components of the numerical fields, for our purposes here, are the presence of intermittent velocity fluctuations and the density-velocity coupling, features that were not present in the models of Paper I. In the analysis that follows, we also consider phaserandomized versions of the velocity fields. For these ‘‘ dephased ’’ fields, we randomize the phases of the wavevectors in Fourier space while retaining the spectral slope previously enforced; the fields thus produced are then equivalent to fBm fields. 2.2. Generalized Scaling Laws The statistical properties of a turbulent velocity field are only partially described by the slope of the energy spectrum EðkÞ in Fourier space, EðkÞ / k ;
ð1Þ
where EðkÞ is the angular integral of the power spectrum [EðkÞ ¼ 4k2 hPðkÞi ]. The energy spectrum description translates into a real space description of the second-order velocity fluctuations as D E ð2Þ jvj2 / L1 ; measured over a scale L. Typically, the exponent is rewritten in terms of , where ¼ ð 1Þ=2 ; so that, in terms of the rms velocity fluctuation, D E1=2 / L : jvj2
ð3Þ
ð4Þ
In this work we consider a more complete statistical description of a turbulent velocity field that can account for the presence of intermittent velocity fluctuations, which manifest themselves as broad (exponential) wings in the probability distribution function of velocities and as modifications to the scaling properties of velocity fluctuations. These properties are also associated with the presence of intermittency, i.e., sparseness in the spatial and temporal dissipation of turbulent energy. Intermittency has a narrower definition than simply the presence of intermittent velocity fluctuations: for the velocity fields considered here, the ‘‘ intermittent ’’ velocity fluctuations arise not only from energy dissipation but also from energy injection by the modeled star formation (which is also sparsely distributed in space and time); obviously the sites of injection and dissipation are related in some way, however (Avila-Reese & Va´zquez-Semadeni 2001). In this paper the adjective ‘‘ intermittent ’’ should be interpreted in the most general sense. Accordingly, a turbulent velocity field is more generally described by the scaling properties of the qth-order velocity fluctuation, denoted as the qth-order structure function, Sq
826
BRUNT ET AL.
Vol. 595
(see Frisch 1995; Boldyrev, Nordlund, & Padoan 2002), Sq ðLÞ ¼ hjvjq i / LðqÞ ;
ð5Þ
parameterized by the function ðqÞ. Since measurements are typically made at a single value of q (or some finite number of values of q), it is convenient to write these measurements as subscripted quantities q [e.g., 2 ¼ ð2Þ]. Similar to the form of the rms velocity fluctuation above, the ‘‘ qth–root mean qth-power ’’ velocity fluctuation is written ðhjvjq iÞ1=q / Lq =q :
ð6Þ
The generalized version of the exponent is defined as q ¼ q =q :
ð7Þ
Note that 2 ¼ ð 1Þ=2 still holds. The variation of q with q provides insight to the flow of energy within a turbulent fluid (Boldyrev et al. 2002). A non- or weakly intermittent field is statistically more homogeneous, i.e., there are less extreme differences in the magnitude of velocity fluctuations sampled over scale L. A nonintermittent velocity field is defined by a linear form for q : q ¼ q ;
ð8Þ
where the unsubscripted is a constant. Then q ¼ ;
ð9Þ
¼ ð 1Þ=2 :
ð10Þ
and of course
The complete moment dependence of the structure functions is necessary to describe the velocity fluctuations within a general velocity field. Moreover, different analysis methods may operate at different effective orders. By analyzing the sensitivity of measurable exponents, e.g., the PCA exponent, to each of the intrinsic q -exponents, we can identify the effective order at which a given analysis method works. 2.3. Statistics of the MHD Velocity Fields The statistics of the original velocity fields provide a basis on which to compare with the output of subsequent analyses and to gauge the transformation of statistics by the observational process. A coarse structure function, Sq ðLÞ, was calculated from each velocity field for q ¼ 0:5, 1.0, 1.5, 2.0, 2.5, and 3.0 from which q is derived. The coarse structure functions were calculated as follows: each velocity field, vz , was partitioned into an ensemble of cubical boxes of L pixels on a side, over which vz is averaged. This procedure was carried out for dyadic values of L (1, 2, 4, 8, 16, 32, 64), where L ¼ 1 corresponds to the original field. For each coarsened version of the field, the absolute differences along the cardinal directions between adjacent coarse pixels was measured. The ensemble of such measurements (i.e., for all possible adjacent pixels) was raised to the power q and averaged at each L. The coarse structure function at qth order is then Sq ðLÞ ¼ hjvL ðrÞ vL ðr þ LÞjq i ;
ð11Þ
where vL is a coarse-grained field. From these measurements
Fig. 1.—Measured coarse structure functions Sq(L) for the MHD field with ¼ 2, in original form (left) and dephased form (right). The absolute scaling of the ordinate is arbitrary, and the measured points have been offset in this direction to improve legibility. From bottom to top, q ¼ 0:5, 1.0, 1.5, 2.0, 2.5, and 3.0.
the exponents q were obtained via the proportionality Sq ðLÞ / Lq :
ð12Þ
This procedure (at q ¼ 2) bears similarities to the more sophisticated D-variance method (Stutzki et al. 1998); the D-variance could easily be extended to operate at multiple q. The reader is referred to Brunt (1999) for more information on the coarse structure function at q ¼ 2. The same measurements were calculated for the dephased version of the fields in which the Fourier space phases have been randomized while retaining the spectral slope (). This dephasing process effectively converts the original field into fBm, which is nonintermittent, because the sharp velocity discontinuities in the original field (encoded in the Fourier space phase structure) are removed. Figure 1 shows example coarse structure functions, Sq ðLÞ, obtained for the original and dephased fields with ¼ 2. The q and q are obtained from the log-log slopes of Sq ðLÞ at the appropriate q, avoiding the turnover region in Sq ðLÞ at large L. The turnover at large lags is characteristic of both periodic fields (increasing separations ultimately decrease) and of the method of measurement (it occurs also in the D-variance technique for example). Its origin is simple to understand: as the outer scale of variability is exceeded, the variance decreases, just as the variance of Gaussian noise decreases upon averaging. Table 1 lists the measured exponents for all the velocity fields used in this study. The variation of q with q for ¼ 2 is shown in Figure 2 for the original and dephased fields. For the original field q is a concave function of q and diagnoses the presence of intermittent velocity fluctuations. In contrast, q for the dephased fields is a linear function of q. The second-order exponents (2 ,2 ) are the same for both original and dephased fields since is unchanged by the dephasing process. Note that for the ¼ 1 field, 2 0 (see Table 1). This scaling exponent describes the fact that the typical fluctuations (at second order) are independent of the scale over which they are measured; i.e., the field is a ‘‘ 1/f noise.’’ The total dispersion of velocities in a ¼ 1 field increases as log L. The exponents q for all original fields are summarized in visual form in Figure 3, along with the histograms of the velocity values in each field, before and after dephasing. Strong intermittent velocity fluctuations detected by the moment-dependence of scaling exponents is directly
No. 2, 2003
INTERSTELLAR TURBULENCE
827
TABLE 1 Scaling Exponents for the Original (O) and Dephased (D) Velocity Fields
O/D
0.5
1
1.5
2
2.5
3
1.0.............. 1.5.............. 2.0.............. 2.5.............. 3.0.............. 1.0.............. 1.5.............. 2.0.............. 2.5.............. 3.0..............
O O O O O D D D D D
0.20 0.02 0.41 0.01 0.61 0.01 0.76 0.02 0.86 0.01 0.01 0.02 0.22 0.01 0.44 0.01 0.65 0.01 0.80 0.01
0.15 0.02 0.37 0.02 0.57 0.02 0.74 0.02 0.86 0.02 0.01 0.03 0.22 0.02 0.45 0.01 0.65 0.01 0.81 0.01
0.08 0.03 0.31 0.02 0.52 0.02 0.70 0.02 0.84 0.02 0.01 0.03 0.22 0.02 0.45 0.01 0.66 0.01 0.81 0.01
0.01 0.03 0.22 0.02 0.45 0.01 0.65 0.01 0.81 0.02 0.01 0.03 0.22 0.02 0.45 0.01 0.65 0.01 0.81 0.01
0.11 0.03 0.13 0.02 0.37 0.01 0.60 0.01 0.78 0.01 0.01 0.03 0.22 0.01 0.45 0.01 0.65 0.01 0.81 0.01
0.19 0.03 0.05 0.02 0.30 0.01 0.54 0.01 0.74 0.01 0.01 0.03 0.22 0.01 0.45 0.01 0.65 0.01 0.81 0.01
associated with broad wings in the velocity histograms. These broad wings are removed by the dephasing operation, as is the moment dependence of the scaling exponents. The fields with lower show more pronounced deviations of q from 2 , as the small-scale structure of the field is emphasized by the shallower spectral slope. The dephased fields all have q ¼ 2 , as can be seen in Table 1. There is a small discrepancy between the measured 2 exponents and the ideal 2 ¼ ð 1Þ=2 expected from the spectral slope of the energy spectrum. This difference can be reduced for fields much larger than those considered here. The measured 2() relationship for the MHD velocity fields is shown in Figure 4. We use this measured 2() relationship to define a ‘‘ metric ’’ for the ensemble of velocity fields, which, as noted
Fig. 2.—Measured scaling exponents q and q ¼ q =q for the coarse structure functions in Fig. 1. The departure of q from a linear relationship with q diagnoses the presence of intermittent velocity fluctuations. The dephasing operation removes the dependence of q on q (i.e., q ¼ qq Þ. The exponents 2 and 2 are unaffected since the (second-order) power-spectrum slope is unmodified.
previously, are identical in their phase structure and differ only in the imposed spectral slope, . Conversion of one velocity field into another, of different , involves powerlaw filtering the Fourier amplitudes only. Conversion of, for example, a velocity field of spectral slope to a new spectral slope þ 1 then simply moves the 2 exponent along this metric. This is important for the analysis of line centroid velocities in x 4. 3. THE OBSERVED SIMULATIONS
3.1. Simulated Molecular Line Observations An observer does not have access to the threedimensional density and velocity fields in the interstellar medium. Instead, the accessible field is the radiation field, Tðx; y; vz Þ, where x and y are projected positions on the sky and vz is a spectral axis. In the case of spectral line emission from the abundant CO isotopes in the molecular ISM or H i 21 cm line emission from the atomic ISM, the data cube, Tðx; y; vz Þ, can be extremely complex. For the simulated molecular spectral line observations, we replicated the input conditions utilized in Paper I; also, the method used in Paper I was used to calculate the emergent intensity fields. In brief, we first determined the excitation temperature of the line and line-center opacity at each grid point from a non-LTE (i.e., large velocity gradient) calculation, which accounts for local radiative trapping (Scoville & Solomon 1974). Following this, the emergent intensity in each velocity channel was calculated for each position by transferring the radiation through the model grid. The reader is referred to Paper I for full details. A pixel scale of 0.156 pc was assumed, giving a total size of 20 pc for the 1283 cube, and we used a mean H2 number density of 100 cm3. The global velocity dispersion of each velocity field was set to 1 km s1. The spectral resolution of the observation was 0.05 km s1. We observed 13CO (J ¼ 1 0), including excitation calculations up to J ¼ 5, assuming a constant abundance value of 1:25 106 relative to H2 and a constant kinetic temperature of Tk ¼ 20 K. We also observed 12CO (J ¼ 1 0), assuming a constant abundance of 104, to investigate the effects of saturation on the PCA method. As an ideal tracer we observe ‘‘ XCO ’’ (J ¼ 1 0), which is simply the 13CO (J ¼ 1 0) transition with the abundance set to 1012 and the H2 density set to 106 cm3 (i.e., the high densities ensure that the transition is thermalized, while the low abundance ensures that the lines do not saturate). For 12CO and 13CO
828
BRUNT ET AL.
Vol. 595
Fig. 3.—(a) Summary of the measured exponents q for the five original velocity fields, in relation to the second-order exponent 2. (b)–( f ) Histograms of velocity fields scaled to the global standard deviation for ¼ 1, 1.5, 2, 2.5, and 3. Filled circles represent the original velocity fields, and open circles denote the dephased velocity fields. Intermittent velocity fluctuations produce the non-Gaussian wings of these distributions.
the excitation was always subthermal, because of the low densities, and the optical depths were generally high because of the resulting overpopulation of the ground state. The density field was either that generated by the MHD simulation [subsequently forced onto a PðkÞ / k4 power spectrum] or a uniform density field. A density power spectrum of PðkÞ / k4 is steeper than the spectrum inferred to be characteristic of molecular clouds (k2.8) as determined by Bensch, Stutzki, & Ossenkopf (2001). (However, measured power spectra of our simulated integrated intensity maps are somewhat shallower, k3:4 for 13CO and k3:0 for 12CO, because of opacity effects, and are thus a little more in accord with the observed values.) Our reason for taking such a steep slope was to avoid ringing artifacts around sharp density contrasts that arise if the spectral slope modification is pushed to too shallow an index. These can cause negative densities directly adjacent to the sharper features in the field.
We also performed simulated observations of the dephased velocity fields for the uniform density case, using 13CO. A summary of the simulated observations carried out under ‘‘ molecular cloud ’’ conditions is given in Table 2. All CO spectral line transitions are J ¼ 1 0. Each set of input TABLE 2 Parameters of the Simulated Observations
Tracer
Abundance
nH2 (cm3)
Density Field
Velocity Field
XCO .................
1012 1012 1.25 106 1.25 106 1.25 106 104
106 106 102 102 102 102
Uniform MHD Uniform Uniform MHD MHD
Original Original Dephased Original Original Original
XCO ................. 13CO
................ ................ 13CO ................ 12CO ................ 13CO
No. 2, 2003
INTERSTELLAR TURBULENCE
829
preferentially lower the abundance of, or may completely destroy, the absolutely less abundant species. No allowance for this is included in our simulated fields. In general, to factor photodissociation effects into simulations would be very difficult: the abundance structure and photodissociation structure will generally be determined by the full history of the cloud. In real clouds, abundance variations will certainly be present, and there may well be no optically thin tracer on which large-scale velocity crowding in low-density material can operate effectively. For self-shielding species, the lines will quickly saturate, masking the inherent multiplicity (e.g., as seen in the higher 12CO observations). Another property of the simulations that is not present in real clouds is the presence of periodic boundary conditions. The velocity field is constrained to reach a maximum fluctuation and then return, causing velocity crowding at the velocity maxima and thus the ‘‘ limb brightening ’’ seen at the line profile edges in the simulated radiation fields. In the real ISM the fluctuations are of course not band-limited at low spatial frequency in this way. The overstructured line profiles of course may alternatively be revealing some crucially important oversight in the physics of interstellar turbulence. 4. VELOCITY LINE CENTROID ANALYSIS Fig. 4.—Variation of the second-order structure function index, 2, with the slope of the energy spectrum, , defining the ‘‘ metric ’’ discussed in the text.
parameters was used to generate simulated observations at each . Examples of the simulated radiation fields are given in Figure 5 (representative position-velocity plots) and Figure 6 (line center channel maps). Clearly, the substantial exponential tails in the lower fields lead to quite unrealistic 12CO and 13CO radiation fields. At higher the positionvelocity maps are reminiscent of those obtained from molecular clouds containing outflows. Note the similarity between the XCO unweighted and density-weighted cases, demonstrating the importance of the velocity field in shaping the observed radiation fields. It is evident also, however, that the structures that are enhanced in XCO when densityweighting is included are also (for the most part) more prominent in the 13CO and 12CO radiation fields. This indicates that when opacity is important, the density field takes on a more important role in shaping the observed fields. On the whole, we believe that the simulated fields contain more spectroscopic structure than is typically seen in real observations. Admittedly, this is a subjective interpretation, and clearly the presence of noise and low spatial dynamic range in real observations may influence our opinion. However, we stress that overstructured line profiles are not particular to the simulation used here: Ossenkopf 2002 commented on this also, for simulations carried out under actual molecular cloud conditions and with no modeled star formation. If turbulence is driven on small scales only, the level of line-profile multiplicity is lowered (Ossenkopf 2002). However, observations otherwise favor large-scale turbulence (e.g., Ossenkopf & Mac Low 2002; Brunt 2003a). Many of the coherent position-velocity features in the simulations are generated over long path lengths in lowdensity material, wherein the projected line-of-sight velocity changes only slowly. In the real ISM, photodissociation will
Projected line centroid velocity fields provide a convenient and potentially powerful way to obtain the intrinsic three-dimensional velocity field statistics in the ISM (e.g., Scalo 1984; Kleiner & Dickman 1985; Miesch & Bally 1994). In this section, we examine some of the observational aspects of this method. 4.1. Line Centroids: Basic Definitions We calculate the true projected mean velocity field as vtrue mean ðx; yÞ ¼
Npix X
vz ðx; y; zÞ=Npix ;
ð13Þ
i¼1
where Npix is the number of spatial pixels along each line of sight (z-direction). We also calculate the density-weighted projected mean velocity: ,N Npix pix X X nH2 ðx; y; zÞvz ðx; y; zÞ nH2 ðx; y; zÞ : vwgt mean ðx; yÞ ¼ i¼1
i¼1
ð14Þ In the analysis below, for comparison, we generate maps of the line centroid velocity from the observed MHD velocity fields: , N N chan chan X X obs vmean ðx; yÞ ¼ Tðx; y; vz Þvz ðx; yÞ Tðx; y; vz Þ ; i¼1
i¼1
ð15Þ where Tðx; y; vz Þ is the observed radiation temperature and Nchan is the number of channels along each line of sight. We analyze the scaling properties of the projected line centroid velocity fluctuations by measuring their power spectra; this allows exploration of the effects of density weighting and other observational factors such as saturation. In what follows, ‘‘ projected mean velocity ’’ and ‘‘ velocity line centroid ’’ are used interchangeably.
830
BRUNT ET AL.
Vol. 595
Fig. 5.—Example position-velocity maps in the simulated radiation fields. The wedges measure radiation temperature in kelvins.
4.2. Overview of Line Centroid Analysis To frame the analysis below, we examine the effects of three-dimensional to two-dimensional projection on velocity field statistics. We distinguish three-dimensional and (projected) two-dimensional properties by appropriate sub-
scripts or brackets. Here we consider only the true projected mean velocity. The velocity field power spectrum in three dimensions has the form P3D ðkÞ / kð3D þ2Þ ;
ð16Þ
No. 2, 2003
INTERSTELLAR TURBULENCE
831
Fig. 6.—Line center channel maps in the simulated radiation fields. The wedges measure radiation temperature in kelvins.
so that the energy spectrum is Z E3D ðkÞ ¼ d P3D ðkÞ / k3D ;
ð17Þ
since d / k2 (i.e., shells in Fourier space). Consider first an unweighted projection of the velocity field over the line-of-sight (z) axis. Generation of a velocity centroid map, by projection over the z-direction, simply selects the wavevectors for which kz ¼ 0, since all others are
oscillatory along the z-direction and will average out (but see below). To see the effect of this we can write, for wavenumbers kx, ky, kz in the x-, y-, z-directions respectively, P3D ðkÞ / ðkx2 þ ky2 þ kz2 Þð3D þ2Þ=2 :
ð18Þ
In two dimensions the line centroid power spectrum is then P2D ðkÞ / ðkx2 þ ky2 Þð3D þ2Þ=2 ;
ð19Þ
832
BRUNT ET AL.
which is obtained simply by setting kz ¼ 0 in equation (18) (see a similar treatment and more discussion in MivilleDescheˆnes, Levrier, & Falgarone 2003). Thus the power spectrum has the same spectral slope in (projected) two dimensions as it does in three dimensions. This is identical to the result derived in real space by Lazarian & Pogosyan (2000) for the projection of the density field (to the column density field). For projected velocity fields it is important to stress that only the transverse component of the velocity field is projected into two dimensions (e.g., Kitamura et al. 1987). This means that velocity fields composed solely of longitudinal modes will result in a vanishing projected line centroid field in two dimensions. To visualize this, consider two types of velocity feature: a rotating shell (purely transverse) and a symmetrically expanding or collapsing shell (purely longitudinal). The rotating shell (unless inclined face-on) always yields a projected variation in centroid velocity as one side is receding while the other side is approaching. Conversely, the expanding/collapsing shell results in no projected line centroid variation. One can consider a more general velocity field being composed of a superposition of these two basic transverse and longitudinal elements with a distribution of spatial scalings and orientations. Evidently, in the case that ISM velocity fields have a significant longitudinal component, use of projected line centroids to infer three-dimensional statistics may be subject to some uncertainty. The projected two-dimensional energy spectrum is Z E2D ðkÞ ¼ d P2D ðkÞ / kð3D þ1Þ / k2D ; ð20Þ since d / k (i.e., annuli in Fourier space). The energy spectra thus steepen, described by 2D ¼ 3D þ 1 ;
ð21Þ
simply because of the change in dimensionality. The modification ! þ 1 corresponds to a power-law smoothing of the velocity fluctuations. The effect on the second-order real-space exponent, ideally [i.e., 2 ð2DÞ ¼ ð2D 1Þ=2], is thus 2 ð2DÞ ¼ 2 ð3DÞ þ 12 (with the above provisos on isotropy understood). See also O’Dell & Castaneda (1987) for discussion of projection effects. For our finite case here, we expect (and show below) that the ! þ 1 transformation moves 2 along the 2() metric (Fig. 4) instead. For an intrinsic pure power law with any 2 ð3DÞ 12, the line centroid structure function will saturate [and yield 2 ð2DÞ ¼ 1], thus providing no discrimination between different intrinsic 2 ð3DÞ 12. ‘‘ Projection smoothing ’’ [2 ð2DÞ ¼ 2 ð3DÞ þ 12] can be understood simply in real space by considering that the small-scale variations in the velocities suffer more suppression (averaging) than the large-scale variations. Conversely, the total dispersion on small scales, e.g., as measured by the observed line widths, will increase upon projection since a dispersion measured over a region l l actually corresponds to a dispersion measured over the much larger scale l l L0 , where L0 is the overall size of the observed region. A particular example of this is seen in the measurements of Ostriker, Stone, & Gammie (1999); their two-dimensional dispersions measured over l l(L0) have a much shallower variation with l than do their three-dimensional dispersions measured over l l l. Recently, Ossenkopf & Mac Low (2002) have
Vol. 595
exploited the differences between projected line centroid dispersions and total dispersions to estimate the effective lineof-sight depth of the Polaris Flare molecular cloud. The importance of depth effects for the interpretation of emission- and absorption-line profiles has been discussed by Pe´ty & Falgarone (2000). The simple power-spectrum projection result is essentially a tautology for the velocity fields considered here, since they are periodic across the cubes. In reality, a given velocity field will be sampled by an observable tracer of some kind, the effective weighting and boundaries of which will not in general lead to exact cancellation of all velocity wavevectors with kz 6¼ 0. For the highest spatial frequency waves (wavelengths much less than the size of the traced region), cancellation is likely to be a reasonable approximation, but the transfer of the lowest spatial frequency waves will clearly be highly dependent on the geometry of the region sampled by the observational tracer. Moreover, real clouds are not band-limited at low frequency as our simulations are, so large-scale (out-of-band) gradients will play a more important role. For highly structured molecular clouds, one may also expect modifications at higher spatial frequencies. As the power spectrum retains the same slope but the energy spectrum, which defines the second-order scaling properties of the field, is modified as a result of the change in dimensionality, velocity line centroid autocorrelation functions and structure functions, which are directly related to the power spectrum, should be interpreted with this in mind; so too should the power-spectrum analysis (or equivalent) of integrated intensity (‘‘ column density ’’) fields. Finally, while only a power-law form of the intrinsic threedimensional power spectrum of velocities has been explicitly considered above, the overall result—that the line centroid power spectrum is identical to a cut through the threedimensional power spectrum with kz ¼ 0—holds in the general case, with our provisos given in the preceding paragraphs.
4.3. Signatures of Intermittent Velocity Fluctuations in Line Centroid Fields From the true velocity line centroid fields (eq. [13]), we have measured 2(2D) using coarse structure functions. We have plotted the measured 2(2D), in relation to the intrinsic three-dimensional exponents, 2(3D), in Figure 7. The results verify the expected projection smoothing [2 ð2DÞ ¼ 2 ð3DÞ þ 12], taking into account our expectations concerning the ‘‘ metric ’’ (Fig. 4) discussed above; i.e., the transformation ! þ 1 that occurs going from three dimensions to two dimensions moves the exponents roughly as 2 ! 2 ð ! þ 1Þ. To search for residual signatures of intermittent velocity fluctuations in the projected line centroid field, we have measured coarse structure functions for q ¼ 0:5, 1.0, 1.5, 2.0, 2.5, and 3.0, as was done for the intrinsic three-dimensional velocity fields. The measured exponents are listed in Table 3; obviously there is a lesser variation of q with q for the projected fields than for the intrinsic fields. As a measure of the amplitude of intermittent velocity fluctuations in the line centroid field compared with the full velocity field, we plot 1(2D)– 2(2D) versus 1(3D)– 2(3D) in Figure 8. The deviations of q from 2, are suppressed during the projection.
No. 2, 2003
INTERSTELLAR TURBULENCE
Fig. 7.—Measured second-order scaling exponents 2(2D) for the centroid velocity fields vs. 2(3D) obtained from the three-dimensional velocity fields. The dashed line shows the idealized relation [2 ð2DÞ ¼ 2 ð3DÞ þ 12, saturating at 2 ð3DÞ ¼ 12] expected from projection smoothing. The solid line shows the relation between 2(2D) and 2(3D) expected from the 2() metric (Fig. 4).
4.4. Observational Aspects of Line Centroid Analysis The above considerations show that use of line centroid statistics has the clear benefit of allowing a relatively straight forward relation to the true intrinsic threedimensional statistics. However, there are a number of observational aspects that must be considered, notably inhomogeneous sampling of the velocity field due to density variations and masking of portions of the velocity field due to opacity. These aspects, which do not generally permit a useful analytic description, may be examined easily with numerical data. Below we use the index of the power spectrum of the velocity fields as a convenient comparator (i.e., P3D / k3D and P2D / k2D ). The indices obey 3D ¼ 3D þ 2 and 2D ¼ 2D þ 1. The indices are conserved upon projection and the indices are not. As a baseline, we first measured the two-dimensional spectral index, 2D, for all the unweighted (true) centroid velocity maps via equation (13). In this case, we expect 2D ¼ 3D , and indeed this is almost identically recovered, as shown in Figure 9. Note that while 2D ¼ 3D , the direct space statistics (see Figure 7) are not equivalent.
833
Fig. 8.—‘‘ Intermittency ’’ indicator 1- 2 in the velocity centroid field (two dimensions) vs. 1- 2 for the full velocity field (three dimensions). The dashed lines have slopes of unity and zero.
Centroid velocity fields where density-weighting (but no radiative transfer) is included (eq. [14]) were then made. The resulting 2D measurements are shown in Figure 9a, in comparison to the true, unweighted case. It is clear that including density weighting results in shallower indices 2D. Examination of the actual power spectra shows that this arises from an increase in the level of small-scale projected velocity structure rather than a suppression of large-scale projected velocity structure. Thus the density field is imprinting small-scale structure on the projected velocity field. We carried out the same procedure using the dephased velocity fields weighted by the original MHD density field; in this case, the physical density-velocity correlations have been removed by the dephasing operation. The measured indices, 2D, are shown in Figure 9b. The 2D closely approximate the unweighted case, apart from a small discrepancy at higher 3D (higher 3D). Thus it is the correlations between density and velocity that are most important—and not the amplitude of the density fluctuations. Finally, we carried out the same analysis for 13CO and 12CO, using the full radiative transfer calculations for the original MHD velocity and density fields. The results of
TABLE 3 Scaling Exponents for the Line Centroid Velocity Fields 3D
2D
0.5
1
1.5
2
2.5
3
1.0............... 1.5............... 2.0............... 2.5............... 3.0...............
2.0 2.5 3.0 3.5 4.0
0.53 0.01 0.71 0.01 0.84 0.01 0.91 0.02 0.95 0.01
0.50 0.01 0.69 0.01 0.82 0.01 0.91 0.01 0.95 0.01
0.47 0.01 0.66 0.01 0.81 0.01 0.90 0.01 0.95 0.01
0.43 0.01 0.64 0.01 0.80 0.01 0.90 0.01 0.95 0.01
0.39 0.01 0.62 0.01 0.79 0.01 0.89 0.01 0.95 0.01
0.36 0.01 0.59 0.01 0.77 0.02 0.89 0.01 0.95 0.01
834
BRUNT ET AL.
Vol. 595
Fig. 9.—Measured indices in three and two dimensions. In each plot, the open circles are measurements obtained from the true mean velocity field (eq. [13]) that are identical to the intrinsic three-dimensional indices. There are error bars on these plots (apart from on the open circles), but the power laws are so well respected that the error bars are not visible. The observed exponents underestimate the intrinsic exponents when density-velocity correlations or significant opacity are present.
these measurements are summarized in Figure 9c and Figure 9d, respectively. Here the measured 2D is further decreased from both the true, unweighted case, and from the density-weighted case (Fig. 9a). This is certainly due to the effects of opacity, as described below. Obviously, 12CO in particular is a poor choice for measuring line centroid fields (and this is in fact universally avoided to our knowledge); use of 13CO is more common. We see here that the moderate opacities in the 13CO lines cause only a slight decrease in 2D relative to the simple densityweighted case. The effects of opacity may perhaps be understood in part by considering that a highly opaque line traces only the front surface of the cloud and thus that the centroid field may be more similar to a two-dimensional cut through the three-dimensional velocity field rather than to an integration through it. Obviously, a cut through the threedimensional field will yield, on average, similar direct-space indices as the full three-dimensional field itself (this corresponds in Fourier space to 2D 3D 1 due to the change in dimensionality only). Indeed, Figure 9d shows that for 12CO, 2D 3D 1 is a reasonable description of the results. This counters the effect of projection smoothing, as now 2 ð2DÞ 2 ð3DÞ for an optically thick line. However,
for opaque lines the real situation is more complicated than this as in fact a range of depths, depending on the velocity structure itself, will be probed. 5. OVERVIEW OF PRINCIPAL COMPONENT ANALYSIS
Principal component analysis is a widely used well-known statistical tool (e.g., Murtagh & Heck 1986). The particular formulation used here is that of Heyer & Schloerb (1997, hereafter HS97). 5.1. Eigenvectors as Basis Vectors Visualization of the radiation fields relies upon projection of the data values onto a series of user-defined basis vectors. Below, we denote the spatial coordinate as ri ¼ ðxi ; yi Þ at the ith position and the jth velocity channel as vj . A simple set of vectors are un ðvj Þ ¼ unj ¼ ðvn vj Þ ;
ð22Þ
where is the usual delta function. The nth vector thus takes a value of unity in the nth channel and is zero elsewhere.
No. 2, 2003
INTERSTELLAR TURBULENCE
Projection of the data onto these vectors via I n ðri Þ ¼
p X
covariance matrix Su ¼ u : ð23Þ
Tij unj
j¼1
produces p images, I n ðri Þ, which are the usual channel maps. The vectors, unj , are orthogonal, so that each projected channel map provides an independent representation of the data. A similar scheme is used to project position-velocity maps. The above projections involve essentially no compression of information, which is desirable in certain respects, but the amount of information present in the channel maps may be difficult to summarize. For turbulent media, the channel map structure is contributed to substantially by the velocity field, in addition than the density and/ or temperature field. Lazarian & Pogosyan (2000) exploit this feature in their emissivity spectral index method. A number of other vectors are in common use that provide a greater degree of compression of information and generate more readily interpretable projected images. When the data are projected onto the vector w0j ¼ 1 (at all j), the zeroth moment is generated I w0 ðri Þ ¼
p X
Tij w0j ¼
j¼1
p X
Tij ;
ð24Þ
j¼1
and I w0 ðri Þ is equal to the integrated intensity when multiplied by the spectroscopic channel width. Projection of the data onto the vector w1j ¼ vj produces I w1 ðri Þ ¼
p X
Tij vj :
ð25Þ
j¼1
When I w1 ðri Þ is normalized by I w0 ðri Þ, the centroid velocity image, vc ðri Þ ¼ I w1 ðri Þ=I w0 ðri Þ ;
ð26Þ
is produced (cf. eq. [15]). A zero mean centroid velocity field may be produced by taking w1j ¼ vj hvi, where hvi is the intensity-weighted mean velocity. The centroid velocity image provides a summary of the kinematic properties of the emission, i.e., the spatial variation (in x and y) of the mean line-of-sight velocity, although some information is lost during the projection. Vectors involving higher order powers of the velocity coordinate, vj , can be used to generate images of the variance of velocities around the centroid velocity (the velocity dispersion), or the skewness or kurtosis of the line profiles. All the above basis vectors are user-defined, although obviously their form is specified such that a physical property of the medium is probed when the data is projected onto them. An alternative set of basis vectors that require no preassignment may be generated by the use of principal component analysis (HS97). In the formulation of HS97, the spectrum, or line profile, at each spatial grid point is taken to be the raw measurable quantity that will be subjected to PCA. From the ensemble of line profiles, the covariance matrix Sjk is calculated as Sjk ¼
835
n 1X Tij Tik : n i¼1
ð27Þ
A set of eigenvectors, umj , and eigenvalues, m , are determined from the solution of the eigenvalue equation for the
ð28Þ
The eigenvalue, m , equals the amount of variance projected onto its corresponding eigenvector, umj . The eigenvalues (and corresponding eigenvectors) are reordered from largest to smallest. In this implementation of PCA, P we do not explicitly subtract the mean value, Tj ¼ Tij =n from each channel image. The eigenvectors, umj , will generally have a nonzero amplitude in each channel, j, but the eigenvectors are orthogonal (HS97). To spatially isolate the sources of variance contained in the mth component, an eigenimage, I m ðri Þ, is constructed from the projected values of the data, Tij , onto the eigenvector, umj , p X I m ðri Þ ¼ Tij umj : ð29Þ j¼1
I m ðr
The symbols umj and i Þ refer to the mth eigenvector and eigenimage, respectively; the coupled eigenvector and eigenimage at order m is referred to as the mth principal component. The eigenvectors and corresponding eigenimages provide a sequence of diagnostic measures of the kinematic structure of the data. A variant of PCA, positive matrix factorization (Juvela, Lehtinen, & Paatero 1998), operates similarly but contains the additional requirement that the basis vectors be positive at all j. Other orthogonal vectors, such as Fourier components, wavelets, or user-designed vectors constructed by the Gram-Schmidt orthonormalization procedure, could be used to project the data onto images but these require input definitions. 5.2. Principal Components The translation from the intrinsic physical fields onto the PCA eigenvectors and eigenimages involves a substantial reordering of information. (The initial transformation of the physical fields onto the raw observational axes is already very complex.) In simple cases the eigenimages have relatively straightforward interpretations. However, when applied to the molecular ISM, the eigenimages are ‘‘ chaotic ’’ with no obvious symmetries, implying that the dominant observational structuring agent is turbulence. In light of this, HS97 adopted a purely statistical approach and characterized the eigenvectors and eigenimages by the e-folding scales of their autocorrelation functions in velocity space for the eigenvectors and in projected x, y space for the eigenimages. The measurements afforded by PCA differ from the more traditional approach of imposing a spatial scale over which, for example, a velocity dispersion is measured that can yield a ‘‘ continuum ’’ of measurements, limited only by the sampling rate. Instead, the recovered information is determined by the structure of the data and the user-interaction is limited to assigning a particular characterization scheme. 5.3. PCA as a Probe of Turbulence The velocity scales, vm , and spatial scales, Lm, obtained via PCA at each order m for a sample of molecular clouds were found to follow a power-law relation (HS97; Brunt & Heyer 2002b) v / L ;
ð30Þ
where the unsubscripted quantities refer to the ensemble
836
BRUNT ET AL.
(vm ; Lm ). The PCA statistical exponent, , is obtained from the set of v; L pairs that are larger than the respective spectral and pixel resolution limits. A bisector fit (Isobe et al. 1990) is always employed in the derivation of . In Brunt (1999) and Paper I a calibration of the measurable exponent to the intrinsic spectral slope of the velocity-field energy spectrum was established. The energy spectrum is the angular integral (over shells of constant wavenumber in the appropriate dimensionality) of the power spectrum of the velocity field. Energy spectra of the power-law form EðkÞ / k were considered for model fBm velocity fields that were subsequently ‘‘ observed,’’ and was measured via PCA. A more general calibration can be established here by considering the full range of real-space exponents q and their relationship to the measurable . Brunt (1999) established a relationship between and the coarse structure function exponent ¼ 2 ¼ q for all q (since the fBm velocity fields are nonintermittent), which is reproduced here as Figure 10. There appears to be some sensitivity of to within a group of equal , which is most noticeable at ¼ 1:5 2:5. The calibration has two regimes and is ð0:59 0:03Þ þ 0:32 0:01 for 0:67 ; ¼ ð31Þ ð1:07 0:08Þ þ 0:03 0:06 for > 0:67 : Clearly, 6¼ in these calibration equations—it is wrong of course to assume that this should be the case. The PCA method affords no analytic relationship between intrinsic and observed scaling exponents, in contrast to line centroid analysis (at least in the optically thin, uniform-density case) and thus must rely on empirical calibration. It would be interesting to examine size–line width relationships on the intrinsic and observed axes, as was done by BallesterosParedes & Mac Low (2003), but with a wider range of intrinsic to provide a similar calibration of that method.
Fig. 10.—Calibration of the PCA exponent to the intrinsic (nonintermittent) exponent for fBm velocity fields. Common symbols denote equal ; filled circles are ¼ 1, stars are ¼ 1:5, triangles are ¼ 2, open circles are ¼ 2:5, asterisks are ¼ 3, and open squares are ¼ 4. The two fitted lines are described in the text.
Vol. 595
In the broader context of the analysis of this paper, we may ask, for an intermittent field, at which moment q is an observed relationship ¼ ðq Þ most robust? 6. PRINCIPAL COMPONENT ANALYSIS
The simulated radiation fields are generated by a complex combination of velocity-field structure and density-field structure. For low-opacity observations, the velocity field is of paramount importance in generating the observed structures; as the level of opacity increases, the density field takes on a more prominent role. To disentangle the statistical properties of the velocity field alone promises to be very difficult. We apply PCA to the fields to evaluate how reliably this can be achieved. The observed spectroscopic data cubes were subjected to PCA with the intention of examining the sensitivity of the derived scaling exponents to each moment q of the original velocity field. Since each velocity field contains a variety of scaling exponents, we refer to each field by , which should be considered a label representing the full range of q. We measured the characteristic scales according to the prescriptions given in Paper I. To calculate the spatial autocorrelation functions, no padding was utilized for the simulated fields, as these are periodic in the spatial coordinates (the translation onto the observational axes breaks the line-of-sight periodicity). 6.1. PCA of the Simulated Radiation Fields First we discuss the results obtained for uniformly weighted 13CO observations. For these low-contrast TABLE 4 Measured PCA -Exponents from the Simulated Observations Density Field Uniform.............. Uniform.............. Uniform.............. MHD.................. MHD.................. MHD.................. Uniform.............. Uniform.............. Uniform.............. MHD.................. MHD.................. MHD.................. Uniform.............. Uniform.............. Uniform.............. MHD.................. MHD.................. MHD.................. Uniform.............. Uniform.............. Uniform.............. MHD.................. MHD.................. MHD.................. Uniform.............. Uniform.............. Uniform.............. MHD.................. MHD.................. MHD..................
Velocity Field Dephased Original Original Original Original Original Dephased Original Original Original Original Original Dephased Original Original Original Original Original Dephased Original Original Original Original Original Dephased Original Original Original Original Original
Tracer
13CO
1.0 1.0 1.0 1.0 1.0 1.0 1.5 1.5 1.5 1.5 1.5 1.5 2.0 2.0 2.0 2.0 2.0 2.0 2.5 2.5 2.5 2.5 2.5 2.5 3.0 3.0 3.0 3.0 3.0 3.0
(0.31) 0.51 0.05 0.43 0.05 0.38 0.04 0.45 0.04 0.48 0.09 0.52 0.01 0.64 0.03 0.61 0.02 0.53 0.02 0.50 0.03 0.59 0.07 0.58 0.02 0.76 0.02 0.73 0.01 0.64 0.02 0.66 0.04 0.75 0.02 0.70 0.02 0.83 0.02 0.80 0.02 0.74 0.03 0.81 0.05 0.80 0.02 0.83 0.02 0.87 0.02 0.82 0.02 0.82 0.03 0.87 0.03 0.85 0.03
13CO XCO XCO 13CO 12CO 13CO 13CO XCO XCO 13CO 12CO 13CO 13CO XCO XCO 13CO 12CO 13CO 13CO XCO XCO 13CO 12CO 13CO 13CO XCO XCO 13CO 12CO
No. 2, 2003
INTERSTELLAR TURBULENCE
837
Fig. 11.—Dependence of the PCA exponent on the intrinsic scaling exponents, 0.5, 1, 1.5, and 2. Open squares: 13CO observations using the dephased velocity field. Filled circles: 13CO observations using the original velocity field. Note that the q are different for the original and dephased fields except at q ¼ 2. The solid lines are the PCA calibration established in Paper I. The PCA exponent is more similar/consistent with q for q ¼ 0:5 1, showing that PCA is a low-order measure of velocity fluctuations. The small points are those from Fig. 10.
intensity fields, the spatial scales of the low-order principal components are poorly measured. The first-order principal component is a nearly uniform field, and its spatial scale is set simply by the overall size of the simulation. The first eigenimage ACF remains above 1/e at all lags. Consequently, we limited the analysis for all uniform density observations to principal components of second order and higher. For the ¼ 1 dephased field no principal components at orders higher than third are recovered. Fits of (eq. [30]) were made to v versus L for all the uniformly weighted 13CO fields, and Table 4 summarizes the measurements. For the ¼ 1 dephased field, we fitted the two recovered components; no error estimate is available. We emphasize that is used here only as a label (see Table 1). It is immediately evident that dephasing the velocity fields has a significant, measurable effect on the -exponents. Since and 2 are unchanged by the dephasing operation, we conclude that the effective order of PCA is not second order (q ¼ 2). The resulting exponents for all uniformly weighted 13CO fields were compared with the intrinsic exponents q
over a range of q. This enables one to identify the effective order at which PCA operates. In Figure 11 we plot the measured for all fields versus q for q ¼ 0:5, 1.0, 1.5, and 2.0. Clearly, is most directly related to q for the lower q values, i.e., PCA operates at a low order, of q 1 or less. At orders q > 2, the discrepancy of measured from the calibration line becomes even more pronounced. The dephased fields adhere closely to the calibration line, since these do not differ statistically from the initial calibration fields. For the original velocity fields there is typically a small increase in measured with respect to the calibration line, even at low q. This is of order 0.1 and is most noticeable for the lower fields. We applied PCA to the XCO simulated observations with the same approach used above. Again we find that the first principal component cannot be used, even in the densityweighted case. The measured -exponents have been added to Table 4. There is only a small effect arising from the influence of the inhomogeneous sampling effected by the density field (up to 0.1, comparable to the effect of densityweighting on the projected velocity line centroid power
838
BRUNT ET AL.
Vol. 595
Fig. 12.—Dependence of the PCA exponent on the intrinsic scaling exponents, 0.5, 1, 1.5, and 2. Open squares: XCO MHD density field. Filled circles: XCO uniform density field. The solid lines are the PCA calibration established in Paper I. The small points are those from Fig. 10. Inhomogeneous sampling by the density field does not significantly affect the PCA -exponent.
spectrum exponent; see Fig. 9). We plot the measured against each q in Figure 12. Again, the lower orders are generally preferred. Fits of (eq. [30]) to the v, L pairs for the densityweighted 13CO and 12CO observations were made, and the resulting exponents have been added to Table 4. Figure 13 shows the variation of with each q in comparison with the calibration line. The lower order exponents are generally more consistent again with the measurements. There is a small effect (up to +0.1 but typically less) arising from opacity effects; this is very much smaller than the effects of opacity on the projected velocity line centroid power spectrum index (Fig. 9). 6.2. The Effective Order of PCA Inspection of Table 4 will reveal that (disregarding the statistically different dephased fields) the -exponent varies only slightly (up to 0.1, which is a little larger than typical measurement errors) between fields observed under substantially different conditions. Given the wide variation in recovered line centroid power spectrum indices, this general
insensitivity to inhomogeneous sampling and saturation is encouraging. A more quantitative statement on this may be made by generating a table of predicted -exponents at each q. We did this only for the original velocity fields, since no distinction is made between q by the dephased velocity fields. Then we calculated, as a function of q, the rms difference between observed and predicted (via the calibration eq. [31]), D ¼ ð1=NÞ
N X
½i ðobservedÞ i ðpredictedÞ2 ;
ð32Þ
i¼1
from the N ¼ 5 measurements at each q. We also calculated D for the whole ensemble of measurements. These Dvalues are plotted versus q in Figure 14. We conclude that PCA operates at low order (q 0:5 1, or perhaps even less). However, it is worth noting that the minimum in D varies with the type of radiation field. There is no reason to discount the possibility that PCA operates at different effective order depending on the details of how the velocity field maps into the radiation-field structure. This is an interesting
No. 2, 2003
INTERSTELLAR TURBULENCE
839
Fig. 13.—Dependence of the PCA exponent on the intrinsic scaling exponents, 0.5, 1, 1.5, and 2. Open squares: 12CO. Filled circles: 13CO. The solid lines are the PCA calibration established in Paper I. The small points are those from Fig. 10. Opacity does not significantly affect the PCA -exponent.
possibility, but not one that can be acted upon with the small amount of available information. We thus accept an effective order of q 0:5 1 and note that details of the physical-to-observed transformation can inject variations up to 0.1 in the retrieved index . Finally, it can be seen that there is a tendency for the measured from the original velocity fields to be larger than can be accounted for even taking into account the operation of the PCA method at low order. In Figure 10 some scatter is evident in the fBm field results; it is not clear at present whether the MHD fields are consistent with this scatter or not (as effectively only one simulation forced onto different has been examined). This possible discrepancy deserves further attention using a wider range of numerical simulations. 6.3. Discussion In the preceding sections we have established calibrations between intrinsic low-order (q 0:5 1) scaling exponents in the three-dimensional velocity fields and the measurable exponents, , available from observations. In Brunt & Heyer (2002b), it was reported that the measured -exponents for molecular regions in the outer
Galaxy are consistent with an energy spectrum of 2 or more for the intrinsic velocity fields. If the velocity fields within these regions are strongly intermittent, then this is not necessarily the case. However, for current PCA measurements to be reconciled with a Kolmogorov spectrum,6 variations of q with q must be very large. All PCA measurements to date are summarized in Figure 15. Here we converted the measured from existing PCA measurements into 1, i.e., assuming a first-order operation of PCA. The 12CO measurements come from molecular clouds in the FCRAO Outer Galaxy Survey (Heyer et al. 1998; Brunt & Heyer 2002b) and additional fields examined in Brunt (2003b). These are from the Gem OB1 molecular complex (Carpenter, Snell, & Schloerb 1995) and from the Polaris Flare molecular cloud (Brunt 2003b). The 13CO measurements (Brunt 2003b) come from the Orion A (Bally et al. 1987), Orion B (Miesch & Bally 1994), Mon R2 (Miesch & Bally 1994), Perseus (Padoan et al. 1999), W3 (Heyer, Carpenter, & Ladd 1996), and NGC 7538 (Deane 2000) molecular clouds in the outer Galaxy, and several inner
6
The Kolmogorov spectrum has 3 ¼ 13.
840
BRUNT ET AL.
Vol. 595
Fig. 14.—rms differences between measured and predicted according to the fBm-based calibration established in Paper I. Each type of radiation field is identified by the symbol key at the top of the plot.
Galaxy molecular clouds (Simon et al. 2001). The two lowest 1 (12CO) occur in the Polaris Flare molecular cloud and the ‘‘ Gem OB1 western front ’’—both clouds that are devoid of star formation and may be in a state of decaying turbulence. The highest 1 (13CO) occurs in the W51 giant molecular cloud in the inner Galaxy. The steepness of this index may arise through line-of-sight blending, however (see Brunt 2003b). The majority of the 1 measurements lie in the range 0.5–0.8, with no obvious pattern otherwise. If a first-order operation of PCA is assumed, then there is a potentially interesting consistency between the range of inferred 1 and the range of 1 recently predicted by Boldyrev et al. (2002). Acceptance of the Boldyrev et al. (2002) theory as an explanation of the inferred range of 1 requires variations in the (multifractal) ‘‘ support ’’ of the turbulence between the clouds represented in Figure 15: 1 varies (for small variation in 2) depending on the dimension of structures on which the turbulence is dissipating. A possible test of this would be to see if larger values of correlate with more pronounced exponential tails in line profiles (see, e.g., Falgarone & Phillips 1990). Absent significant variation of q with q (e.g., 1 2 0:1 or more) the PCA results favor steep energy spectra. If 1 2 then the typical range of of 0.5–0.8 corresponds to in the range 2–2.6. Indices in this range are steeper than Kolmogorov ( 5=3) and Burgers turbulence ( 2). The preliminary measurements of from three-dimensional numerical simulations of turbulence (Ossenkopf & Mac Low 2002; Vestuto et al. 2003) are generally consistent with the range found observationally by the PCA method. If the level of intermittency is strong, however, then could be substantially lower than 2–2.6, as 1 only provides an upper limit to 2. For example, the Boldyrev et al. (2002) theory, if correct, suggests that almost our entire range of 1 could be consistent with 1:74 < < 1:89.
Fig. 15.—Histograms of the inferred first-order scaling exponents, 1, obtained from PCA. The heavy solid vertical lines delineate the theoretically predicted range of 1 (0.42–0.78), and the vertical dotted lines delineate the theoretically predicted range of 2 (0.37–0.445) from Boldyrev et al. (2002).
No. 2, 2003
INTERSTELLAR TURBULENCE
In principle, PCA operating at q 1 and line centroid analysis operating at q ¼ 2 could be combined to estimate two of the intrinsic velocity field exponents. However, the line centroid analysis would require high-sensitivity observations of an optically thin tracer, and there are still some uncertainties arising from the effect of density-velocity correlations. Both these uncertainties and the uncertainties of the PCA method will almost certainly preclude the application of this idea. Finally, we point out that is reliable to 0.1 under the conditions in which PCA has been tested to date—both with models and with multiple isotope data in single clouds. The inferred intrinsic exponents are less secure. The translation from to a q is only reliable within the conditions we have examined. 7. SUMMARY
We have analyzed the statistical properties of spectrally modified MHD simulations using differential scaling of velocity fluctuations at different orders to better describe the velocity fields. Intermittent velocity fluctuations are associated both with differential scaling of moments and with the appearance of broad, exponential tails in the velocity histograms. We have used simulated observations of numerical models of interstellar turbulence to verify and extend a previously established calibration of the PCA method that used simple molecular cloud models. PCA is a tracer (via the index ) of an intrinsic low-order scaling exponent q, with q 0:5 1, of the three-dimensional velocity field, as shown by analysis of simulated observations of modified MHD simulations. For nonintermittent fields, the distinction between moments (orders) is not important. We have outlined a prescription for translating back into the scaling properties of the intrinsic velocity field. The measured indices, , are predictable to approximately 0.1 of the previ-
841
ously established calibration, with a possible bias to steeper recovered indices (), even allowing for a low-order operation of the PCA method. The PCA method is not strongly affected by inhomogeneous sampling of the velocity field or by high opacities in the spectral lines. We have also analyzed the statistical properties of projected velocity line centroid fields, including the effects of density-weighting and opacity. We find that (1) in the case of unweighted projection of the mean velocity field, the measured power-spectrum index is the same as that of the three-dimensional field, as can be established analytically; (2) correlations between density and velocity result in shallower line centroid power spectra (i.e., with more smallscale power) in the density-weighted case; (3) it is the correlations between density and velocity rather than the amplitude of density fluctuations that are important in the production of shallower line centroid power spectra; (4) signatures of intermittent velocity fluctuations, as traced by the differential scaling of moments, are suppressed by the projection from three dimensions to two dimensions but are still discernible; and (5) significant line opacity also results in shallower line centroid power spectra, probably because of the velocity field being traced over the front surface of the cloud only, resulting, effectively, in a cut through the three-dimensional velocity field rather than an integration through it. The Dominion Radio Astrophysical Observatory is a National Facility operated by the National Research Council. The Canadian Galactic Plane Survey is a Canadian project with international partners and is supported by the Natural Sciences and Engineering Research Council (NSERC). C. B. is supported by a grant to the Canadian Galactic Plane Survey from NSERC. FCRAO is supported by NSF grant AST 01-00793. E. V.-S. and B. P. acknowledge financial support from Conacyt grant 27752-E.
REFERENCES Avila-Reese, V., & Va´zquez-Semadeni, E. 2001, ApJ, 553, 645 Lazarian, A., Pogosyan, D., Va´zquez-Semadeni, E., & Pichardo, B. 2001, Ballesteros-Paredes, J., & Mac Low, M.-M. 2002, ApJ, 570, 734 ApJ, 555, 130 Bally, J., Stark, A. A., Wilson, R. W., & Langer, W. D. 1987, ApJ, 312, Miesch, M. S., & Bally, J. 1994, ApJ, 429, 645 L45 Miville-Descheˆnes, M.-A., Levrier, F., & Falgarone, E´, 2003, ApJ, 593, 831 Bensch, F., Stutzki, J., & Ossenkopf, V. 2001, A&A, 366, 636 Murtagh, F., & Heck, A. 1986, Multivariate Data Analysis (Dordrecht: ˚ ., & Padoan, P. 2002, ApJ, 573, 678 Boldyrev, S., Nordlund, A Kluwer) Brunt, C. M. 1999, Ph.D. thesis, Univ. Massachusetts at Amherst O’Dell, C. R., & Castaneda, H. O. 1987, ApJ, 317, 686 ———. 2003a, ApJ, 583, 280 Ossenkopf, V. 2002, A&A, 391, 295 ———. 2003b, ApJ, 584, 293 Ossenkopf, V., & Mac Low, M.-M. 2002, A&A, 390, 307 Brunt, C. M., & Heyer, M. H. 2002a, ApJ, 566, 276 (Paper I) Ostriker, E. C., Stone, J. M., & Gammie, C. F. 2001, ApJ, 546, 980 ˚ . 1999, ApJ, ———. 2002b, ApJ, 566, 289 Padoan, P., Bally, J., Billawala, Y., Juvela, M., & Nordlund, A Burgers, J. M. 1974, The Nonlinear Diffusion Equation (Dordrecht: 525, 318 Reidel) Padoan, P., Goodman, A. A., & Juvela, M. 2003, ApJ, 588, 881 ˚ . 1998, ApJ, 504, 300 Carpenter, J. M., Snell, R. L., & Schloerb, F. P. 1995, ApJ, 445, 246 Padoan, P., Juvela, M., Bally, J., & Nordlund, A Deane, J. 2000, Ph.D. thesis, Univ. Hawaii Padoan, P., Rosolowsky, E. W., & Goodman, A. A. 2001, ApJ, 547, 862 ´ . Lis, D. C., Phillips, T. G., Pouquet, A., Porter, A., & Falgarone´, E Passot, T., Va´zquez-Semadeni, E., & Pouquet, A. 1995, ApJ, 455, 536 ´ . 2000, A&A, 356, 279 Woodward, P. R. 1994, ApJ, 436, 728 Pe´ty, J., & Falgarone, E Falgarone´, E´., & Phillips, T. G. 1990, ApJ, 359, 344 Pichardo, B., Va´zquez-Semadeni, E., Gazol, A., Passot, T., & Frisch, U. 1995, Turbulence (Cambridge: Cambridge Univ. Press) Ballesteros-Paredes, J. 2000, ApJ, 532, 353 Heyer, M. H., Brunt, C., Snell, R. L., Howe, J., Schloerb, F. P., & Rosolowsky, E. W., Goodman, A. A., Wilner, D. J., & Williams, J. P. 1999, Carpenter, J. C. 1998, ApJS, 115, 241 ApJ, 524, 887 Heyer, M. H., Carpenter, J., & Ladd, E. F. 1996, ApJ, 463, 630 Scalo, J. M. 1984, ApJ, 277, 556 Heyer, M. H., & Schloerb, F. P. 1997, ApJ, 475, 173 (HS97) Scoville, N. Z., & Solomon, P. M. 1974, ApJ, 187, L67 Isobe, T., Feigelson, E. D., Akritas, M. G., & Babu, G. J. 1990, ApJ, 364, Simon, R., Jackson, J. M., Clemens, D. P., Bania, T. M., & Heyer, M. H. 104 2001, ApJ, 551, 747 Juvela, M., Lehtinen, K., & Paatero, P. 1996, MNRAS, 280, 616 Stutzki, J., Bensch, F., Heithausen, A., Ossenkopf, V., & Zielinsky, M. Kitamura, Y., Sunada, K., Hayashi, M., & Hasegawa, T. 1993, ApJ, 413, 1998, A&A, 336, 697 221 Va´zquez-Semadeni, E. 2002, in ASP Conf. Ser. 276, Seeing Through the Kleiner, S. C., & Dickman, R. L. 1985, ApJ, 295, 466 Dust: The Detection of H i and the Exploration of the ISM in Galaxies, Kolmogorov, A. N. 1941, Dokl. Akad. Nauk S.S.R., 30, 301 ed. A. R. Taylor, T. L. Landecker, & A. G. Willis ( San Francisco: ASP), Kothes, R., & Kerton, C. R. 2002, A&A, 390, 337 155 Lazarian, A., & Pogosyan, D. 2000, ApJ, 537, 720 Vestuto, J. G., Ostriker, E. C., & Stone, J. M. 2003, ApJ, 590, 858
