Turbulence in the Ionized Gas of the Orion Nebula
Abstract
In order to study the nature, origin, and impact of turbulent velocity fluctuations in the ionized gas of the Orion Nebula, we apply a variety of statistical techniques to observed velocity cubes. The cubes are derived from high resolving power () longslit spectroscopy of optical emission lines that span a range of ionizations. From Velocity Channel Analysis (VCA), we find that the slope of the velocity power spectrum is consistent with predictions of Kolmogorov theory between scales of 8 and 22 arcsec (0.02 to 0.05 pc). The outer scale, which is the dominant scale of density fluctuations in the nebula, approximately coincides with the autocorrelation length of the velocity fluctuations that we determine from the second order velocity structure function. We propose that this is the principal driving scale of the turbulence, which originates in the autocorrelation length of dense cores in the Orion molecular filament. By combining analysis of the nonthermal line widths with the systematic trends of velocity centroid versus ionization, we find that the global champagne flow and smaller scale turbulence each contribute in equal measure to the total velocity dispersion, with respective rootmeansquare widths of 4–5 km s. The turbulence is subsonic and can account for only one half of the derived variance in ionized density, with the remaining variance provided by density gradients in photoevaporation flows from globules and filaments. Intercomparison with results from simulations implies that the ionized gas is confined to a thick shell and does not fill the interior of the nebula.
keywords:
hydrodynamics — \ionHii regions — ISM: individual objects: M 42 — ISM: kinematics and dynamics — turbulence1 Introduction
Detailed kinematic studies of \ionHii regions have shown that the velocity structure is extremely complex. The spectral lines are broadened in excess of thermal broadening, suggesting that the gas has disordered motion along the line of sight. After systematic motions such as expansion and rotation are accounted for, a random velocity component still remains. This has been attributed to turbulence in the photoionized gas. There have been numerous attempts to confirm and characterize this turbulence in both Galactic and extragalactic \ionHii regions based on statistical methods that examine the variation of the pointtopoint radial velocities with scale (for example, Münch, 1958; Roy & Joncas, 1985; O’Dell & Castañeda, 1987; MivilleDeschênes, Joncas, & Durand, 1995; Medina Tanco et al., 1997; Lagrois & Joncas, 2011, and references cited by these papers).
Structure functions of velocity centroids have become a standard statistical tool since von Hoerner (1951) investigated the projection of a threedimensional correlation function onto the plane of the sky. They measure the variation of the mean velocity integrated along the line of sight as a function of the planeofsky separation. However, the construction of structure functions over a wide range of scales requires high velocity resolution data with extensive spatial coverage. This can be realized by multiple highvelocityresolution longslit spectroscopic observations of an optically thin emission line at a variety of positions across an \ionHii region (O’Dell, Townsley, & Castañeda, 1987; Castañeda & O’Dell, 1987; Wen & O’Dell, 1993), or by FabryPerot interferometry (Roy & Joncas, 1985; MivilleDeschênes, Joncas, & Durand, 1995; Lagrois et al., 2011), which has much lower velocity resolution.
The observational requirements of high velocity resolution coupled with high spatial resolution and coverage mean that M42, the Orion Nebula, is ideally suited for this type of study. It is the closest and brightest \ionHii region ( pc, see e.g., O’Dell & Henney, 2008) and its geometry, main kinematic features and stellar population are well known (see O’Dell O’Dell, 2001 for a review). The inner part of the nebula measures some ( pc at the quoted distance) and is bright in optical emission lines such as [\ionSii] 6716,6731, [\ionNii] 6583, H and [\ionOiii] 5007. The principal ionizing star is Ori C (spectral type O7), which also possesses a fast stellar wind (O’Dell et al., 2009), and there are also 3 Btype stars. These 4 early type stars are the Trapezium cluster. M42 is a site of ongoing star formation on the nearside of the Orion Molecular Cloud (OMC1) and there is a large population of young stars, some of which are the sources of stellar jets and HerbigHaro objects and some have associated protoplanetary discs (proplyds, O’Dell, Wen, & Hu, 1993). M42 is an example of a blister \ionHii region (Zuckerman, 1973; Balick, Gammon, & Hjellming, 1974; TenorioTagle, 1979; O’Dell et al., 2009) and the photoionized gas is streaming away from the background molecular cloud with blueshifted velocities of order 10 km s.
Early attempts at characterizing the turbulence in M42 using the [\ionOiii] 5007 emission line were inconclusive and suffered from poor spatial coverage (von Hoerner, 1951; Münch, 1958). Detailed kinematic studies of the [\ionOiii] 5007 emission line (Castañeda, 1988) and the [\ionSiii] 6312 emission line (Wen & O’Dell, 1993) in the inner region of the Orion Nebula, which decomposed the line profiles into individual components, found evidence for the presence of turbulence. They obtained statistical correlations in the structure functions of the principal line component at spatial scales . However, the powerlaw indices do not agree with those expected from the standard Kolmogorov theory for isothermal, homogeneous turbulence nor with those predicted by von Hoerner (1951), taking into account projection smoothing. O’Dell & Wen (1992) carried out a similar study for [\ionOi] 6300. Although their structure functions appear to agree with Kolmogorov theory for a wide range of spatial scales, there is a problem with the line widths, which are much larger than they should be for the calculated structure function. O’Dell & Wen (1992) conclude that the [\ionOi] emission comes from the vicinity of the highly irregular ionization front, and that the large line widths are due to the acceleration of the gas away from the corrugated surface.
Other statistical techniques have been used to study turbulence in the interstellar medium. In particular, velocity channel analysis (Lazarian & Pogosyan, 2000) is optimally applied to spatially resolved PPV data from optically thin emission lines and can be used when the turbulence is supersonic, unlike the velocity centroid statistics, which are best employed when the turbulence is subsonic or at most mildly supersonic (Lazarian, 2004). It has been applied to the HI gas in the Galaxy (Lazarian & Pogosyan, 2000) and in the SMC (Stanimirović & Lazarian, 2001). There are no previous VCA studies of the ionized gas in the Orion Nebula (or any other ionized gas VCA studies of which we are aware).
In a previous paper (Medina et al., 2014, hereafter Paper I) we investigated the scale dependence of fluctuations inside a realistic model of an evolving turbulent \ionHii region. We calculated structure functions of velocity centroid maps for a variety of simulated emission lines but found that these were not a reliable way to recover the intrinsic 3D power spectrum of the velocity fluctuations of our numerical simulation. We found that the velocity channel analysis technique was a more promising approach, despite being intrinsically limited by thermal broadening. Using this method, we successfully recovered the logarithmic slope of the underlying velocity power spectrum to a precision of from simulated high resolution optical emissionline spectroscopy. Furthermore, our simulations revealed that multiple scales of energy injection due to champagne flows and the photoionization of clumps and filaments result in a flatter spectrum of fluctuations than would be expected from Kolmogorov theory of topdown turbulence driven at the largest scales.
In this paper, we apply both velocity channel analysis of positionvelocity arrays and structure functions of velocity centroid maps to high spatial and velocity resolution observations of low and highionization emission lines of the central region of the Orion Nebula. These observations were presented in the form of an atlas of emission lines by GarcíaDíaz & Henney (2007) and GarcíaDíaz et al. (2008). The longslit observations cover the full area at separation and allow maps to be made of the velocity moments. These maps are used to calculate the structure functions. In addition, we have the positionvelocity arrays of the individual longslit observations, with which we can perform the velocity channel analysis. We investigate whether the structure function method and the velocity channel analysis give the same information about the nature of the turbulence in the photoionized gas of the Orion Nebula. We also compare the results to our work in Paper I to assess to what extent numerical simulations can produce the detailed kinematics of \ionHii regions.
The paper is structured as follows: in Section 2 we describe the observational data and the statistical methods used in this work. Section 3 presents the results of the velocity channel analysis and the structure function calculations. In Section 4 we discuss our findings with respect to previous observational work and with regard to the numerical simulations of Paper I. We combine these results with complementary information derived from studies of linewidths and surface brightness fluctuations to form a picture of the causal links between physical processes and velocity and brightness fluctuations in the Orion Nebula. The conclusions are summarized in § 6.
2 Methods
2.1 Observational Data
The spectral line data was presented in the emission line atlas of GarcíaDíaz et al. (2008). In this paper, we use the H, [\ionSii] 6716, 6731 Å, [\ionNii] 6583 Å and [\ionOiii] 5005 Å data obtained with the echelle spectrograph attached to the 4 m telescope at Kitt Peak National Observatory (KPNO; for observational details see Doi, O’Dell, & Hartigan, 2004), with supplementary [\ionSii] data observed with the Manchester Echelle Spectrometer attached to the 2.1 m telescope at the Observatorio Astronómico Nacional at San Pedro Mártir (OANSPM), Mexico. The data covers the central (Huygens) region of the Orion Nebula^{1}^{1}1 The total observed area is from to in right ascension and from to in declination. The position of Ori C is R.A. Dec. (J2000) and the H, [\ionNii] and [\ionOiii] data consists of 96 approximately NorthSouth orientated slits at intervals, where the slit width corresponds to . The velocity resolution of these observations is 8 km s. The [\ionSii] 6716, 6731 Å dataset has a total of 92 NorthSouth pointings covering the same region, consisting of 37 positions observed at KPNO and 55 positions observed at OANSPM. The KPNO [\ionSii] data consists of two disjoint regions, one in the east and one in the west of the nebula, while the OANSPM data consists of 35 pointings with a slit width corresponding to and 12 km s velocity resolution, and 20 pointings with a slit width and 6 km s velocity resolution. The slit length at OANSPM is .
In the present paper, we use the calibrated positionvelocity (PV) arrays obtained by GarcíaDíaz et al. (2008) from each individual longslit spectrum for the velocity channel analysis of each emission line. For the secondorder structure functions, we work with the velocity moment maps, particularly the mean velocity map (i.e., map of velocity centroids) reconstructed from the PV arrays by GarcíaDíaz et al. (2008).
2.2 Statistical Methods
To analyze the velocity fluctuations due to turbulence in the ionized gas, we make use of two statistical tools: velocity channel analysis (VCA) of the PV data arrays and secondorder structure functions of the velocity centroid maps.
2.2.1 Velocity Channel Analysis
VCA consists of taking the spatial power spectrum of the emission intensity (or brightness) in velocity channels of spectroscopic PV data. First, the PV data is binned into velocity channels of width . The relative contribution of velocity fluctuations to fluctuations in the total intensity decreases as the width of the velocity slices increases, because thicker velocity slices average out the contributions of many velocity fluctuations. The very thickest velocity slice gives information only on the density spectral index. A velocity slice has width , where is the number of channels and for the spectra used in this work km s and km s.
The thickest velocity channel () corresponds to the total velocity integrated emissionline surface brightness at each pixel along the slit. The thinnest velocity slices for the present work are chosen to have km s, to give a good sampling of the spectrograph velocity resolution (FWHM) of 6–8 km s for all the emission lines. To use thinner velocity slices would be to amplify the noise (see Medina et al., 2014, hereafter Paper I).
The power spectra for each velocity slice of each PV array for a given emission line are summed and normalized by the total power. Since each spectroscopic observation corresponds to a finite slit length, we use a Welch window function to reduce the edge effects in the estimated power spectrum (Press et al., 1992) before applying the fast Fourier transform (FFT).
2.2.2 Secondorder Structure Functions
The secondorder structure function of the velocity centroids is (see Paper I)
(1) 
In this definition, is the twodimensional position vector in the plane of the sky, while is the separation vector. The normalization is by the number of pairs of points at each separation, , and the variance of centroid velocity fluctuations, , defined by
(2) 
Here, is the mean centroid velocity
(3) 
The summation in Equation 1 is over all data pairs for each separation, , while the summations in the centroid variation and mean (Eqs. 2 and 3) are over the total number of array elements, i.e., valid pixels in the plane.
We also define an intensity weighted structure function by
(4) 
where is the sum of the weights for each separation and , are the weights (i.e., intensities) of each pair of pixels. This form of the structure function obviously favours bright structures and reduces the contribution of fainter regions. This is one way to reduce the contribution of noise to the structure function.
The structure function is affected by smallscale, high velocity features such as jets and HerbigHaro objects and we aim to eliminate these pixels since they are not associated with the underlying turbulence. The first step is to examine the probability density function (pdf) of the velocity centroid maps and the corresponding velocity dispersion maps of the H 6583, [\ionNii] 6853, [\ionOiii] 5007, [\ionSii] 6716 and [\ionSii] 6731 lines from the emission line atlas published by GarcíaDíaz et al. (2008). In practice, a 2% threshold was uniformly applied to the velocity centroid pdf binned at 1 km s resolution. The same threshold was also applied to one side of the velocity dispersion pdf binned at 1 km s resolution to eliminate the small number of pixels with anomolously high values. Figure 1 shows the initial pdfs of the centroid velocities and the velocity dispersions of the 5 emission lines studied in this paper. Note that the velocity dispersions have not been corrected for the instrumental width (around 3 km s) nor thermal broadening. Thermal broadening most affects the H line, and is responsible for the peak of the velocity dispersion distribution occurring at km s, instead of the lower value of km s seen for the other emission lines. Figure 1 also highlights in grey the range of values used to create masks of “good” pixels for the subsequent structure function calculations. Pixels having or values outside these ranges are excluded from the structure function analysis.
Figures 2 and 3 show the 2D interpolated velocity centroid maps together with the masks constructed from pixels with extreme values of the centroid velocity or velocity dispersion identified through consideration of the pdfs. It can be seen that these pixels correspond to high velocity features, such as HH 201, 202, 203, 204, and 529, and also to numerical artefacts introduced during the original construction of the maps from the individual slits.
2.3 Application to the Orion Nebula
Although the standard theory assumes that the turbulent velocity fluctuations are homogeneous and isotropic, the real Orion Nebula is a complex, inhomogeneous object with illdefined boundaries. We expect variation in the powerlaw indices obtained from both the velocity channel analysis and the secondorder structure functions depending on position in the nebula due to largescale inhomogeneities. We estimate the confidence bounds on the powerlaw indices by resampling.
For the velocity channel analysis, this takes the form of a bootstrap Monte Carlo method, that is, resampling with replacement of the set of PV arrays. The powerlaw indices of the power spectra resulting from 10 different resamplings of the set of PV arrays (96 arrays in the case of the [\ionNii], H and [\ionOiii] emission lines and 20 arrays for [\ionSii] — see Fig. 17 for the positions of the highresolution [\ionSii] slits) are obtained using a leastsquares fitting procedure. The variation of these sample powerlaw indices provides an estimate of the confidence bounds for the powerlaw index across the Orion Nebula.
Furthermore, our longslit spectra are all oriented NorthSouth and so in order to check whether this has any effect on the powerlaw index, we analyze a supplementary PV dataset of [\ionOiii] 5007 observations from 18 slit positions perpendicular to the main data set (i.e., oriented EastWest). The slit positions and orientations are indicated in Figure 18 (O’Dell et al., 2015).
We estimate the effects of largescale spatial inhomogeneity on the secondorder structure function by evaluating (see Eqs. 1 and 4) for 100 distinct, randomly selected rectangular frames within each 2D map. Each frame has twothirds the and dimension of the initial map (see Fig. 4). The variation in powerlaw indices of the smaller maps give us an estimate of the spread of powerlaw behaviour due to location.
3 Results
We present the results of the velocity channel analysis in the form of power spectra plotted against wavenumber in units of number of waves per arcsecond and also number of waves per parsec. The conversion to parsec assumes a distance of 440 pc to the Orion Nebula (O’Dell & Henney, 2008). The secondorder structure functions are plotted against separation scale in both arcseconds and parsec.
3.1 Velocity Channel Analysis
The 1D, normalised, compensated power spectra for thin ( km s) and thick velocity channels are shown for all emission lines in Figures 5 and 6 and for the [\ionOiii] horizontal slits in Figure 7. The coloured points represent the average power spectrum of 96 distinct slits, in the case of the [\ionNii] 6583, H, and [\ionOiii] 5007 emission lines, where the slits cover uniformly the region shown in the velocity centroid maps. In the case of the two [\ionSii] emission lines, only the 20 highest resolution slits are used, since the other slits are too affected by noise at high wavenumber (see Fig. 17 for the positions of the 20 slits and GarcíaDíaz et al., 2008 for observational details). The 18 [\ionOiii] horizontal slit positions are shown in Figure 18 and the small number of observations in this dataset means that the compensated power spectra are much noisier than those of the other datasets.
By plotting it becomes very apparent that there is a break in the powerlaw behaviour around wavenumber arcsec ( arcsec for the [\ionSii] lines). This corresponds to a size scale of – arcsec. Indeed, the behaviour of the power spectrum for all the emission lines can be divided into four wavenumber ranges, which we have indicated by I, II, III and IV in the first panel of Figure 5. Regimes I and II correspond to wavenumbers smaller than the break point, while regimes III and IV correspond to wavenumbers larger than the break point.
For wavenumbers in regime II, the powerlaw indices of all the power spectra for all the emission lines are all for both the thin velocity channels and the thick channels, with the exception of the [\ionSii] emission lines, where the power law index is . In all cases the power laws for the thin channels are less steep than those of the thick channels. For wavenumbers in regime III, the powerlaw indices are all steeper than the critical value , indicating that there is very little power at small spatial scales. Noise dominates the spectra in regime IV. The spectral indices in regimes I and IV are similar and are both shallow indicating that the velocity fluctuations at the largest and smallest scales are uncorrelated.
The overall shape of the compensated power spectra are qualitatively similar for all lines and for both thin and thick slices. In particular, the break in the slope always occurs at the same scale: 7–8 arcsec, indicating that it must be a feature of the emissivity fluctuations in the nebula. The power spectrum of the thin slices is generally shallower (less negative) than that of the thick slices, which is indicative of the additional effect of velocity fluctuations (discussed in more depth in § 3.3).
The power spectra resulting from 10 different slit resamplings are shown in grey in Figures 5 and 6. The spread in powerlaw indices we obtain from these resamplings gives an estimate of the confidence bounds on the spectral index of the average power spectrum. This spread is not primarily due to observational uncertainties, which are negligible for all but the very smallest size scales (largest wavenumbers) but rather to largescale inhomogeneity of the power spectrum depending on the position of the observations.
3.2 Secondorder Structure Functions
3.2.1 Shape and Slope of Structure Function for each emission line
We use Equation 1 to calculate the secondorder structure functions of the velocity centroid maps for the [\ionNii] 6583, H, [\ionOiii] 5007, [\ionSii] 6716 and [\ionSii] 6731 emission lines. Figures 8 and 9 (top panels) show the complete structure function of the pdfselected pixels, together with the structure functions obtained from 100 smaller rectangular frames within the velocity centroid map (see Fig. 4). The fit to the powerlaw index of the structure function is performed over the spatial separations corresponding to regime II of the velocity channel analysis (roughly arcsec for the [\ionNii], H and [\ionOiii] emission lines, and arcsec for [\ionSii] emission lines; see Fig 5). The noise, estimated to be equal to the onepixel separation value , has been subtracted from the structure function. In the figures, we only plot the structure function up to spatial separations arcsec. Above this value, the structure function falls away steeply.
Unlike in the case of the power spectra, the structure functions show no clear break at a scale of arcsec. The powerlaw fit to scales arcsec is not too bad a fit to smaller scales also. The structure function has a slight negative curvature, giving gradually steeper slopes at smaller scales. The more pronounced steepening below arcsec is due to the spatial resolution of the observations and corresponds to the region where the seeing and interslit separation become important.
The steepest powerlaw indices over the fit range are for the H and [\ionOiii] emission lines. The form of the structure function curves is the same for both emission lines and the value of the powerlaw index is in both cases. The two [\ionSii] emission lines have very similar structure functions to each other. Their slopes are shallower than those of the H and [\ionOiii] structure functions, having . There is considerable variation at the smallest scales in the structure function among the different randomly selected frames. Finally, the structure functions for the [\ionNii] line can be divided into two distinct populations: one set corresponds to frames selected from the upper part of the velocity centroid map, the other to frames from the lower half (see Figure 4). One set has a powerlaw index similar to that of the [\ionSii] lines, while the other is intermediate between the [\ionSii] and H cases.
In the lower panels of Figures 8 and 9 we show alternative structure functions for the same emissionline velocity centroid maps. The solid red circles show the same pdfselected structure functions as in the upper panels, while the solid blue triangles show the structure functions obtained when all the pixels are used. There is little difference between the two structure functions for the [\ionNii], H and [\ionOiii] emission lines but the [\ionSii] structure functions show considerable differences. This is because more pixels are eliminated from the velocity centroid maps on the basis of the pdf analysis for the [\ionSii] lines: some of these pixels correspond to highvelocity features, such as HH objects, while others correspond to numerical noise. Evidently, these pixels influence the structure function.
Also in the lower panels of Figures 8 and 9, we plot the intensity weighted structure functions of the corresponding cases. The figures show that once again, for the H and [\ionOiii] emission lines, there is little difference between these new structure functions (open symbols) and the numberofpointsweighted structure functions (solid symbols). However, this time the [\ionNii] structure functions are different, with the intensity weighting producing structure functions with shallower slopes (smaller powerlaw indices). For the shorter wavelength [\ionSii] 6716 case, the intensity weighting makes little difference to the pdfselected structure function (red circles) but does reduce the noise contribution to the full structure function (blue triangles) and leads to a similar powerlaw index to the pdfselected cases over the fitted range. For the longer wavelength [\ionSii] 6731 case, the intensity weighting produces a shallower power law, even for the pdfselected case (red circles). It is not possible to plot the intensity weighted structure function for the full set of pixels on the same scale because the contributions of regions of bad pixels both in the velocity centroid map and the intensity map dominate at all scales.
3.3 Analysis of powerlaw indices
Description  Powerlaw  Relationship  Kolmogorov  

Index  Value  
3D emissivity fluctuations  
3D velocity fluctuations  (a)  
3D secondorder structure function  (a)  
2D secondorder structure function  
Projection Smoothing  (b)  
Sheetlike Emission  (c)  
Intensity fluctuations  
Very Thick Velocity Slice  (d)  
Thin Velocity Slice  
Shallow density  (d)  
Steep density  (d)  
(a) Kolmogorov (1941)  
(b) von Hoerner (1951)  
(c) Castañeda & O’Dell (1987)  
(d) Lazarian & Pogosyan (2000) 
Emission  Thin  Thick  Steep or  
Line  Shallow  
6716  Steep  …  
6731  Steep  …  
6583  Shallow  …  
H 6563  Shallow  …  …  …  
5007  Shallow  …  
5007H  Shallow  … 
Emission  Thin  Thick  Steep or  

Line  Shallow  
6716  Steep  NEG  N/A  
6731  Steep  NEG  N/A  
6583  Steep  NEG  N/A  
H 6563  Steep  NEG  N/A  
5007  Steep  NEG  N/A  
5007H  Steep  NEG  N/A 
Emission  
Line  
6716  …  
6731  …  
6583  …  
H 6563  …  
5007  … 
We want to use the powerlaw indices from the velocity channel analysis and the secondorder structure functions obtained from the spectroscopic observations to recover the threedimensional velocity statistics of the ionized gas in the Orion Nebula. There are various relationships between the observationally derived powerlaw indices and the underlying 3D density and velocity fluctuations, which are summarized in Table 1 and described below.
The velocity channel analysis gives us the powerlaw indices of the average power spectra of thin and thick velocity slices. The thickest velocity slice corresponds to the velocityintegrated surface brightness and the powerlaw indices of the spectra of the thick slices are predicted to be equal to the those of the 3D power spectra of the respective emissivities, i.e., (Lazarian & Pogosyan, 2000). Thus, we can use our calculated value of as an estimate of . The critical value divides ‘steep’ from ‘shallow’ power spectra. A ‘steep’ emissivity spectrum () means that the emissivity is dominated by fluctuations at intermediate to large scales, while a ‘shallow’ spectrum () indicates that the emissivity is dominated by fluctuations at small scales. As the velocity slices become thinner, velocity fluctuations dominate the spectra. By using a combination of thick slice and thin slice spectral indices, we can, in theory recover the velocity powerlaw index.
There are relations between the spectral index of the 3D velocity fluctuation power spectrum and the powerlaw indices of the VCA power spectra, depending on whether the power spectrum of the emissivity fluctuations is ‘steep’ or ‘shallow’ (Lazarian & Pogosyan, 2000). In the steep case, the powerlaw index of the average power spectrum of the thin isovelocity channels is given by
(5) 
where is the spectral index of the 3D velocity power spectrum (see Table 1). For the ‘shallow’ case, the relation is with the difference in the powerlaw indices of the average powerspectra of the thin and thick velocity slices:
(6) 
There are also relationships between the powerlaw indices of the twodimensional structure function of the emissionline velocity centroid map and the threedimensional structure function of the velocity field. For projection from 3 to 2 dimensions we have , which is known as projection smoothing (von Hoerner, 1951; Münch, 1958; O’Dell & Castañeda, 1987; Brunt & Mac Low, 2004). If, however, the distribution of emitters is sheetlike, i.e., essentially twodimensional, then (Castañeda & O’Dell, 1987; MivilleDeschênes, Joncas, & Durand, 1995). The threedimensional structure function powerlaw index is related to the threedimensional power spectrum spectral index () of the underlying velocity field through , where for homogeneous, incompressible turbulence (i.e., Kolmogorov spectrum) (see Table 1).
3.3.1 VCA: Regime II
For the low wavenumber range, arcsec, corresponding to spatial scales arcsec, the powerlaw indices of the thick slice power spectra for regime II are and are listed in Table 2. We consider each emission line in turn, beginning with the [\ionSii] lines, which represent the lowest ionization state of our data. and for the [\ionSii] 6716 and [\ionSii] 6731 lines, respectively. This is the critical value separating ‘steep’ from ‘shallow’ power spectra. Using the relations between powerlaw indices listed in Table 1 (see also Eq. 5), we obtain and for the short and long wavelength [\ionSii] lines, respectively. The corresponding powerlaw indices of the 3D velocity fluctuations are then and , respectively. These values of are consistent with the Kolmogorov value of . Note that the [\ionSii] lines are least affected by thermal broadening and, moreover, the 20 spectra used to calculate the power spectra have the highest velocity resolution of all our observational data.
Moving to higher ionization lines, the calculated ‘thick’ and ‘thin’ slice powerlaw indices of the [\ionNii] 6583 spectra fall into the ‘shallow’ regime. Following the same analysis as above, we recover , which is again close to the Kolmogorov value. For the H power spectra, there is little difference between the thin and thick powerlaw indices since thermal broadening is important for this line (see, e.g., Paper I). It is therefore not possible to recover a meaningful value of . Finally, the [\ionOiii] 5007 power spectra again fall into the ‘shallow’ regime and we obtain .
It is gratifying that the [\ionNii], [\ionOiii] and [\ionSii] power spectra powerlaw indices are consistent with each other and with the Kolmogorov value, within the errors.
3.3.2 VCA: Regime III
For the high wavenumber range arcsec, corresponding to spatial scales arcsec, only the [\ionSii] 6716 and [\ionSii] 6731 power spectra have ‘thick’ slice powerlaw indices steeper than their ‘thin’ slice powerlaw indices. The [\ionNii], H and [\ionOiii] have curiously opposite behaviour and their powerlaw indices cannot be accommodated by the theory (see Table 3). For the [\ionSii] power spectra, the value of falls into the ‘steep’ regime. However, for both [\ionSii] lines, the value of is less than and so this would lead to negative values for on applying Equation 5.
3.3.3 Structure function power indices
We fit the structure functions over the separation range corresponding to VCA regime II. The structure function graphs themselves do not strongly suggest that any particular range is preferable, although in the range arcsec the structure function falls off very steeply due to the effects of seeing and the interslit separation, which has to be interpolated to produce the velocity centroid maps. At the largest separations, the structure function powerlaw index flattens off, i.e for arcsec. Over the separation range we fit, we find for the H and [\ionOiii] lines, and for the [\ionSii] and [\ionNii] lines.
The values we obtain for allow us to determine a range for . The 3D structurefunction index must lie between the two limits given by projection smoothing and a sheetlike distribution of emitters, that is (see Table 1). Our lowest value is and the highest is , therefore must lie in the range . This is consistent with the results of the VCA, for which we found but is such a wide range of values that it is not a useful diagnostic.
3.4 Lineofsight velocity dispersions
The spectral lines are broadened beyond the thermal component by a combination of ordered and disordered motions along the line of sight. Moreover, the lines often do not have Gaussian profiles and can sometimes be separated into distinct components that could originate from different physical components of the ionized gas in the nebula having different kinematic characteristics (Castañeda, 1988). Contributions to ordered motions include expansion, contraction and rotation, while random motions are often ascribed loosely to “turbulence”.
The nonthermal lineofsight velocity dispersion may be estimated by correcting the observed line widths for the contributions from the spectrograph resolution, finestructure splitting, and thermal Doppler broadening, as in Eq. (2) of GarcíaDíaz et al. (2008). A striking fact about the nebular line profiles is that the lineofsight velocity dispersion is several times larger than the planeofsky dispersion in mean velocities. This is illustrated in Figure 10, which shows fluxweighted histograms of the nonthermal RMS line widths versus mean velocity for all lines of the GarcíaDíaz et al. (2008) spectral atlas. The lineofsight RMS velocity dispersion is 9–10 , whereas the RMS planeofsky dispersion of the mean velocities is only 2–4 .
However, the observed line widths are affected by additional broadening due to dustscattering (Castañeda, 1988; Henney, 1998), which gives an extended red wing to the line profile that contains 10–20% of the total line flux (see Fig. 11). By means of fitting multiple Gaussian components to each line profile, it is possible to effectively remove this scattered component and calculate statistics for the line core alone. We have done this for the [\ionOiii] line, with results shown in the inset box of Figure 11. It can be seen that the lineofsight velocity dispersion is reduced by almost a factor of two.
It is more difficult to perform the Gaussian decomposition for the lower ionization lines since the backscattered component is not as cleanly separated from the core component. On the other hand, an alternative method of suppressing the effects of scattering on the line widths is to use the Full Width at Half Maximum (FWHM) instead of the RMS width. Even a weak component with a velocity that is very different from the line centroid can have a large effect on the due to the dependence of the second velocity moment. But the same component will have almost no effect on the FWHM if its amplitude is less than half that of the line core. We therefore define a FWHMbased effective as , with the constant of proportionality being chosen so that for a Gaussian profile. The results are shown in Figure 12, where it can be seen that nearly all the lines show substantial reductions in the lineofsight line widths with respect to the momentderived values of Figure 10. The exception is , where larger thermal broadening means that is still contaminated by the backscattered component due to blending. Note that is insensitive to not only the scattered component, but also to any other weak component that is not blended with the line core. For the lowionization lines such as [\ionNii] and [\ionSii], this includes the kinematic component known as the Diffuse Blue Layer (Deharveng, 1973), which produces line splitting in the SE and N regions of our observed maps.
(a)  (b) 

