Turbulence in the Ionized Gas of the Orion Nebula

# Turbulence in the Ionized Gas of the Orion Nebula

S. J. Arthur, S.-N. X. Medina, W. J. Henney
Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, Campus Morelia, 58090 Morelia, Michoacán, México.
Max Planck Institüt für Radioastronomie, Bonn, Germany.
E-mail: j.arthur@crya.unam.mx
###### 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 non-thermal 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 root-mean-square 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 — turbulence
pagerange: Turbulence in the Ionized Gas of the Orion NebulaEpubyear: 2016

## 1 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 point-to-point radial velocities with scale (for example, Münch, 1958; Roy & Joncas, 1985; O’Dell & Castañeda, 1987; Miville-Deschê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 three-dimensional 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 plane-of-sky 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 high-velocity-resolution 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 Fabry-Perot interferometry (Roy & Joncas, 1985; Miville-Deschê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 B-type 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 (OMC-1) and there is a large population of young stars, some of which are the sources of stellar jets and Herbig-Haro 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; Tenorio-Tagle, 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 power-law 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 emission-line 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 top-down turbulence driven at the largest scales.

In this paper, we apply both velocity channel analysis of position-velocity arrays and structure functions of velocity centroid maps to high spatial and velocity resolution observations of low- and high-ionization 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ía-Díaz & Henney (2007) and García-Dí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 position-velocity 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ía-Dí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 (OAN-SPM), Mexico. The data covers the central (Huygens) region of the Orion Nebula111 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 North-South 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 North-South pointings covering the same region, consisting of 37 positions observed at KPNO and 55 positions observed at OAN-SPM. The KPNO [\ionSii] data consists of two disjoint regions, one in the east and one in the west of the nebula, while the OAN-SPM 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 OAN-SPM is .

In the present paper, we use the calibrated position-velocity (PV) arrays obtained by García-Díaz et al. (2008) from each individual longslit spectrum for the velocity channel analysis of each emission line. For the second-order 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ía-Dí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 second-order 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 emission-line 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 Second-order Structure Functions

The second-order structure function of the velocity centroids is (see Paper I)

 S2(l)=∑pairs[Vc(r)−Vc(r+l)]2σ2vcN(l) . (1)

In this definition, is the two-dimensional 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

 σ2vc≡∑pixels[Vc(r)−⟨Vc⟩]2N . (2)