4 Discussion
4.1 Comparison with previous structure function determinations
Figure 14 illustrates how various effects modify the purely turbulent structure functions for a highly idealized case. The most promising range of length scales for measuring the turbulent velocity spectrum is between a few times the seeing width and about half the correlation length, , of the turbulence. At scales larger than , the structure function flattens as it tends towards the asymptotic value of 2 for a homogeneous random field. If there is a linear velocity gradient across the map, then the structure function will steepen again at the largest scales. Alternatively, if the turbulent velocity dispersion is inhomogeneous, being larger in the center of the map than in the periphery, then the structure function slope will become negative at the largest scales (not illustrated). The figure does not include the effects of noise, but that is easily dealt with in the case that the noise is “white” (spatially uncorrelated, such as shot noise), since the effect is to simply add a constant value to the structure function at all scales. This can be estimated either from the variance of velocity differences at separations significantly less than the seeing width or from independent observations of the same spatial point, and then subtracted from the observed structure function prior to analysis. Spatially correlated “noise” due to systematic instrumental effects is more difficult to deal with, since it will add a spurious scaledependent term to the structure function.
Range  

Reference  Ion  Method  ()  ()  Slope 
O’Dell & Wen 1992  [\ionOi]  Mean  3:  6–85  : 
This paper  [\ionSii]  Mean  5.4  7–32  
This paper  [\ionNii]  Mean  5.6  8–22  
Wen & O’Dell 1993  [\ionSiii]  Comp A  13.8  5–20  : 
This paper  Mean  9.4  8–22  
This paper  [\ionOiii]  Mean  10.2  8–22  
Castañeda 1988  [\ionOiii]  Comp A  13.7  3–15  
This paper  [\ionOiii]  Comp A  15.3  3–15 
Previous studies of the velocity structure function in Orion have been carried out based on slit spectra (Castañeda, 1988; O’Dell & Wen, 1992; Wen & O’Dell, 1993). Table 5 compares these results with our own for different emission lines, ordered from lower to higher ionization. In spite of the differences in methodology, a broad agreement is seen, with both the magnitude of the velocity dispersion and the steepness of the structure function slope increasing with ionization. The most directly comparable methodology to our own is that of O’Dell & Wen 1992 who used the fluxweighted mean velocity of the [\ionOi] 6300 line. The fitted range of to extends to larger scales than in studies of other lines, and this biases the slope determination towards lower values due to the slight curvature of the structure function.
The other studies are harder to compare with our own since they are based on multicomponent Gaussian fits to the line profiles, such as shown in Figure 11. In order to check the effects of this methodological difference, we have calculated the [\ionOiii] structure function for the strongest component of the threeGaussian decomposition of our line profiles (component A, see § 3.4), with results that are also included in Table 5. We find a slightly larger planeofsky velocity dispersion and shallower structure function than we found using the mean velocity of the entire profile, which is consistent with the results of Castañeda (1988) for the same line. Experimentation shows that this is partly due to underdetermination of the Gaussian fits, particularly for components A and B, which are severely blended at most positions. This results in a fitting degeneracy between the velocity separation and the flux ratio of the two components, which spuriously contributes to the variation in .
The structure function slopes obtained by Mc Leod et al. (2016), based on integral field spectroscopy with the MUSE instrument (Weilbacher et al., 2015) are significantly flatter than all other studies, and we have decided not to include them in the comparison table. For example, they obtain a slope of for [\ionOiii] and for [\ionO1]. This appears to be the result of uncorrected fixedpattern noise in their mean velocity maps, which can be seen as tartanlike horizontal and vertical stripes in their Fig. 10. We will show in a following paper that once these instrumental artifacts have been removed, the MUSE results are consistent with other studies.
4.2 Comparison with previous simulation results
In our previous paper (\al@{2014MNRAS.445.1797M}), we used velocity channel analysis and secondorder structure functions to analyse our 3D numerical simulations of an evolving \ionHii region in a turbulent molecular cloud. The spectral index of the 3D power spectrum of the underlying velocity fluctuations for our simulations was , which is shallower than that predicted by Kolmogorov theory. We found that for our velocity channel analysis of the emission lines of the heavier ions, [\ionOiii], [\ionNii] and [\ionSii], the measured values of the thin velocity slice spectral index agree to within with the values derived using the formulae highlighted in Table 1 using the simulation value of . On the other hand, the secondorder structure functions we calculated from our simulations proved less capable of reliably recovering the 3D velocity power spectrum. These results motivated the work in the current paper, that is, by applying velocity channel analysis to high resolution spectroscopic observations we have the possibility of recovering the spectral index of the underlying 3D velocity fluctuations.
In the present work, we do not have access to the underlying 3D velocity distribution. Our VCA spectral indices show some similarities but also some differences with respect to those obtained from the simulations. We will discuss only the heavier ions, since the H is seriously affected by thermal broadening. To begin with, for the simulations, we found that the [\ionNii] and [\ionSii] values of (thick velocity slice spectral index) could be considered “shallow” (), while the [\ionOiii] had a “steep” spectral index (). For the compensated power spectra , a powerlaw index equal to the critical value would appear horizontal. We remark that the 2D power spectra of the simulations are not at all noisy, unlike the 1D power spectra presented in Section 3.1 and, moreover, a single power law is sufficient to cover most of the wavenumber range. In contrast, our observationally derived power spectra can be split into different regimes where different power laws apply with sharp breaks between them (see Figs. 5 and 6). In the principal wavenumber range (regime II in Section 3.1), the spectral index of the thick velocity slices, , is “shallow” for [\ionNii] and [\ionOiii] but apparently steep for [\ionSii]. Thus, while there is agreement in the [\ionNii] spectral index between observational and simulation results, the [\ionOiii] spectral index is shallower in real life. The spectral index of [\ionSii] is misleading, since this emission line is affected by collisional deexcitation for electron density greater than about cm, which will suppress power at higher (smaller length scales). This was not an issue in the numerical simulations because the electron densities were lower than in the Orion Nebula.
With regard to the powerlaw indices of the secondorder structure functions, our simulations found a clear sequence . The observational results, summarized in Table 4, show a similar sequence but the absolute values are slightly different. We note that the offset of approximately 0.3 in values between simulation and observational results is no larger than the spread in observationally derived results obtained with different methodologies discussed in Section 4.1. The simulation structure function follows a clear power law over a wide range of spatial separations (straight line in loglog space) whereas the observational structure function has a slowly varying power law as a function of separation scale.
Both the structure function powerlaw indices and the VCA values are potential diagnostics of the underlying velocity power spectrum. The differences we find in these quantities between the simulations and the observations are consistent with a steeper velocity spectrum for the real nebula. On the other hand, the thick slice VCA spectral index , which is simply the spectral index of the surface brightness power spectrum, is a diagnostic of the underlying emissivity power spectrum. The simulation results give very different values for the [\ionNii] and [\ionOiii] thickslice spectral indices, suggesting a very different spatial distribution for these two emission regions. The observational values are more similar for [\ionNii] and [\ionOiii].
We remark here that in Paper I, we used the projected emission maps to estimate the mean radius of the nebula in a given emission line and used this as the upper limit of the separation range for the powerlaw fit to the structure functions. The lower limit was set to be 10 computational cells, i.e. the scale at which numerical dissipation effects cease to be important. The corresponding wavenumber ranges were used in the VCA fits. This is in contrast to the present paper, where the VCA graphs clearly indicate the wavenumber ranges of interest and the corresponding separation scales are used to fit the structurefunction power laws.
How can we explain the differences between the observational results and those obtained from the numerical simulations? To begin with, the simulations are not a model for the Orion Nebula: the size of the simulation box was 4 pc, which, although similar in size to the Extended Orion Nebula (O’Dell & Harris, 2010), is an order of magnitude larger than the observed region studied in this paper. The simulation ran for 300,000 yrs, while the Orion Nebula is of uncertain age ( million years). The simulations omitted several physical processes. The most important one is probably the stellar wind from Ori C (O’Dell et al., 2009), which will evacuate a cavity inside the distribution of photoionized gas such that even the high ionization emission lines will come from a thick shell rather than a filled volume. This can explain why the observed [\ionOiii] surface brightness spectral index is shallower than in the simulations. The wind bubble will also act as an obstacle around which the photoionized gas streams away from the ionization front in a champagne flow (see Arthur & Hoare, 2006 for examples of numerical simulations of this scenario). Furthermore, the limited dynamical range of the simulations may have prevented the development of a fully resolved turbulent energy cascade.
Finally, the real Orion Nebula contains a large population of lowmass stars, some of which are sources of jets and HerbigHaro objects. Although the effect of such structures is spatially limited and we have excised extreme velocity features from our sample used to calculate the structure function (see Section 2.2.2), there still could be a small contribution from bright features that lie moreorless in the plane of the sky.
4.3 Turbulent contribution to spectral line broadening
So far we have concentrated on the slopes of the structure function and power spectra, since they can be theoretically related to the spectrum of underlying velocity fluctuations in the nebula. However, other questions require consideration of the magnitude of the velocity fluctuations, as measured by different techniques. Such questions include whether turbulence alone is sufficient to account for the nonthermal line broadening observed in the nebula, and whether that turbulence is driven primarily by the large scale thermal expansion of the nebula, or by smaller scale photoevaporation flows, and to what extent stellar winds play a role.
The planeofsky dispersion in centroid velocities is the RMS width of the marginal distribution along the axis of Figure 10, which is the same as the quantity used to normalize the structure functions in § 2. Figure 13(a) summarizes the observational determination of and , while Figure 13(b) shows the same quantities calculated for the turbulent \ionHii region simulation of Paper I, as presented in Appendix D. As discussed above, the FWHMderived and reflectioncorrected observational values are the most reliable, and these fall broadly in the same region of the graph as do the evolutionary tracks, with and .
The fact that the lineofsight velocity dispersion is roughly twice the planeofsky velocity dispersion can be interpreted in at least two different ways. In the case of a homogeneous turbulent velocity field with characteristic correlation length , the projection from three to two dimensions over a lineofsight depth reduces the planeofsky amplitude of fluctuations if , a phenomenon known as “projection smearing” (von Hoerner, 1951; Scalo, 1984). Our and correspond to and in Fig. 1 of Scalo (1984), from where it can be seen that a value of requires –, depending on the steepness of the velocity fluctuation spectrum. The results from our structure function analysis (§ 3.2) imply a correlation length – pc for all lines, which would require a very large lineofsight depth pc in order to explain the observed by projection smoothing. This is inconsistent with independent evidence (Baldwin et al., 1991; O’Dell, 2001; GarcíaDíaz & Henney, 2007) that the emitting layer thickness is much smaller than this in the region covered by our maps: – pc, being thinner on the West side and for the lower ionization lines.
The same discrepancy arises when the projection smearing argument is applied to our simulated \ionHii region. In this case, – pc for evolutionary times later than Myr (Fig. A4 of Paper I ), thus requiring pc, which is clearly unacceptable since it is bigger than the computational box of the simulation. For both the simulation and the observations, it is clear that for the highionization lines and for the lowionization lines. Therefore, projection smearing of the largescale fluctuations is negligible and cannot explain the difference between and .
A second, contrasting interpretation of the evidence would be in terms of largescale, ordered motions. Consider an emission shell that expands at velocity . If the shell emissivity is homogeneous, then the integrated line profile is rectangular, with mean velocity and velocity width . Furthermore, the spatially resolved line profile at any point will also have , so that . However, if there are emissivity fluctuations between different parts of the shell, then will fluctuate on the plane of the sky, according to the relative brightness of the redshifted and blueshifted hemispheres. The required RMS fractional variation in the emissivity on the scale of the shell diameter is found to be .
The observed largescale () brightness fluctuations are illustrated in Figure 15, which shows lognormal fits to the PDFs of surface brightness, , after normalizing by the mean, , and binning the maps at pixels. The RMS width of the lognormal PDF, is seen to be in the range 0.45–0.6 for all lines. This is related to the RMS fractional brightness fluctuation as , or if . The relationship between and depends on both lineofsight projection (Brunt, Federrath, & Price, 2010), which tends to make , and fluctuations in the foreground dust extinction, which have the opposite effect of increasing . The first effect dominates, so that the surface brightness PDFs imply (see § 4.5 for details). This is larger than the value derived in the previous paragraph, which implies that emissivity fluctuations combined with an ordered velocity field are entirely sufficient to explain the observed planeofsky variation in mean velocities, without requiring any fluctuations in the velocity field itself.
Although it is a priori unlikely that there are no velocity fluctuations in the ionized gas, this is yet another reason why the structure function of the mean velocity is not an effective diagnostic of these fluctuations in the presence of strongly inhomogeneous emissivity and largescale velocity gradients. In Orion, the ordered large scale expansion of the nebula is an asymmetrical champagne flow away from the background molecular cloud (Zuckerman, 1973), which produces systematically larger blueshifts with increasing ionization (e.g., Fig. 11 of Baldwin et al., 2000). This offers a simple method for estimating the relative contribution of ordered versus turbulent motions to the total velocity dispersion. The mean systematic difference between the [\ionOi] and [\ionOiii] centroid velocities is (Table 2 of GarcíaDíaz et al., 2008), which gives a champagneflow contribution to the velocity dispersion of . The turbulent contribution to the velocity dispersion is then . The uncertainties in this analysis are large, so that all that can be confidently asserted is that the ordered and turbulent velocity dispersions are roughly equal with –.
4.4 What is the significance of the 22 arcsec and 8 arcsec length scales?
The velocity channel analysis suggests that two length scales are important in the Orion Nebula, corresponding to the limits of regime II (see Fig. 5). The break in power law at 22 arcsec and 8 arcsec occurs for both thin and thick velocity slices. This indicates that it is a feature of the emissivity power spectrum, and not the velocity power spectrum. Below 8 arcsec (regime III), the emissivity power spectrum is very steep in all the lines, indicating that smallscale fluctuations are relatively unimportant. Above 22 arcsec (regime I), the power spectrum is very flat, similar to the noisedominated spectrum in regime IV, suggesting that fluctuations are relatively uncorrelated on larger scales. This is underlined by the analysis in Appendix E, where it is shown that a combination of Gaussian brightness peaks with widths (FWHM) from and can capture the broad features of the observed power spectra.^{2}^{2}2The factor of two difference between these widths and the length scales discussed earlier in the paragraph is because a fluctuation consists of both a peak and a trough, and so the width of the peak is only one half of the fluctuation wavelength.
The outer scale of coincides with scale where the structure functions reach a value of unity (Fig. 8), which corresponds to the correlation length, , of the velocity fluctuations (see Fig. 14). It is therefore plausible to associate this scale, which corresponds to a physical size of parsec, with the driving scale of turbulence in the nebula. The inner scale of , corresponding to a physical size of parsec, is harder to associate with any particular process since the structure functions (Fig. 8) show no apparent feature at this scale.
4.5 Does velocity turbulence cause the surface brightness fluctuations?
The surface brightness fluctuations on the plane of the sky are primarily caused by emissivity fluctuations within the nebular volume, which are in turn caused by fluctuations in electron density, temperature, and ionization. The temperature and ionization dependence of the emissivity is very different for each line, but the electron density dependence is similar in all cases, being in the low density limit, which is appropriate for all but the [\ionSii] lines. It therefore seems likely that any commonalities in the statistics between all the different emission lines will give us information about the electron density fluctuations within the nebula.
In § 4.3 it was shown that the rms fractional surface brightness variation in 2D is for all lines, and the rms emissivity variation in 3D is predicted to be , where the “deprojection factor” is – (Brunt, Federrath, & Price, 2010). On the other hand, if the emissivity fluctuations are due to variations in the density squared, then the rms fractional 3D density variation is , which approximately cancels out the deprojection factor so that . If the density fluctuations are caused by the turbulent velocity fluctuations, then numerical simulations (Konstandin et al., 2012) show that there is a linear relationship between and the rms Mach number, , of the turbulence: , where to , depending on whether the turbulent driving is primarily solenoidal or compressive. The rms Mach number is the ratio of the velocity dispersion to the ionized isothermal sound speed , where (see § 3.4) and . Thus, so that, given , an upper limit to the turbulent contribution to the density fluctuations is .
Furthermore, the slopes of the surface brightness fluctuation spectra in regime II are significantly shallower than expected from a topdown turbulent cascade in the subsonic limit (Konstandin et al., 2015). A similar result is seen in our numerical simulations (see § 4.2), with the important difference that in the simulation the velocity spectrum is also shallow, whereas the thinslice VCA analysis for the observations (§ 3.3.1) is consistent with a Kolmogorov slope for the velocity fluctuation spectrum in regime II.
We therefore require a further mechanism to explain the roughly 50% of the variance in ionized density that cannot be accounted for by turbulent velocity fluctuations. This could plausibly be provided by the brightrimmed structure of the photoevaporation flows away from dense molecular globules and filaments (e.g., (Bertoldi & McKee, 1990; Henney et al., 2009)), which are responsible for driving the turbulence. We have calculated the emissivityweighted density PDF for a simple model of a single spherically divergent, isothermal evaporation flow from a Dcritical ionization front (Dyson, 1968) and find . For an an ensemble of such flows with varying peak densities the would be even higher, so that in order for their global contribution to rival that of the velocity fluctuations it is sufficient that a fraction – of the total emission should come from such flows.
5 Speculation
We offer a speculative account of the complex web of physical processes that give rise to to the velocity and brightness fluctuations that we observe in the Orion Nebula. This is illustrated in Figure 16, where the most important causal links are shown by thick arrows and secondary processes by thin arrows.
The principal origin of all structure in the \ionHii region is the highly filamentary and clumpy density structure in the molecular cloud from which it is emerging, which in turn has its origin in some combination of thermal and gravitational instability and supersonic turbulence (Padoan & Nordlund, 2002; BallesterosParedes et al., 2011). In the molecular gas, thermal pressure is negligible compared with magnetic pressure, turbulent ram pressure, and the gravitational potential. However, the large temperature increase that accompanies photoionization means that thermal pressure dominates in the \ionHii region, so that density gradients are converted into pressure gradients that can accelerate the gas. The fractal nature of the molecular density means that gas acceleration occurs on multiple scales, from the global outward radial expansion of the \ionHii region (which in Orion is a highly onesided champagne flow) down to photoevaporation flows from individual globules. One piece of evidence for a direct connection between molecular density fluctuations and ionized velocity fluctuations is that Kainulainen et al. (2016) find correlation lengths of order pc for the separations of molecular cores along the ridge that lies behind the Orion Nebula, which is similar to the correlation lengths we find for the velocity fluctuations in the nebula.
Ionized density fluctuations can arise directly from the molecular density fluctuations, such as the bright rims at the edges of photoionized globules (Henney et al., 2009), and this is most important in the lower ionization zones near the ionization front where the [\ionSii] and [\ionNii] emission is strong. In the more highly ionized interior of the nebula, it is collisions between opposing velocity streams that produce the ionized density fluctuations, but these fluctuations are less extreme than those seen in molecular gas because the turbulence is subsonic.
The ionized density fluctuations are the primary determinant of the emission line surface brightness fluctuations (§ 4.5), although ionization and temperature structure can make a contribution for particular lines and there is also a direct contribution from foreground molecular density fluctuations via dust extinction (O’Dell & YusefZadeh, 2000).
Finally, a variety of other processes, such as O star winds, radiation pressure, and bipolar jets from young stars can play a secondary role in stirring up gas motions. In the case of the Orion Nebula, evidence for the influence of stellar wind interactions is restricted to the central pc (GarcíaArredondo, Henney, & Arthur, 2001) and the lowdensity western outskirts (Güdel et al., 2008), and they seem to have little influence on the bulk of the nebular gas. Stellar wind effects are more important in older and more massive regions that contain LBV and WolfRayet stars (e.g., Smith & Brooks, 2007). Similarly, radiation pressure, although unimportant in Orion, becomes much more important in higher luminosity regions (Krumholz & Matzner, 2009). HerbigHaro jets and bowshocks dominate the far wings () of the velocity distribution in Orion (Henney et al., 2007), but the total kinetic energy of these high velocity flows is relatively low, so that the effect on the global velocity statistics is minor.
6 Summary
We have used statistical analysis of highresolution spectroscopic observations of optical emission lines in the central pc () of the Orion Nebula in order to characterize the turbulence in the ionized gas. The analysis has been guided and informed by radiation hydrodynamic simulations of \ionHii region evolution. The techniques that we have applied are:

Secondorder structure function of velocity centroids (§ 2.2.2), which gives the variation as a function of planeofsky separation of the differences in average lineofsight velocity.

Velocity channel analysis (VCA; § 2.2.1), which compares the spatial power spectrum slope of velocityresolved and velocityintegrated emission profiles of the same line.

Linewidth analysis (§ 3.4), which is sensitive to velocity fluctuations along the line of sight

Probability density function (PDF; § 4.3) of the surface brightness in different lines
Our principal empirical findings are as follows:

The VCA technique is the most reliable means of determining the spectrum of velocity fluctuations in the ionized gas (§ 3.3.1), and we find consistent evidence from both low and high ionization lines for a Kolmogorovtype spectrum () for length scales, , between pc () and pc (). Unfortunately, VCA can not be applied if the thermal or instrumental line width is larger than the velocity differences of interest (Appendix C), which rules out its application to the line and to scales smaller than pc.

The structure functions show systematic trends with degree of ionization (§ 3.3.3). Higher ionization lines tend to show larger autocorrelation scales, larger total planeofsky velocity dispersions, and steeper slopes than lower ionization lines. The changes in slopes are difficult to interpret because of the influence of projection smearing and sensitivity to details of the observational methodology (§ 4.1).

The characteristic length of pc is special in at least two ways, corresponding to both the autocorrelation scale of velocity differences for lowionization lines (Figs. 8, 6, 14) and also a break in the power spectrum of surface brightness fluctuations in all lines (Figs. 5–7). We suggest that this is the dominant scale for density fluctuations in the nebula (§ 4.4) and is also the main driving scale of the turbulence. A further break in the surface brightness power spectra occurs at the smaller scale of pc (), but there is no obvious feature in the structure functions at this scale.

Comparison of the application of turbulent diagnostics to numerical simulations (Paper I) with application of the same diagnostics to Orion leads us to conclude (§ 4.2) that even the highionization line emission (eg [OIII]) is confined to a thick shell and does not fill the interior of the nebula. Furthermore, the underlying power spectrum is shallower in the simulations, implying that smallscale turbulent driving is less important in the nebula than it is in the simulations.

There are three lines of evidence suggesting that the velocity fluctuations are not homogeneous on the largest scales, but rather that the turbulent conditions themselves vary, both across the sky and along the line of sight, on scales larger than the velocity autocorrelation length of – pc:

The structure function slope of the [\ionNii] line is significantly steeper in the southern half of our observed field than in the northern half (Fig. 8).

The planeofsky velocity dispersion increases with increasing ionization (Table 5), implying an increasing amplitude of fluctuations towards the interior of the nebula

The lineofsight nonthermal velocity dispersion (after removing the confounding effect of dust scattering; § 3.4) is typically twice as large () as the planeofsky velocity dispersion (). In order to explain this ratio in terms of a homogeneous turbulent layer, the lineofsight depth of the layer would need to be at least 10 times the velocity autocorrelation length, which is unrealistically large (§ 4.3). Instead, the result is more naturally explained by largescale velocity gradients (such as radial expansion), combined with emissivity fluctuations along the line of sight.


The turbulent and ordered components of the velocity dispersion ( and , respectively) are of similar magnitude: estimated to be (§ 4.3).

The PDF of surface brightness fluctuations is approximately lognormal with a fractional width of – in all lines (Fig. 4.3). Turbulent velocity fluctuations can only account for half of the variance in surface brightness. The remaining part may be due to density gradients in photoevaporation flows (§ 4.5).
Acknowledgements
SNXM acknowledges IMPRS for a PhD research scholarship. SJA, WJH and SNXM would like to thank DGAPAUNAM for financial support through projects PAPIIT IN101713 and IN112816. This work has made use of NASA’s Astrophysics Data System.
References
 Arthur et al. (2011) Arthur S. J., Henney W. J., Mellema G., de Colle F., VázquezSemadeni E., 2011, MNRAS, 414, 1747
 Arthur & Hoare (2006) Arthur S. J., Hoare M. G., 2006, ApJS, 165, 283
 Baldwin et al. (1991) Baldwin J. A., Ferland G. J., Martin P. G., Corbin M. R., Cota S. A., Peterson B. M., Slettebak A., 1991, ApJ, 374, 580
 Baldwin et al. (2000) Baldwin, J. A., Verner, E. M., Verner, D. A., Ferland, G. J., Martin, P. G., Korista, K. T., Rubin R. H., 2000, ApJS, 129, 229
 Balick, Gammon, & Hjellming (1974) Balick B., Gammon R. H., Hjellming R. M., 1974, PASP, 86, 616
 BallesterosParedes et al. (2011) BallesterosParedes J., VázquezSemadeni E., Gazol A., Hartmann L. W., Heitsch F., Colín P., 2011, MNRAS, 416, 1436
 Bertoldi & McKee (1990) Bertoldi F., McKee C. F., 1990, ApJ, 354, 529
 Brunt, Federrath, & Price (2010) Brunt C. M., Federrath C., Price D. J., 2010, MNRAS, 403, 1507
 Brunt, Federrath, & Price (2010) Brunt C. M., Federrath C., Price D. J., 2010, MNRAS, 405, L56
 Brunt & Heyer (2002) Brunt C. M., Heyer M. H., 2002, ApJ, 566, 276
 Brunt et al. (2003) Brunt C. M., Heyer M. H., VázquezSemadeni E., Pichardo B., 2003, ApJ, 595, 824
 Brunt & Mac Low (2004) Brunt C. M., Mac Low M.M., 2004, ApJ, 604, 196
 Castañeda (1988) Castañeda H. O., 1988, ApJS, 67, 93
 Castañeda & O’Dell (1987) Castañeda H. O., O’Dell C. R., 1987, ApJ, 315, L55
 Chandrasekhar (1949) Chandrasekhar S., 1949, ApJ, 110, 329
 Deharveng (1973) Deharveng L., 1973, A&A, 29, 341
 Dyson (1968) Dyson J. E., 1968, Ap&SS, 1, 388
 Doi, O’Dell, & Hartigan (2004) Doi T., O’Dell C. R., Hartigan P., 2004, AJ, 127, 3456
 Esquivel & Lazarian (2005) Esquivel A., Lazarian A., 2005, ApJ, 631, 320
 Esquivel et al. (2007) Esquivel A., Lazarian A., Horibe S., Cho J., Ossenkopf V., Stutzki J., 2007, MNRAS, 381, 1733
 Esquivel et al. (2003) Esquivel A., Lazarian A., Pogosyan D., Cho J., 2003, MNRAS, 342, 325
 GarcíaArredondo, Henney, & Arthur (2001) GarcíaArredondo F., Henney W. J., Arthur S. J., 2001, ApJ, 561, 830
 GarcíaDíaz & Henney (2007) GarcíaDíaz M. T., Henney W. J., 2007, AJ, 133, 952
 GarcíaDíaz et al. (2008) GarcíaDíaz M. T., Henney W. J., López J. A., Doi T., 2008, RMxAA, 44, 181
 Güdel et al. (2008) Güdel M., Briggs K. R., Montmerle T., Audard M., Rebull L., Skinner S. L., 2008, Sci, 319, 309
 Henney (1998) Henney W. J., 1998, ApJ, 503, 760
 Henney et al. (2009) Henney W. J., Arthur S. J., de Colle F., Mellema G., 2009, MNRAS, 398, 157
 Henney et al. (2007) Henney W. J., O’Dell C. R., Zapata L. A., GarcíaDíaz M. T., Rodríguez L. F., Robberto M., 2007, AJ, 133, 2192
 Kainulainen et al. (2016) Kainulainen J., Stutz A. M., Stanke T., AbreuVicente J., Beuther H., Henning T., Johnston K. G., Megeath T., 2016, arXiv, arXiv:1603.05688
 Kolmogorov (1941) Kolmogorov A., 1941, DoSSR, 30, 301
 Konstandin et al. (2012) Konstandin L., Girichidis P., Federrath C., Klessen R. S., 2012, ApJ, 761, 149
 Konstandin et al. (2015) Konstandin L., Schmidt W., Girichidis P., Peters T., Shetty R., Klessen R. S., 2015, ArXiv eprints
 Kritsuk et al. (2007) Kritsuk A. G., Norman M. L., Padoan P., Wagner R., 2007, ApJ, 665, 416
 Krumholz & Matzner (2009) Krumholz M. R., Matzner C. D., 2009, ApJ, 703, 1352
 Lagrois & Joncas (2011) Lagrois D., Joncas G., 2011, MNRAS, 413, 721
 Lagrois et al. (2011) Lagrois D., Joncas G., Drissen L., Arsenault R., 2011, MNRAS, 413, 705
 Lazarian (2004) Lazarian A., 2004, ASSL, 319, 601
 Lazarian & Pogosyan (2000) Lazarian A., Pogosyan D., 2000, ApJ, 537, 720
 Mc Leod et al. (2016) Mc Leod A. F., Weilbacher P. M., Ginsburg A., Dale J. E., Ramsay S., Testi L., 2016, MNRAS, 455, 4057
 Medina et al. (2014) Medina S.N. X., Arthur S. J., Henney W. J., Mellema G., Gazol A., 2014, MNRAS, 445, 1797 (Paper I)
 Medina Tanco et al. (1997) Medina Tanco G. A., Sabalisck N., JatencoPereira V., Opher R., 1997, ApJ, 487, 163
 Mellema et al. (2006a) Mellema G., Arthur S. J., Henney W. J., Iliev I. T., Shapiro P. R., 2006, ApJ, 647, 397
 Mellema et al. (2006b) Mellema G., Iliev I. T., Alvarez M. A., Shapiro P. R., 2006, NewA, 11, 374
 MivilleDeschênes, Joncas, & Durand (1995) MivilleDeschênes M.A., Joncas G., Durand D., 1995, ApJ, 454, 316
 MivilleDeschênes, Levrier, & Falgarone (2003) MivilleDeschênes M.A., Levrier F., Falgarone E., 2003, ApJ, 593, 831
 Münch (1958) Münch G., 1958, RvMP, 30, 1035
 O’Dell (1986) O’Dell C. R., 1986, ApJ, 304, 767
 O’Dell (2001) O’Dell C. R., 2001, ARA&A, 39, 99
 O’Dell & Castañeda (1987) O’Dell C. R., Castañeda H. O., 1987, ApJ, 317, 686
 O’Dell et al. (2015) O’Dell C. R., Ferland G. J., Henney W. J., Peimbert M., GarcíaDíaz M. T., Rubin R. H., 2015, AJ, 150, 108
 O’Dell & Harris (2010) O’Dell C. R., Harris J. A., 2010, AJ, 140, 985
 O’Dell & Henney (2008) O’Dell C. R., Henney W. J., 2008, AJ, 136, 1566
 O’Dell et al. (2009) O’Dell C. R., Henney W. J., Abel N. P., Ferland G. J., Arthur S. J., 2009, AJ, 137, 367
 O’Dell & Townsley (1988) O’Dell C. R., Townsley L. K., 1988, A&A, 198, 283
 O’Dell, Townsley, & Castañeda (1987) O’Dell C. R., Townsley L. K., Castañeda H. O., 1987, ApJ, 317, 676
 O’Dell & Wen (1992) O’Dell C. R., Wen Z., 1992, ApJ, 387, 229
 O’Dell, Wen, & Hu (1993) O’Dell C. R., Wen Z., Hu X., 1993, ApJ, 410, 696
 O’Dell & YusefZadeh (2000) O’Dell C. R., YusefZadeh F., 2000, AJ, 120, 382
 Padoan & Nordlund (2002) Padoan P., Nordlund Å., 2002, ApJ, 576, 870
 Pogge, Owen, & Atwood (1992) Pogge R. W., Owen J. M., Atwood B., 1992, ApJ, 399, 147
 Press et al. (1992) Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numerical recipes in FORTRAN. The art of scientific computing, Cambridge: University Press
 Roy & Joncas (1985) Roy J.R., Joncas G., 1985, ApJ, 288, 142
 Scalo (1984) Scalo J. M., 1984, ApJ, 277, 556
 Schlafly et al. (2014) Schlafly E. F., et al., 2014, ApJ, 786, 29
 Smith & Brooks (2007) Smith N., Brooks K. J., 2007, MNRAS, 379, 1279
 Stanimirović & Lazarian (2001) Stanimirović S., Lazarian A., 2001, ApJ, 551, L53
 TenorioTagle (1979) TenorioTagle G., 1979, A&A, 71, 59
 Tobin et al. (2009) Tobin J. J., Hartmann L., Furesz G., Mateo M., Megeath S. T., 2009, ApJ, 697, 11031118
 VázquezSemadeni, Kim, & BallesterosParedes (2005) VázquezSemadeni E., Kim J., BallesterosParedes J., 2005, ApJ, 630, L49
 von Hoerner (1951) von Hoerner S., 1951, ZA, 30, 17
 Weilbacher et al. (2015) Weilbacher P. M., et al., 2015, A&A, 582, A114
 Wen & O’Dell (1993) Wen Z., O’Dell C. R., 1993, ApJ, 409, 262
 Zuckerman (1973) Zuckerman B., 1973, ApJ, 183, 863
Appendix A Positions of highresolution [S Ii] observations
The [\ionSii] 6716, 6731Å observations consist of 92 NorthSouth pointings, 37 observed at KPNO with the same characteristics as the H, [\ionNii] and [\ionOiii] observations used in this paper and 55 observed at OANSPM. Of this latter group, 20 pointings have a high velocity resolution of 6 km s and the remaining 35 pointings have a lower velocity resolution of 12 km s. All 92 datasets are used to construct the velocity moment maps shown in Figure 3. However, only the 20 highest resolution observations are useful for the VCA calculations, since the other observations are too affected by noise at high wavenumber. The positions of the longslits for these 20 observations are shown in Figure 17
Appendix B Positions of [O Iii] horizontal slits
The main dataset consists of NorthSouth oriented longslit spectra. We can test whether orientation affects our velocity channel analysis by examining a supplementary data set for the [\ionOiii] 5007 emission line consisting of observations perpendicular to the main data set (i.e., oriented EastWest) made at OANSPM. Fifteen spectra were obtained in steps of starting at south of Ori C and proceeding south. The slit positions and orientations are indicated in Figure 18 and the observations are described in O’Dell et al. (2015).
Appendix C Are the thin slices really thin?
Our analysis of the powerlaw indices in Section 3.3 relates the difference between the spectral index obtained from the thin slices and that of the thickest velocity slice to the underlying 3D powerspectrum index of the velocity fluctuations. We based our choice of thin slice ( km s) on the thinnest slice not affected by instrumental broadening or thermal broadening. Esquivel et al. (2003) discuss the criteria for a velocity channel to be thin or thick, which depends not only on the number of channels into which the total velocity range is divided, but also on the scale. For channels spanning the velocity range a channel is “just thin enough” if
(7) 
where is a reference scale, is the dispersion of velocities of points separated by this distance, is the scale of interest and is the spectral index of the underlying velocity fluctuations (3D structure function). The problem here is that both and refer to the 3D velocity distribution, while we have access to only the lineofsight velocities and can calculate only the 2D structure function of the velocity centroids. We do, however, also have the lineofsight velocity linewidths and, as discussed in Section 3.4 we can estimate the turbulent linewidth (actually, we use , the mean lineofsight velocity dispersion shown in Figure 13). Assuming that the lineofsight turbulent linewidth is related to the velocity dispersion at the largest scales before there is decorrelation, then we have and the scale corresponds to the local maximum in the secondorder structure function, i.e the turnover point. For the spectral index, , we use the Kolmogorov theoretical value , which we also derive from the VCA.
In particular, for the [\ionOiii] 5007 emission line we obtain km s and arcsec. Putting everything together, we can rearrange Equation 7 to find the minimum scale for which our km s velocity slices are indeed thin:
(8) 
giving arcsec. This is consistent with our VCA (see Fig. 5), where the spectral indices of the thin velocity slices cease to make sense (they become steeper than the thick slice slopes) for scales smaller than arcsec. For the range corresponding to regime II of the VCA, our thin velocity slices with km s are in fact “thin enough” for the [\ionOiii] emission line.
For the [\ionSii] emission lines, the value of km s and the turnover in the structure function occurs at arcsec and we obtain a minimum scale arcsec, which again is consistent with our VCA (see Fig. 6). Finally, for the [\ionNii] 6583 emission line, we find km s and the turnover scale is arcsec, from which we obtain arcsec, which again is entirely consistent with our velocity channel analysis (see Fig. 5).
Appendix D Planeofsky versus lineofsight variations in velocity from simulated H Ii region
For two representative times, Figure 19 shows the joint distribution of mean velocity and nonthermal line width , calculated from the line profiles of synthetic positionpositionvelocity cubes for the simulated turbulent \ionHii region of Medina et al. (2014). Figure 20 shows the temporal evolution of the planeofsky average and standard deviation of the same two quantities for the [\ionOiii] line. Results from the upperright and lowerleft panel of Figure 20 respectively provide the data that go into the horizontal and vertical axes of Figure 13(b) in § 3.4.
At late times (Fig. 19c,d), when dust absorption is relatively unimportant, the average line centroid velocity (upper left panel of Fig. 20) reflects the champagne flow due to the largestscale density gradients in our simulation box. This leads to both blue and red shifts, since opposite viewing directions (e.g., and ) have roughly equal but opposite mean velocities, so that pairs of PDFs in Fig. 19c,d are rough mirror images. At earlier times (Fig. 19a,b), radial expansion dominates and the dust optical depth is larger, which leads to selectively greater absorption of receding regions of the nebula. This produces predominantly blueshifted mean velocities from all viewing directions.
Appendix E Toy model of surface brightness profiles
To try and gain insight into the observed power spectra, we have determined the simplest emission model that can reproduce the qualitative features of regimes I–IV. The emission model is illustrated in Figure 21 and consists of the following spatial components:

A narrow Gaussian with peak brightness 800 () and FWHM ().

A broader Gaussian with peak brightness 400 () and FWHM ().

A very broad Gaussian with peak brightness 240 () and FWHM ().

A sinusoidal variation with wavelength () and relative amplitude 15%, which multiplies the sum of components (i)–(iii).

Poisson noise, assuming that the brightness is in units of counts/pixel.
The wavenumber of each component is indicated by vertical lines on the power spectrum shown in the upper panel of Figure 21 (dashed lines for the Gaussian components (i)–(iii), dotted lines for the sinusoidal component (iv)). For the Gaussian components, we assume since the Gaussian peak represents only half of a fluctuation wavelength.
Components (i) and (ii) are chosen so as to reproduce regime II in the power spectra of § 3.1, with a relatively shallow slope () between scales of 8 and 22. As is, the model best represents the [\ionNii] power spectrum (upper right panel of Fig. 5), but the precise value of the slope can be controlled by varying the relative strength of the narrow () and the broader () component. For example, the [\ionOiii] power spectrum requires nearly equal amplitudes for these two components.
With only components (i) and (ii), the model power spectrum is too shallow in regime I (with ) and is too steep in regime III, showing a deep minimum in the curve before arriving at the noisedominated regime IV. These minor deficiencies are ameliorated by the introduction of the secondary components (iii) and (iv), respectively.
The average power spectrum from 100 realizations of the toy model is calculated, where the brightness and width of each component is chosen according to normal distributions with central value ( RMS % variation) as given above. The result is in remarkably good agreement with the observed power spectra, except that the observations tend to show sharper breaks at the boundaries between the different regimes. It should be emphasized that the toy model is only an idealization of the spatial variations present in the real surface brightness maps, which typically show prominent peaks in each onedimensional slit profile, rather than the single peak of the model.