Here, is the mean centroid velocity

 ⟨Vc⟩≡∑pixelsVc(r)N . (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

 S2(l)=∑pairs[Vc(r)−Vc(r+l)]2I(r)I(r+l)σ2vcW(l) , (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 small-scale, high velocity features such as jets and Herbig-Haro 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 H6583, [\ionNii] 6853, [\ionOiii] 5007, [\ionSii] 6716 and [\ionSii] 6731 lines from the emission line atlas published by García-Dí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 ill-defined boundaries. We expect variation in the power-law indices obtained from both the velocity channel analysis and the second-order structure functions depending on position in the nebula due to large-scale inhomogeneities. We estimate the confidence bounds on the power-law 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 power-law 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 high-resolution [\ionSii] slits) are obtained using a least-squares fitting procedure. The variation of these sample power-law indices provides an estimate of the confidence bounds for the power-law index across the Orion Nebula.

Furthermore, our longslit spectra are all oriented North-South and so in order to check whether this has any effect on the power-law index, we analyze a supplementary PV dataset of [\ionOiii] 5007 observations from 18 slit positions perpendicular to the main data set (i.e., oriented East-West). The slit positions and orientations are indicated in Figure 18 (O’Dell et al., 2015).

We estimate the effects of large-scale spatial inhomogeneity on the second-order structure function by evaluating (see Eqs. 1 and 4) for 100 distinct, randomly selected rectangular frames within each 2D map. Each frame has two-thirds the and dimension of the initial map (see Fig. 4). The variation in power-law indices of the smaller maps give us an estimate of the spread of power-law 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 second-order 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ía-Dí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 power-law 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 power-law 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 power-law 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 power-law 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 large-scale inhomogeneity of the power spectrum depending on the position of the observations.

### 3.2 Second-order Structure Functions

#### 3.2.1 Shape and Slope of Structure Function for each emission line

We use Equation 1 to calculate the second-order 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 pdf-selected pixels, together with the structure functions obtained from 100 smaller rectangular frames within the velocity centroid map (see Fig. 4). The fit to the power-law 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 one-pixel 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 power-law 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 inter-slit separation become important.

The steepest power-law 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 power-law 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 power-law 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 emission-line velocity centroid maps. The solid red circles show the same pdf-selected 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 high-velocity 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 number-of-points-weighted structure functions (solid symbols). However, this time the [\ionNii] structure functions are different, with the intensity weighting producing structure functions with shallower slopes (smaller power-law indices). For the shorter wavelength [\ionSii] 6716 case, the intensity weighting makes little difference to the pdf-selected structure function (red circles) but does reduce the noise contribution to the full structure function (blue triangles) and leads to a similar power-law index to the pdf-selected cases over the fitted range. For the longer wavelength [\ionSii] 6731 case, the intensity weighting produces a shallower power law, even for the pdf-selected 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 power-law indices

We want to use the power-law indices from the velocity channel analysis and the second-order structure functions obtained from the spectroscopic observations to recover the three-dimensional velocity statistics of the ionized gas in the Orion Nebula. There are various relationships between the observationally derived power-law 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 power-law indices of the average power spectra of thin and thick velocity slices. The thickest velocity slice corresponds to the velocity-integrated surface brightness and the power-law 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 power-law index.

There are relations between the spectral index of the 3D velocity fluctuation power spectrum and the power-law 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 power-law index of the average power spectrum of the thin isovelocity channels is given by

 γt=−3+12m3D , (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 power-law indices of the average power-spectra of the thin and thick velocity slices:

 γt−γT=12m3D . (6)

There are also relationships between the power-law indices of the two-dimensional structure function of the emission-line velocity centroid map and the three-dimensional 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 sheet-like, i.e., essentially two-dimensional, then (Castañeda & O’Dell, 1987; Miville-Deschênes, Joncas, & Durand, 1995). The three-dimensional structure function power-law index is related to the three-dimensional 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 power-law 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 power-law indices listed in Table 1 (see also Eq. 5), we obtain and for the short and long wavelength [\ionSii] lines, respectively. The corresponding power-law 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 power-law 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 power-law 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 power-law 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 power-law indices steeper than their ‘thin’ slice power-law indices. The [\ionNii], H and [\ionOiii] have curiously opposite behaviour and their power-law 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 power-law 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 structure-function index must lie between the two limits given by projection smoothing and a sheet-like 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 Line-of-sight 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 non-thermal line-of-sight velocity dispersion may be estimated by correcting the observed line widths for the contributions from the spectrograph resolution, fine-structure splitting, and thermal Doppler broadening, as in Eq. (2) of García-Díaz et al. (2008). A striking fact about the nebular line profiles is that the line-of-sight velocity dispersion is several times larger than the plane-of-sky dispersion in mean velocities. This is illustrated in Figure 10, which shows flux-weighted histograms of the non-thermal RMS line widths versus mean velocity for all lines of the García-Díaz et al. (2008) spectral atlas. The line-of-sight RMS velocity dispersion is 9–10 , whereas the RMS plane-of-sky dispersion of the mean velocities is only 2–4 .

However, the observed line widths are affected by additional broadening due to dust-scattering (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 line-of-sight 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 back-scattered 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 FWHM-based 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 line-of-sight line widths with respect to the moment-derived values of Figure 10. The exception is , where larger thermal broadening means that is still contaminated by the back-scattered 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 low-ionization 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.

## 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 scale-dependent term to the structure function.

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 flux-weighted 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 multi-component 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 three-Gaussian decomposition of our line profiles (component A, see § 3.4), with results that are also included in Table 5. We find a slightly larger plane-of-sky 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 under-determination 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 fixed-pattern noise in their mean velocity maps, which can be seen as tartan-like 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 second-order 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 second-order 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 power-law 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 power-law indices of the second-order 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 log-log space) whereas the observational structure function has a slowly varying power law as a function of separation scale.

Both the structure function power-law 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] thick-slice 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 power-law 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 structure-function 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 low-mass stars, some of which are sources of jets and Herbig-Haro 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 more-or-less 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 non-thermal 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 plane-of-sky 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 FWHM-derived and reflection-corrected 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 line-of-sight velocity dispersion is roughly twice the plane-of-sky 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 line-of-sight depth reduces the plane-of-sky 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 line-of-sight 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ía-Dí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 high-ionization lines and for the low-ionization lines. Therefore, projection smearing of the large-scale fluctuations is negligible and cannot explain the difference between and .

A second, contrasting interpretation of the evidence would be in terms of large-scale, 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 red-shifted and blue-shifted hemispheres. The required RMS fractional variation in the emissivity on the scale of the shell diameter is found to be .

The observed large-scale () brightness fluctuations are illustrated in Figure 15, which shows log-normal fits to the PDFs of surface brightness, , after normalizing by the mean, , and binning the maps at pixels. The RMS width of the log-normal 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 line-of-sight 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 plane-of-sky 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 large-scale 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 blue-shifts 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ía-Díaz et al., 2008), which gives a champagne-flow 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 small-scale fluctuations are relatively unimportant. Above 22 arcsec (regime I), the power spectrum is very flat, similar to the noise-dominated 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.222The 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 “de-projection 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 de-projection 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 top-down 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 thin-slice 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 bright-rimmed 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 emissivity-weighted density PDF for a simple model of a single spherically divergent, isothermal evaporation flow from a D-critical 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; Ballesteros-Paredes 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 one-sided 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 & Yusef-Zadeh, 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ía-Arredondo, Henney, & Arthur, 2001) and the low-density 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 Wolf-Rayet stars (e.g., Smith & Brooks, 2007). Similarly, radiation pressure, although unimportant in Orion, becomes much more important in higher luminosity regions (Krumholz & Matzner, 2009). Herbig-Haro 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 high-resolution 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:

1. Second-order structure function of velocity centroids (§ 2.2.2), which gives the variation as a function of plane-of-sky separation of the differences in average line-of-sight velocity.

2. Velocity channel analysis (VCA; § 2.2.1), which compares the spatial power spectrum slope of velocity-resolved and velocity-integrated emission profiles of the same line.

3. Line-width analysis (§ 3.4), which is sensitive to velocity fluctuations along the line of sight

4. Probability density function (PDF; § 4.3) of the surface brightness in different lines

Our principal empirical findings are as follows:

1. 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 Kolmogorov-type 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.

2. The structure functions show systematic trends with degree of ionization (§ 3.3.3). Higher ionization lines tend to show larger autocorrelation scales, larger total plane-of-sky 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).

3. The characteristic length of  pc is special in at least two ways, corresponding to both the autocorrelation scale of velocity differences for low-ionization lines (Figs. 8, 6, 14) and also a break in the power spectrum of surface brightness fluctuations in all lines (Figs. 57). 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.

4. 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 high-ionization 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 small-scale turbulent driving is less important in the nebula than it is in the simulations.

5. 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:

1. 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).

2. The plane-of-sky velocity dispersion increases with increasing ionization (Table 5), implying an increasing amplitude of fluctuations towards the interior of the nebula

3. The line-of-sight non-thermal velocity dispersion (after removing the confounding effect of dust scattering; § 3.4) is typically twice as large () as the plane-of-sky velocity dispersion (). In order to explain this ratio in terms of a homogeneous turbulent layer, the line-of-sight 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 large-scale velocity gradients (such as radial expansion), combined with emissivity fluctuations along the line of sight.

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

7. The PDF of surface brightness fluctuations is approximately log-normal 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 DGAPA-UNAM 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ázquez-Semadeni 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
• Ballesteros-Paredes et al. (2011) Ballesteros-Paredes J., Vázquez-Semadeni 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ázquez-Semadeni 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ía-Arredondo, Henney, & Arthur (2001) García-Arredondo F., Henney W. J., Arthur S. J., 2001, ApJ, 561, 830
• García-Díaz & Henney (2007) García-Díaz M. T., Henney W. J., 2007, AJ, 133, 952
• García-Díaz et al. (2008) García-Dí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ía-Díaz M. T., Rodríguez L. F., Robberto M., 2007, AJ, 133, 2192
• Kainulainen et al. (2016) Kainulainen J., Stutz A. M., Stanke T., Abreu-Vicente 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 e-prints
• 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., Jatenco-Pereira 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
• Miville-Deschênes, Joncas, & Durand (1995) Miville-Deschênes M.-A., Joncas G., Durand D., 1995, ApJ, 454, 316
• Miville-Deschênes, Levrier, & Falgarone (2003) Miville-Deschê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ía-Dí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 & Yusef-Zadeh (2000) O’Dell C. R., Yusef-Zadeh 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
• Tenorio-Tagle (1979) Tenorio-Tagle 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, 1103-1118
• Vázquez-Semadeni, Kim, & Ballesteros-Paredes (2005) Vázquez-Semadeni E., Kim J., Ballesteros-Paredes 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 high-resolution [S Ii] observations

The [\ionSii] 6716, 6731Å observations consist of 92 North-South pointings, 37 observed at KPNO with the same characteristics as the H, [\ionNii] and [\ionOiii] observations used in this paper and 55 observed at OAN-SPM. 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 North-South 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 East-West) made at OAN-SPM. 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 power-law 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 power-spectrum 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

 N>ΔvσL(Lr)m/2 (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 line-of-sight velocities and can calculate only the 2D structure function of the velocity centroids. We do, however, also have the line-of-sight velocity linewidths and, as discussed in Section 3.4 we can estimate the turbulent linewidth (actually, we use , the mean line-of-sight velocity dispersion shown in Figure 13). Assuming that the line-of-sight 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 second-order 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:

 r>(δvσL)2/mL (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 Plane-of-sky versus line-of-sight variations in velocity from simulated H Ii region

For two representative times, Figure 19 shows the joint distribution of mean velocity and non-thermal line width , calculated from the line profiles of synthetic position-position-velocity cubes for the simulated turbulent \ionHii region of Medina et al. (2014). Figure 20 shows the temporal evolution of the plane-of-sky average and standard deviation of the same two quantities for the [\ionOiii] line. Results from the upper-right and lower-left 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 largest-scale 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 blue-shifted 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:

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

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

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

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

5. 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 noise-dominated 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 one-dimensional slit profile, rather than the single peak of the model.

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters