---
author:
- 'Ma.  T. García-Díaz, W. J. Henney, J. A. López, and T. Doi'
bibliography:
- 'All-Sorted.bib'
- 'Raga.bib'
- 'Extras.bib'
- 'ads-orion.bib'
title: |
    Velocity Structure in the Orion Nebula.\
    II.  Emission Line Atlas of Partially Ionized\
    to Fully Ionized Gas
---

Introduction {#sec:introduction}
============

The Orion nebula, M42, is the nearest and best studied high-mass star-forming region in our galaxy .  It is being ionized by a compact group of high-mass stars known as the Trapezium, which was created in the first of three recent epsiodes of star formation that have occured within a region of size one-tenth of a parsec at the center of the Orion Nebula Cluster.  The other two star formation episodes, associated with the BN/KL and Orion South regions, are still ongoing and are ocurring in dense molecular gas that lies behind the ionized nebula.  The O7 V star [$\theta^1$ Ori C]{} is the most luminous and hottest star of the Trapezium and it is responsible for the main ionization front.  Many morphological features in Orion that are apparent in emission lines, such as arcs and filaments are associated with collimated gas flows from star-forming activity, such as Herbig-Haro (HH) objects .  Jets, bright bars, proplyds and a myriad of other features are conspicuous in optical emission lines making of Orion an ideal laboratory for studying photoionization and hydrodynamic phenomena.  Previous high spectral resolution studies over restricted regions in Orion [e.g., @1988ApJS...67...93C; @1993ApJ...409..262W; @1998MNRAS.295..401E; @2000ApJS..129..229B] indicate the presence of complex line profiles; the changing complexity of the emission line profiles over distinct regions of Orion can now be fully appreciated and clearly identified with various phenomena over the entire region covered by our observations. 

In this paper we present a comprehensive atlas of high-resolution echelle slit spectra in multiple emission lines, covering a $3'\times
5'$ region in the center of the Orion nebula.  The atlas combines results from two separate datasets, one obtained at Kitt Peak National Observatory (KPNO), and the other from the Observatorio Astronómico Nacional at San Pedro Mártir, B.C., Mexico (SPM).  Between them, these datasets, which have comparable spatial and spectral resolution, span a wide range of ionization stages, going from [[[\[\]]{}6300Å]{}]{}, which traces partially ionized gas at the ionization front, up to [[[\[\]]{}5007Å]{}]{}, which traces the most highly ionized gas in the interior of the nebula.  Figures and data from the atlas are publicly available from the web page <http://www.astrosmo.unam.mx/~w.henney/orionatlas>. 

Various aspects of the data presented here have been discussed in previous papers: @2004AJ....127.3456D [hereafter Paper I] concentrated on the high-velocity features associated with stellar jets; @2005ApJ...627..813H analysed the large-scale nebular structure and kinematics in terms of numerical simulations of champagne flows; @2007AJ....133..952G [hereafter Paper I] investigated the global structures found in low-velocity emission from low-ionization gas; @2007AJ....133.2192H analysed the structure and dynamics of the low-excitation jet, [HH 528]{}.  However the full set of individual long-slit spectra are shown here for the first time and the usefulness and applicability of this comprehensive data set goes beyond the analysis of the previous papers.  As example applications of the atlas data, we calculate velocity moment maps of the nebula in each emission line and present a comprehensive analysis of the line profile statistics. 

The structure of the paper is as follows.  In § \[sec:atlas\] we present the atlas of slit spectra.  In § \[sec:moment-maps\] we summarise the spatial variations in line profiles across the face of the nebula by means of maps of integrated line quantities, such as centroid velocity and linewidth.  In § \[sec:velocity-statistics\] we analyse the statistical properties of the velocity field in the nebula, including correlations between the kinematics in different emission lines.  In § \[sec:deriv-mean-electr\] we derive the nebular electron temperature from the difference in linewidths of [$\mathrm{H}\alpha$]{} and [\[\]]{}.  In § \[sec:discussion\] we discuss our results in the light of previous work.  in § \[sec:conclusions\] we present our conclusions. 

The Atlas {#sec:atlas}
=========

![ image](siil_01x)

![ image](siil_02x)

![ image](oi_01x)

![ image](oi_02x)

![ image](siii_01x)

![ image](ha_01x)

![ image](ha_02x)

![ image](nii_01x)

![ image](nii_02x)

![ image](oiii_01x)

![ image](oiii_02x)

  ------------------------------------- ------ ----- -------
                                                     
  Line                                   KPNO   SPM   Total
  [[[\[\]]{}6300Å]{}]{}                         60     60
  [[[\[\]]{}6716,6731Å]{}]{}              37    55     92
  [[[\[\]]{}6584Å]{}]{}                   96           96
  [[[\[\]]{}6312Å]{}]{}                         60     60
  [[[$\mathrm{H}\alpha$]{}6563Å]{}]{}     96           96
  [[[\[\]]{}5007Å]{}]{}                   96           96
  ------------------------------------- ------ ----- -------

  : Emission lines in the atlas[]{data-label="tab:atlas"}

The spectral line data used to construct the atlas come from observations at two different observatories, KPNO and SPM.  The number of north-south oriented slits from each observatory that contribute to the dataset for each emission line is given in Table \[tab:atlas\]. 

The SPM observations were obtained with the Manchester echelle spectrometer [MES; @2003RMxAA..39..185M] combined with the f/7.8 focus on the 2.1 m San Pedro Mártir telescope of the Universidad Nacional Autónoma de México, further details of the observations are described in ).  The KPNO observations were obtained with the echelle spectrometer attached to the f/8 Cassegrain focus of the Mayall 4 m telescope at Kitt Peak National Observatory, further details are given in @2001ApJ...556..203O and @2004AJ....127.3456D.

\[recalibration\] The SPM-MES spectra were obtained with a slit width of 150 $\mu$m (2.0) for a resolution of 12 [$\mathrm{km\ s}^{-1}$]{}, although a few positions were observed with a narrower, 70 $\mu$m (0.9) slit and a 6 [$\mathrm{km\ s}^{-1}$]{}, resolution.  The slit length at SPM is 312.  For KPNO the slit width is 130 $\mu$m (0.8), the slit length 300 and 8 [$\mathrm{km\ s}^{-1}$]{} resolution.  Further details of the data reduction steps are given in and .  Since the original analysis of was oriented towards fast-moving jets and shocks, the tolerances employed for their velocity calibration ($\simeq 4~{\ensuremath{\mathrm{km\ s}^{-1}}}{}$) did not need to be as tight as in .  Furthermore, although employed an E-W slit in order to tie together the velocity calibrations of the multiple N-S slits, this was not done by .  Instead, it was assumed that the mean velocity, integrated over each slit, was constant as a function of right ascension. 

We have therefore internally re-calibrated each velocity cube in two ways.  First, we have applied small rectification corrections in order to remove obvious discontinuities between adjacent slits, which become apparent when examining the kinematics of slow-moving gas.  For [$\mathrm{H}\alpha$]{} and [\[\]]{}, the affected slits are all those lying between $x =
-55$ and $x = 47$, whereas for [\[\]]{} no corrections were necessary.  A common velocity correction, which is linear in distance along the slit, is applied to all the affected slits in each emission line: $\delta V = a + b y$ ($\delta V$ is in [$\mathrm{km\ s}^{-1}$]{}, while $x$ and $y$ are in arcsec with respect to [$\theta^1$ Ori C]{}), with $a = -0.7$, $b = -0.013$ for [$\mathrm{H}\alpha$]{} and $a = -0.4$, $b = -0.015$ for [\[\]]{}.  The velocity corrections to individual pixels vary between $\delta V \simeq
-3~{\ensuremath{\mathrm{km\ s}^{-1}}}{}$ in the far south to $\delta V \simeq +1~{\ensuremath{\mathrm{km\ s}^{-1}}}{}$ in the far north.  Although this procedure greatly reduced the discontinuities, it did not remove them completely in the case of [$\mathrm{H}\alpha$]{}. 

For the second recalibration, in order to correct large-scale trends in the mean velocity as a function of right ascension, we have obtained new observations in February 2007 through an E-W oriented slit, positioned 22S of [$\theta^1$ Ori C]{}.  These observations were obtained with SPM-MES, using a 70 $\mu$m slit and an instrumental configuration identical to that of , except for the filters employed.  The mean velocities of [\[\]]{}, [$\mathrm{H}\alpha$]{}, and [\[\]]{} were calculated as a function of RA along the slit, and these were compared with the mean velocities extracted from the velocity cubes for the same slit.  For [\[\]]{}, we found an rms difference $< 1~{\ensuremath{\mathrm{km\ s}^{-1}}}{}$, so no correction was applied.  However, for [$\mathrm{H}\alpha$]{} and [\[\]]{}, systematic discrepancies of $\sim 5~{\ensuremath{\mathrm{km\ s}^{-1}}}{}$ are seen in the west of the nebula.  We have therefore fitted a third order polynomial to the velocity differences as a function of RA and used it to correct each slit of the cube. 

In order to check the calibration of the absolute velocity scale for each emission line, we have compared our results with those from previous studies [@1980PASP...92...22B; @1988ApJS...67...93C; @1993ApJ...409..262W; @1999AJ....118.2350H; @2000ApJS..129..229B].  In all cases, we found agreement to within $2~{\ensuremath{\mathrm{km\ s}^{-1}}}$.

In summary, we estimate that over most of the [\[\]]{}, [\[\]]{}, [$\mathrm{H}\alpha$]{}, and [\[\]]{} maps, the relative velocities of each emission line are accurate to $\sim 0.5~{\ensuremath{\mathrm{km\ s}^{-1}}}{}$ over relatively short distances ($<
50''$) in the declination direction, but only accurate to $\simeq
1~{\ensuremath{\mathrm{km\ s}^{-1}}}{}$ when comparing points lying in different slits or widely separated in declination.  For the [\[\]]{} and [\[\]]{} lines, the low signal-to-noise ratio over much of the map means that the accuracy is much worse than this for individual spatial pixels, and therefore the same velocity precision can only be retrieved by integrating the line profile over an extended region and thus losing spatial resolution.  For the absolute value of the mean velocity of the different lines, we also estimate an uncertainty of $1$–$2~{\ensuremath{\mathrm{km\ s}^{-1}}}{}$.

Figures \[fig:mosaico\_a\]–\[fig:mosaico\_e\] show the individual long-slit spectra in the form of calibrated position-velocity (P-V) arrays.  The HST image of Orion at the left of the figures shows with a highlighted strip the region covered by the sixteen slit positions shown to the right.  The strip advances to the right (west) in the following rows indicating the corresponding regions for the spectra.  Figures 1 & 2 correspond to the [[[\[\]]{}6731Å]{}]{} emission line.  We do not show here the data for the [[[\[\]]{}6716Å]{}]{} component, although the [\[\]]{} line ratio has been used by to derive the spatially resolved electron density for the whole region.  Strips A–F in Figures \[fig:mosaico\_a\] and \[fig:mosaico\_b\] cover the full survey area for [\[\]]{}.  Figures \[fig:mosaico\_c\] and \[fig:mosaico\_d\] (Strips A–C) correspond to [[[\[\]]{}6300Å]{}]{} and finally Figures \[fig:mosaico\_d\] (Strip D) and \[fig:mosaico\_e\] (Strips A–B) correspond to [[[\[\]]{}6731Å]{}]{}. 

The spectra are presented within a heliocentric velocity range of $-40$ to $+60~{\ensuremath{\mathrm{km\ s}^{-1}}}{}$ in all cases.  We choose this range because it allows one to appreciate the complex structure of the nebula at these relatively low velocities, although some high velocity features, as those from HH jets fall outside the range displayed, but still can be clearly identified.  The declination scale and right ascension scale set the zero referred with respect to [$\theta^1$ Ori C]{}.  The stellar continua have been subtracted from the spectra.  Note that the mean velocities of the background molecular cloud and of the stars in the Orion Nebula Cluster are both ${\ensuremath{V_\odot}}\simeq 27 \pm 2~{\ensuremath{\mathrm{km\ s}^{-1}}}$ and that the conversion to “local standard of rest” (LSR) velocities is $V_\mathrm{LSR} =
{\ensuremath{V_\odot}}- 18.1~{\ensuremath{\mathrm{km\ s}^{-1}}}{}$.

Maps of line profile parameters {#sec:moment-maps}
===============================

![ image](moments-bymoment-sum)

![ image](moments-bymoment-mean)

![ image](moments-bymoment-peak)

![ image](moments-bymoment-sigma)

![ image](moments-bymoment-fwhm)

In order to provide a global overview of the nebular structure and kinematics, we show in Figures \[fig:moments-sum\] to \[fig:moments-fwhm\] maps of derived parameters of the line profile of each emission line at each point in the nebula.  These offer a complementary approach to the isovelocity channel maps that have already been published [@2004AJ....127.3456D; @2007AJ....133..952G; @2007AJ....133.2192H].  The parameters shown are total line surface brightness, $S$, mean velocity, [$\langle V\rangle$]{}, root-mean-square (RMS) linewidth, $\sigma$, peak velocity, [$V_\mathrm{p}$]{}, and full width at half maximum (FWHM) linewidth, $W$. The first three of these quantities can be defined in terms of the first three velocity moments, $M_k$, of the line profiles $I(V)$: $$\label{eq:moments}
  M_k = \int_{V_1}^{V_2} V^k I(V) \, dV ,$$ such that $S = M_0$, ${\ensuremath{\langle V\rangle}}= M_1/M_0$, and $\sigma^2 = (M_2/M_0) -
{\ensuremath{\langle V\rangle}}^2$. The limits of integration are chosen to be $V_1 = -40~{\ensuremath{\mathrm{km\ s}^{-1}}}$ and $V_2 = +70~{\ensuremath{\mathrm{km\ s}^{-1}}}$ for all lines, which represents a compromise between being a wide enough range to include most emission of interest, but not so wide as to introduce too much noise into the maps.  This is is a particular issue for $\sigma$, which is very sensitive to noise in the far wings of the $I(V)$ profile.  The remaining quantities, [$V_\mathrm{p}$]{}, and $W$ are directly measured from the line profile.  In order to determine [$V_\mathrm{p}$]{} to a precision greater than the size of the velocity pixels, we use parabolic interpolation between the peak pixel and its two neighbours.  Similarly, $W$ is calculated as the difference between the velocities of two half-power points, which are found from linear interpolation between adjacent pixels.  In the case that the profile has more than two half-power points, $W$ is defined as the difference between the highest and lowest of these. 

Figure \[fig:moments-sum\] shows the reconstructed surface brightness maps in the six lines.  Comparison with published direct imaging [e.g., @1992ApJ...399..147P] shows generally good agreement, although some artefacts remain from the process of combining the individual spectra.  These can easily be identified as brightness discontinuities across vertical lines in the maps.  The most serious of these artefacts are in the [\[\]]{} map at $x= +10$, $x=+45$, in the [\[\]]{} and [$\mathrm{H}\alpha$]{} maps at $x=+45$, in the [\[\]]{} map at $x=0$, and in the [\[\]]{} map at $x=-5$.

The mean and peak velocities are shown in Figures \[fig:moments-mean\] and \[fig:moments-peak\], where it can be seen that the two quantities are generally well-correlated for a given line, although there are some important differences.  The mean velocity varies smoothly across the face of the nebula, whereas the peak velocity shows frequent discontinuities, which are caused by one particular line component becoming brighter than another.  The mean velocity is also more sensitive to high velocity components, such as HH objects.  This is particularly notable in the case of [\[\]]{}, where the HH objects do not show up at all in the [$V_\mathrm{p}$]{} image.  Since [\[\]]{} and [\[\]]{} are very much weaker than the other lines, the velocity maps are rather noisy, especially in the faint regions of the nebula.  We have therefore masked out from the velocity maps all regions with a surface brightness below a certain threshold, which appear as solid gray areas in the figures. 

The non-thermal contributions to the RMS and FWHM line widths are shown in Figure \[fig:moments-sigma\] and \[fig:moments-fwhm\].  These are calculated by subtracting in quadrature the instrumental and thermal widths from the raw measured widths, assuming a gas temperature of $9100$ K (see § \[sec:deriv-mean-electr\]).  For a Gaussian profile, the relationship between the FWHM and RMS widths is $W = (8 \ln 2)^{1/2} \sigma \simeq 2.355 \sigma$. We therefore multiply the values of $\sigma$ by 2.355 to allow direct comparison between the $\sigma$ and $W$ maps.  There is generally a much greater difference between the $\sigma$ and $W$ maps of a given line than between the maps of [$\langle V\rangle$]{} and [$V_\mathrm{p}$]{}.  The $\sigma$ maps are very much dominated by high-velocity flows, which are much less apparent in the $W$ maps.  As with the [$V_\mathrm{p}$]{} maps, the $W$ maps show discontinuities between patches of the nebula with high and low line widths.  An exception is the [$\mathrm{H}\alpha$]{} line, where the increased thermal broadening leads to a more Gaussian-like line profile and a closer resemblance between the $W$ and $\sigma$ maps.  As with the velocity maps, we have masked out regions of low signal-to-noise ratio for the [\[\]]{} and [\[\]]{} lines. 

An empirical description of the nebular structure and kinematics was given in § 3 of and will not be repeated here.  The reader is referred in particular to Fig.  7 of for orientation in identifying nebular features mentioned in the remainder of the paper. 

One-point velocity statistics {#sec:velocity-statistics}
=============================

                                                line   $\langle V \rangle$   $V_\mathrm{p}$        $\sigma^2$               $W$
  -------------------------------------------------- --------------------- ---------------- ----------------- -----------------
                                    ${[\ion{O}{3}]}$        $16.3 \pm 2.8$   $15.5 \pm 4.0$      $113 \pm 29$   $14.7 \pm 5.1 $
                   ${\ensuremath{\mathrm{H}\alpha}}$        $16.8 \pm 3.0$   $16.2 \pm 4.1$      $195 \pm 28$   $18.6 \pm 5.1 $
                                    ${[\ion{S}{3}]}$        $19.7 \pm 3.0$   $19.2 \pm 4.0$      $99 \pm 46 $   $14.6 \pm 4.4 $
                                    ${[\ion{N}{2}]}$        $20.5 \pm 2.9$   $21.1 \pm 3.7$      $106 \pm 24$   $15.4 \pm 3.7 $
                                    ${[\ion{S}{2}]}$        $21.2 \pm 2.4$   $21.8 \pm 3.2$      $114 \pm 21$   $15.2 \pm 3.5 $
                                    ${[\ion{O}{1}]}$        $25.7 \pm 3.4$   $26.3 \pm 2.5$      $106 \pm 65$   $11.5 \pm 2.6 $
    ${\ensuremath{\mathrm{H}\alpha}}-{[\ion{O}{3}]}$        $0.5 \pm 1.4 $   $0.7 \pm 3.2 $   $82.2 \pm 11.8$     $3.9 \pm 3.0$
                                          $r       $               $0.889$          $0.700$          $0.917 $           $0.825$
                     ${[\ion{N}{2}]}-{[\ion{O}{3}]}$        $4.2 \pm 2.8 $   $5.6 \pm 4.5 $   $-7.0 \pm 36.8$     $0.7 \pm 6.5$
                                          $r       $          $0.511     $     $0.274     $     $0.287      $     $-0.067     $
                     ${[\ion{S}{2}]}-{[\ion{S}{3}]}$        $1.6 \pm 2.5 $   $2.6 \pm 3.3 $   $18.0 \pm 46.8$     $0.2 \pm 5.6$
                                          $r       $          $0.617     $     $0.588     $     $0.229      $          $0.065 $
                     ${[\ion{O}{1}]}-{[\ion{S}{2}]}$        $4.4 \pm 3.0 $   $4.5 \pm 2.7 $   $-7.6 \pm 64.8$   $-3.7 \pm 4.2 $
                                          $r       $          $0.507     $     $0.580     $     $0.175      $     $0.078      $

\
\

\
\

\
\

\
\

In Figures \[fig:meanvsmean\] and \[fig:peakvspeak\] we show joint distributions of the mean and peak velocities, respectively, for selected pairs of emission lines, which illustrate the correlations that exist between the kinematics of the different emission zones of the nebula.  Figure \[fig:sigmavssigma\] shows the joint distributions of the raw velocity variance, $\sigma^2$, which is the square of the RMS velocity dispersion, and Figure \[fig:fwhmvsfwhm\] shows distributions of the non-thermal FWHM linewidth, $W$. The values of $W$ have been approximately corrected for thermal, instrumental and fine-structure broadening in the same way as in Figure \[fig:moments-fwhm\], whereas the values of $\sigma^2$ in Figure \[fig:sigmavssigma\] are uncorrected and measured directly from the data. 

The statistics for individual lines are summarised in the upper part of Table \[tab:stats\], which shows the mean and standard deviation ([s.d.]{})  of each quantity.  These are calculated with an equal weighting of each pixel in the map after masking out any points lying outside the ranges shown in Figures \[fig:meanvsmean\]–\[fig:fwhmvsfwhm\].  For [\[\]]{} and [\[\]]{}, any pixel with surface brightness below the thresholds shown in Figures \[fig:moments-mean\]–\[fig:moments-fwhm\] were also masked out.  The lower part of Table \[tab:stats\] show the mean and [s.d.]{}  of the differences between the kinematic quantities for the line pairs shown in the figures, together with the linear correlation coefficient, $r$, of the joint distributions.  In cases where the quantities are well correlated ($r \gtrsim 0.5$), the [s.d.]{}  of the difference can be significantly less than the [s.d.]{} ’s of the individual values. 

The general trend of increasing blueshifts with increasing degree of ionization is readily apparent from the table, which is arranged from highest to lowest ionization.  The mean and peak velocities always show a positive correlation between different lines.  The correlation is highest ($r \simeq 0.9$) between H$\alpha$ and [\[\]]{}, which is not surprising since there is a large overlap of the emission zones of these two lines.  However, even pairs of lines that should not greatly overlap (such as ${[\ion{N}{2}]}{} - {[\ion{O}{3}]}{}$ and ${[\ion{S}{2}]}{} - {[\ion{S}{3}]}{}$) also show significant correlation ($r > 0.5$). 

The [s.d.]{}  of 1.4 [$\mathrm{km\ s}^{-1}$]{} for ${\ensuremath{\langle V\rangle}}_{{\ensuremath{\mathrm{H}\alpha}}} - {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$ is a hard upper limit on the uncorrelated average wavelength calibration errors within the individual lines.  Given that ${\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}} -
{\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$ has a much larger [s.d.]{} ( 2.8 [$\mathrm{km\ s}^{-1}$]{}) and that some fraction of the [$\mathrm{H}\alpha$]{} emission comes from the [\[\]]{} zone, it follows that at least some of the scatter in ${\ensuremath{\langle V\rangle}}_{{\ensuremath{\mathrm{H}\alpha}}} -
{\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$ must be real.  This is consistent with our estimate above in § \[sec:atlas\] that the relative velocity calibration within each emission line is good to $\simeq
1~{\ensuremath{\mathrm{km\ s}^{-1}}}{}$. Our worst-case estimate for the uncertainty in the *absolute* velocity calibration for each line is 2 [$\mathrm{km\ s}^{-1}$]{} and this translates into uncertainties in the average values of the differences in mean and peak velocities.  This issue is returned to in § \[sec:deriv-mean-electr\] below, where we argue that the absolute velocity calibration must be better than this conservative estimate. 

An extremely tight correlation ($r \simeq 0.9$) is seen between the mean square velocity widths of H$\alpha$ and [\[\]]{}, whereas no other pair of lines show a significant correlation in this quantity.  The large offset $\bigl\langle \sigma^2_{{\ensuremath{\mathrm{H}\alpha}}} - \sigma^2_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}
\bigr\rangle = (82.2 \pm 18.2)~{\ensuremath{\mathrm{km^2\ s}^{-2}}}$ is due mainly to the extra thermal broadening of [$\mathrm{H}\alpha$]{}, as is discussed further in § \[sec:deriv-mean-electr\] below.  Other pairs of lines show much smaller offsets in $\sigma^2$, with no obvious trends with ionization.  The [s.d.]{}  of $\sigma^2$ is very large in the case of [\[\]]{} and [\[\]]{} because of the low signal-noise ratio of these weak lines. 

A similar picture is seen for the corrected full width half maxima, $W$: only in the case of ${\ensuremath{\mathrm{H}\alpha}}{} - {[\ion{O}{3}]}{}$ does a correlation exist.  Note that even though the effects of thermal and fine structure broadening have been subtracted in quadrature to obtain $W$, a residual difference remains between [$\mathrm{H}\alpha$]{} and [\[\]]{} of $\simeq
4~{\ensuremath{\mathrm{km\ s}^{-1}}}$. We believe that this is because the FWHM is such a mathematically ill-behaved quantity that quadrature subtraction cannot be expected to work reliably in the presence of strongly non-Gaussian lineshapes.  For example, consider a line that consists of two components, A and B, separated by $\delta v$ and each of width $w$, where the intensity of B is less than half that of A. The measured FWHM will not be affected at all by the presence of component B unless the components are blended: $\delta v \lesssim w$. The larger thermal broadening of [$\mathrm{H}\alpha$]{} gives a larger $w$, so that there will be cases where the FWHM of [$\mathrm{H}\alpha$]{} can “see” component B, but the FWHM of [\[\]]{} cannot.  Evidence for this can be seen in Figure \[fig:moments-fwhm\], where the [\[\]]{} map shows multiple dark patches of low FWHM, which are absent from the [$\mathrm{H}\alpha$]{} map. 

All the metal lines show very similar mean values of $\langle W\rangle
\simeq 15~{\ensuremath{\mathrm{km\ s}^{-1}}}{}$, except for [\[\]]{}, which is significantly lower, $\langle W\rangle \simeq 12~{\ensuremath{\mathrm{km\ s}^{-1}}}{}$. The FWHM is affected by noise to a much smaller extent than $\sigma^2$ and the metal lines all have low thermal widths, so should not be affected by the problem discussed in the previous paragraph. 

Derivation of the mean electron temperature from linewidths {#sec:deriv-mean-electr}
===========================================================

In this section, we use the statistics of the observed line profiles of [\[\]]{}, [$\mathrm{H}\alpha$]{}, and [\[\]]{} (§ \[sec:velocity-statistics\]) to derive the mean electron temperature in the nebula and to show that the “non-thermal” component to the line broadening is not significantly larger in [$\mathrm{H}\alpha$]{} than in the metal lines.  The method we employ is essentially an old one , which estimates the gas temperature by comparing the width of a metal line with that of a hydrogen line, under the assumption that the excess width of the hydrogen line is principally due to its greater thermal broadening on account of its lower atomic weight.  A problem with many applications of this method has been that the metal line and hydrogen line (usually [\[\]]{} and [$\mathrm{H}\alpha$]{}) are not necessarily emitted by the same volume of gas.  Indeed, previous applications of this method to the Orion nebula have shown that the [$\mathrm{H}\beta$]{}-[\[\]]{} pair gives results that are very different from the [$\mathrm{H}\alpha$]{}-[\[\]]{} pair and that are seemingly more reliable.  Our own version of the method improves on previous versions by using all three emission lines [\[\]]{}, [$\mathrm{H}\alpha$]{}, and [\[\]]{} in order to control for the imperfect overlap of the emission regions. 

Components of the observed line widths {#sec:comp-observ-line}
--------------------------------------

For each emission line, the observed mean square velocity width at a certain point on the nebula can be broken down schematically into four components: $$\label{eq:sigsq}
  \sigma^2 = \sigma^2_\mathrm{th} + \sigma^2_\mathrm{fs}
  + \sigma^2_\mathrm{ins} + \sigma^2_\mathrm{nt} .$$ The first term in this equation is the thermal Doppler width, $\sigma^2_\mathrm{th} = 82.5 (T_4/A)$ [$\mathrm{km^2\ s}^{-2}$]{}, where $T_4 =
T/10^4~\mathrm{K}$ and $A$ is the atomic weight of the atom or ion.  The second term is the fine structure broadening, which is important for hydrogen and helium recombination lines, but not for metal lines such as [\[\]]{} and [\[\]]{}.  In Appendix \[sec:fine-struct-comp\], we show that its value for [$\mathrm{H}\alpha$]{} is $\sigma^2_\mathrm{fs} = 10.233~{\ensuremath{\mathrm{km^2\ s}^{-2}}}$. The third term is the instrumental width, which for the KPNO observations is approximately $
\sigma^2_\mathrm{ins} = 11.5~{\ensuremath{\mathrm{km^2\ s}^{-2}}}$, and which is the same for the three lines.  The fourth term, $\sigma^2_\mathrm{nt}$, is the “non-thermal” contribution, which is a catch-all term for any other broadening process.  This might include Doppler broadening due to velocity gradients or turbulence within the emitting regions, broadening due to scattering by dust particles, or any other process. 

Applying the above equation to differences in $\sigma^2$ measured in § \[sec:velocity-statistics\] gives $$\begin{gathered}
  \label{eq:sigsqdiff1}
  \sigma^2_{{\ensuremath{\mathrm{H}\alpha}}} - \sigma^2_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}} = 77.34 T_4 + 10.233\\ 
  + \sigma^2_\mathrm{nt,{\ensuremath{\mathrm{H}\alpha}}} - \sigma^2_\mathrm{nt,{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}\end{gathered}$$ and $$\label{eq:sigsqdiff2}
  \sigma^2_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}} - \sigma^2_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}} = 0.74 T_4 
  + \sigma^2_\mathrm{nt,{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}} - \sigma^2_\mathrm{nt,{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}} .$$ Note that the instrumental contribution cancels from these equations.  The simplest way of extracting the temperature, and what has been done in the past, is then to assume that the difference in non-thermal widths is negligible, so that $T_4$ can be calculated directly from the measured $\sigma^2_{{\ensuremath{\mathrm{H}\alpha}}} - \sigma^2_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$ as $$\label{eq:Tnaive}
  T'_4 =  \bigl(\sigma^2_{{\ensuremath{\mathrm{H}\alpha}}} - \sigma^2_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}} - 10.233\bigr)/77.34$$ However, by including the data on the [\[\]]{} line, it is possible to derive a lower bound on the value of $\sigma^2_\mathrm{nt,{\ensuremath{\mathrm{H}\alpha}}} -
\sigma^2_\mathrm{nt,{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$ and so arrive at a more reliable temperature determination than the naive estimate $T'_4$.

Correction for differences in temperature and non-thermal widths between high and low ionization zones {#sec:corr-diff-non}
------------------------------------------------------------------------------------------------------

We adopt a simplified two-zone model for the emission structure along a given line of sight through the nebula.  The first zone, A, contributes all the [\[\]]{} emission plus a fraction, $f$, of the [$\mathrm{H}\alpha$]{} emission.  The second zone, B, contributes all the [\[\]]{} emission plus the remaining fraction, $1-f$, of the [$\mathrm{H}\alpha$]{} emission.  We further assume that the non-thermal broadening is the same for all emission lines that arise from the same volume (this assumption should hold true for turbulent or other kinematic broadening, an issue we will return to later).  A fundamental property of the moments of the line profile given in equation (\[eq:moments\]) is that they are all *additive*, with consequences that are elucidated in Appendix \[sec:prop-veloc-moments\].  From equation  it follows that $$\begin{gathered}
  \label{eq:twozone}
  \sigma^2_\mathrm{nt,{\ensuremath{\mathrm{H}\alpha}}} - \sigma^2_\mathrm{nt,{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}} = 
  (1 - f) (\sigma^2_\mathrm{nt,{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}} -
  \sigma^2_\mathrm{nt,{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}})\\
  + f (1 - f) \bigl( {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}} - {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}
  \bigr)^2 .\end{gathered}$$ Note that the above equation is only valid if the effects of dust extinction are the same for zones A and B, but this is probably a good approximation since @2000AJ....120..382O showed that the majority of the dust extinction in Orion arises in the neutral veil, in the foreground of the ionized gas. 

We also allow the two emission zones to have different temperatures, $T_{\mathrm{A}}$ and $T_{\mathrm{B}}$, such that the mean [$\mathrm{H}\alpha$]{} temperature is $T = f T_{\mathrm{A}} + (1-f) T_{\mathrm{B}}$. We can now combine equation (\[eq:twozone\]) with equations (\[eq:sigsqdiff1\]) and (\[eq:sigsqdiff2\]), modified to account for the temperature difference between the zones, to obtain $$\label{eq:t4}
  T_4 = 
  \eta T'_4 + \epsilon - \zeta - \xi ,$$ where $$\begin{gathered}
  \label{eq:eta}
  \eta = \bigl(1 - 0.00957\, (1 - f)\bigr)^{-1}, \\
  \label{eq:epsilon}
  \epsilon = 0.00957\, \eta\, f (1 - f) \bigl( T_{\mathrm{B}} - T_{\mathrm{A}} \bigr)/
  10^4~\mathrm{K}, \\
  \label{eq:zeta}
  \zeta = 0.01293\, \eta\, (1 - f) \bigl(\sigma^2_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}} - \sigma^2_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}} \bigr), \\
  \label{eq:xi}
  \xi = 0.01293\, \eta\, f (1 - f) \bigl({\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}} -
    {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}\bigr)^2 .  \end{gathered}$$

Techniques for estimating the ionization factor  {#sec:techn-estim-fact}
------------------------------------------------

Equation (\[eq:t4\]) gives the temperature at any point in the nebula in terms of observable quantities and the factor $f$, which can be estimated in two ways: from surface brightness ratios, or from ratios of differences in mean velocity.  We shall describe each of these techniques in turn. 

We define for each emission line a per-zone mean effective emission coefficient, ${\ensuremath{\beta}}$, with units $\mathrm{erg\ cm^3\ s^{-1}\
  sr^{-1}}$, such that the intrinsic surface brightness of a zone is $S = {\ensuremath{\beta}}{\ensuremath{E}}$, where ${\ensuremath{E}}= \int n_\mathrm{p} n_\mathrm{e}\,dz$ is the emission measure of the zone.  Since [[[\[\]]{}6584Å]{}]{} and [[[$\mathrm{H}\alpha$]{}6563Å]{}]{} are so close in wavelength, they are both affected in the same way by dust extinction, so the observed ratio of their surface brightnesses is equal to the intrinsic value.  Therefore, given that [\[\]]{} by construction arises only in Zone B, we have $$\begin{gathered}
  \label{eq:sratio}
  R_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}} \equiv \frac{S({[\ion{N}{2}]})}{S({\ensuremath{\mathrm{H}\alpha}})} = \frac{ {\ensuremath{\beta}}{\ensuremath{_{\mathrm{B}}}}({[\ion{N}{2}]})\,
    {\ensuremath{E}}{\ensuremath{_{\mathrm{B}}}}}{{\ensuremath{\beta}}{\ensuremath{_{\mathrm{A}}}}({\ensuremath{\mathrm{H}\alpha}})\, 
    {\ensuremath{E}}{\ensuremath{_{\mathrm{A}}}}+ {\ensuremath{\beta}}{\ensuremath{_{\mathrm{B}}}}({\ensuremath{\mathrm{H}\alpha}})\, {\ensuremath{E}}{\ensuremath{_{\mathrm{B}}}}} \\ = (1 - f) \frac{{\ensuremath{\beta}}{\ensuremath{_{\mathrm{B}}}}({[\ion{N}{2}]})
  }{{\ensuremath{\beta}}{\ensuremath{_{\mathrm{B}}}}({\ensuremath{\mathrm{H}\alpha}})}, \end{gathered}$$ where the last equality follows from the definition of $f$. Hence, if the emission coefficients ${\ensuremath{\beta}}{\ensuremath{_{\mathrm{B}}}}({[\ion{N}{2}]})$ and ${\ensuremath{\beta}}{\ensuremath{_{\mathrm{B}}}}({\ensuremath{\mathrm{H}\alpha}})$ are known, then $f$ can be calculated from the observed surface brightness ratio as $f = 1 - R_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}} {\ensuremath{\beta}}{\ensuremath{_{\mathrm{B}}}}({\ensuremath{\mathrm{H}\alpha}}) / {\ensuremath{\beta}}{\ensuremath{_{\mathrm{B}}}}({[\ion{N}{2}]})$. To make use of this equation, it is first necessary to photometrically calibrate the surface brightness maps and we have done this by comparison with published spectrophotometry of the region [@1992ApJ...399..147P; @2000ApJS..129..229B; @2003AJ....125.2590O].  It is then necessary to estimate the ratio of emission coefficients for conditions typical of Zone B: $T \simeq
10^4$ K, $n_{\mathrm{e}} \simeq 5000$, which we have done using the Cloudy plasma code [@2000RMxAC...9..153F], assuming an N/H abundance ratio of $6\times 10^{-5}$, finding that the ratio is roughly unity.  This is consistent with the N[$^+$]{}/N abundance found by @1998MNRAS.295..401E.

![ image](tmaps-figure)

Assuming that ${\ensuremath{\beta}}{\ensuremath{_{\mathrm{B}}}}({\ensuremath{\mathrm{H}\alpha}}) / {\ensuremath{\beta}}{\ensuremath{_{\mathrm{B}}}}({[\ion{N}{2}]}) = 1$, we find a mean value of $\langle f \rangle = 0.76$ for the entire area covered by our observations.  The upper-left panel of Figure \[fig:tmaps\] shows the map of $f$ calculated by this technique, which varies between 0.26 and 0.91.  In reality, ${\ensuremath{\beta}}{\ensuremath{_{\mathrm{B}}}}({\ensuremath{\mathrm{H}\alpha}}) / {\ensuremath{\beta}}{\ensuremath{_{\mathrm{B}}}}({[\ion{N}{2}]})$ will vary with position due to variations in electron density ([\[\]]{} is affected by collisional deexcitation but [$\mathrm{H}\alpha$]{} is not).  A better way of calculating $f$ would be to use the ${[\ion{O}{3}]}/{\ensuremath{\mathrm{H}\beta}}$ surface brightness ratio, since the critical density for [\[\]]{} is higher than that typically found in the nebula (similar arguments to those above give $f = R_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}} {\ensuremath{\beta}}{\ensuremath{_{\mathrm{A}}}}({\ensuremath{\mathrm{H}\beta}}) / {\ensuremath{\beta}}{\ensuremath{_{\mathrm{A}}}}({[\ion{O}{3}]})$, where $R_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}} = S({[\ion{O}{3}]})/S({\ensuremath{\mathrm{H}\beta}})$).  However, our dataset does not include [$\mathrm{H}\beta$]{} and the ${[\ion{O}{3}]}/{\ensuremath{\mathrm{H}\alpha}}$ ratio would have to be corrected for the highly variable foreground extinction [@2000AJ....120..382O].  A comparison of the ${[\ion{N}{2}]}/{\ensuremath{\mathrm{H}\alpha}}$ and ${[\ion{O}{3}]}/{\ensuremath{\mathrm{H}\beta}}$ ratios over a limited area (using data from Figs.  11 and 12 of [@2003AJ....125.2590O]) gives an rms scatter of $\pm0.05$ between the values of $f$ derived from the $R_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}}$ and $R_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$ techniques, which is roughly half of the total rms variation in $f$ from our maps.  The $R_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$-derived values of $f$ have a greater range than the $R_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}}$-derived values.  This is to be expected since, all else being equal, higher density regions will have lower ionization. 

![ Joint distribution of differences in mean velocities: ${\ensuremath{\langle V\rangle}}_{{\ensuremath{\mathrm{H}\alpha}}} - {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$ versus ${\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}} -
    {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$. The values of ${\ensuremath{\langle V\rangle}}_{{\ensuremath{\mathrm{H}\alpha}}} -
    {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$ have been shifted by 0.5 [$\mathrm{km\ s}^{-1}$]{} so as to make the distribution symmetric about a line through the origin (see discusion in text).  The dashed lines show the boundaries of the two triangular regions (upper right and lower left) that give physically possible values of $f$ in the two-zone model.  The solid line shows the relationship ${\ensuremath{\langle V\rangle}}_{{\ensuremath{\mathrm{H}\alpha}}} - {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}} =
    (1 - \langle f\rangle)({\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}} - {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}})$, where $\langle f\rangle$ is determined from the $R_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}}$ method.  The solid gray circle indicates the effect of $\pm 1~{\ensuremath{\mathrm{km\ s}^{-1}}}$ uncertainties in the velocities. []{ data-label="fig:vmean-diff"}](vmean-diffs)

A completely independent way of estimating $f$, and one that does not depend on knowing the emission coefficients, is to use the observed mean velocities of the three lines.  A corollary of the additive property of the line profile velocity moments (Appendix \[sec:prop-veloc-moments\]) is that in the two-zone model we have from equation  that $$\label{eq:addvmean}
  {\ensuremath{\langle V\rangle}}_{{\ensuremath{\mathrm{H}\alpha}}} = f {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}} + (1 - f) {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}},$$ so that $f$ can be estimated as $$\label{eq:fvmean}
  f = 1 - \frac{{\ensuremath{\langle V\rangle}}_{{\ensuremath{\mathrm{H}\alpha}}} - {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}}
  {{\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}} - {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}}.$$ Despite its advantages over the surface brightness method, this technique has two serious limitations.  First, equation (\[eq:fvmean\]) has a singularity when ${\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}}
= {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$. If the two-zone model were strictly valid, then any pixel with ${\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}} = {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$ would also have ${\ensuremath{\langle V\rangle}}_{{\ensuremath{\mathrm{H}\alpha}}} = {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$, so the singularity would be benign.  However, in reality the two-zone model is only an approximation and, additionally, the observed values of [$\langle V\rangle$]{} will always include a certain measurement error.  As a result, applying equation (\[eq:fvmean\]) to real data will inevitably result in discontinuities in the $f$ map and unphysical values ($f <0$ or $f >
1$) for some pixels.  Second, the observed values of ${\ensuremath{\langle V\rangle}}_{{\ensuremath{\mathrm{H}\alpha}}}
-{\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$ are rather small (typically $-2$ to $+3$ [$\mathrm{km\ s}^{-1}$]{}, see Fig.  \[fig:meanvsmean\]*a*), and are of the same order as the systematic uncertainties in the wavelength calibration. 

Given these caveats, it is perhaps surprising that this method works at all.  However, the average *slope* of the relationship between ${\ensuremath{\langle V\rangle}}_{{\ensuremath{\mathrm{H}\alpha}}} - {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$ and ${\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}} -
{\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$ should be independent of any inaccuracies in the absolute velocity scales of the individual lines, and by equation (\[eq:fvmean\]) should be equal to $1 - \langle
f\rangle$. This is illustrated in Figure \[fig:vmean-diff\], where it can be seen that the value $\langle f\rangle = 0.76$, calculated from the $R_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}}$ method above and shown by the solid line, is indeed close to the mean trend in the velocity differences.  Furthermore, to the extent that the two-zone approximation is valid, the symmetry axis of the joint distribution of velocity differences ought to pass through the origin.  This is only true when we shift all our ${\ensuremath{\langle V\rangle}}_{{\ensuremath{\mathrm{H}\alpha}}} - {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$ values by $+0.5~{\ensuremath{\mathrm{km\ s}^{-1}}}{}$ and this is what is done in the figure.  Such a small shift is well within the estimated uncertainties in the velocity zero point of each line.  Of course, there may also be a similarly small constant error in the ${\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}} - {\ensuremath{\langle V\rangle}}_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$ values, but we have no information that would allow us to determine the correction in this case. 

In order to correspond to a value of $f$ between 0 and 1, the velocity differences must lie in one of the two triangular areas of the graph in Figure \[fig:vmean-diff\], bounded by dashed lines.  It can be seen from the graph that the majority of pixels do indeed lie in these areas, mostly on the positive side of the origin.  The small fraction of pixels lying in unphysical areas of the graph can be explained by residual errors in the velocity determination of order $1~{\ensuremath{\mathrm{km\ s}^{-1}}}$ (indicated by the solid gray disk).  The lower left panel of Figure \[fig:tmaps\] shows the map of $f$ calculated using this velocity difference technique.  For much of the map, the agreement with the surface brightness ratio technique (upper left panel) is rather good, and the low ionization regions associated with the linear bright bars are readily apparent in both maps.  The unphysical values of $f$ show up as solid black or white regions in the map.  The fact that these regions occur predominantly at the edges of the map, with white regions ($f >1$) concentrated in the West and black regions ($f < 0$) in the East, suggests that the residual velocity errors are largely systematic, rather than random. 

Corrected value of the electron temperature {#sec:corr-value-electr}
-------------------------------------------

                    $\langle T_4 \rangle$   $\langle T_4 \rangle{_{\scriptscriptstyle\!S}}$   $t^2_{\mathrm{A}}$   $N_{\mathrm{pix}}$
  --------------- ----------------------- ------------------------------------------------- -------------------- --------------------
            naive                   0.930                                             0.942                0.025               178432
        corrected                   0.913                                             0.916                0.014               178432
     low $\sigma$                   0.916                                         **0.919**            **0.012**               173180
    high $\sigma$                   0.824                                             0.792                0.104                 5252

  : Summary of temperature determinations[]{data-label="tab:temp"}

We are finally in a position to calculate the correction factors, $\epsilon$, $\zeta$, and $\xi$ needed to estimate the electron temperature.  The mean values are $\langle\epsilon\rangle = 2.9 \times
10^{-4}$, $\langle\zeta\rangle = -0.0218$ and $\langle\xi\rangle =
0.0417$, so that the corrected mean temperature from equation (\[eq:t4\]) is $\langle T_4 \rangle = 0.913$, whereas the naive estimate would have been $\langle T'_4 \rangle = 0.930$ (note that all these averages are calculated with equal weighting for each pixel).  Maps of the spatial variation of the correction factors $\zeta$ and $\xi$ are shown in the central panels of Figure \[fig:tmaps\], while the right-hand panels show maps of the naive and corrected temperatures.  The correction factor $\epsilon$ is uniformly very small and will not be considered further.  The effect of the corrections on the mean temperature is rather small, merely lowering it by 200 K, but, as can be seen from the maps, the correction factors can be much larger in localized regions ($>
2000~\mathrm{K}$).  Some of the variation seen in the uncorrected $T'_4$ map is removed by the correction factors, with the most notable example being at the bright bar, which is seen as a dark diagonal band in the $T'_4$ map, but not at all in the $T_4$ map, principally due to $\zeta$. Another example is the “Red Fan” [@2007AJ....133..952G], seen as a light feature in the lower right of the $T'_4$ map, but which is effectively removed by the $\xi$ correction (the thin vertical streaks seen superposed on this feature are observational artifacts, which can also be seen in the $\sigma_{{\ensuremath{\mathrm{H}\alpha}}}$ map of Fig.  \[fig:moments-sigma\]).  Other striking features seen in the temperature maps are those associated with the high-velocity jets emanating from the Orion South region [@2007AJ....133.2192H]: [HH 202]{} (northeast), [HH 203/204]{} (southwest), and [HH 529/269]{} (east-west).  These features also tend to have high absolute values of $\zeta$ or $\xi$, but the application of the correction factors is not successful in removing them from the temperature map.  Despite this, we do not believe that they correspond to real temperature variations.  Since all these features have large rms velocity widths (Fig.  \[fig:moments-sigma\]), we have masked out of the corrected map all pixels with $\sigma^2_{\mathrm{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}}$ or $\sigma^2_{\mathrm{{\mathrm{[{\setcounter{ionstage}{2}  \ensuremath{\mathrm{N\,\Roman{ionstage}}}}]}}}}$ greater than 180 [$\mathrm{km^2\ s}^{-2}$]{} (visible as small solid gray areas in the figure). 

In order to quantify the variation in the naive and corrected temperature maps, we calculate the plane-of-sky Peimbert $t^2$ parameter [@1967ApJ...150..825P; @2003AJ....125.2590O]: $t_{\mathrm{A}}^2 = \bigl\langle (T - \langle T\rangle{_{\scriptscriptstyle\!E}})^2
\bigr\rangle{_{\scriptscriptstyle\!E}}/ \langle T\rangle{_{\scriptscriptstyle\!E}}$. The averages in this expression, $\langle\cdots\rangle{_{\scriptscriptstyle\!E}}$, should theoretically be weighted by the emission measure of each pixel, but we approximate this by weighting by the [$\mathrm{H}\alpha$]{} surface brightness: $\langle\cdots\rangle{_{\scriptscriptstyle\!S}}$. The results are summarised in Table \[tab:temp\].  It can be seen that the [$\mathrm{H}\alpha$]{}-weighted averages, $\langle T_4 \rangle{_{\scriptscriptstyle\!S}}$, are hardly any different from the uniform-weighted averages, $\langle T_4 \rangle$. The fractional variance, $t^2_{\mathrm{A}}$, is significantly reduced after applying the temperature corrections (second row of table).  A further small reduction in $t^2_{\mathrm{A}}$ occurs (third table row) when the high-$\sigma$ pixels are eliminated as discussed above.  The high-$\sigma$ pixels themselves (fourth table row) show a significantly lower mean temperature and a very large $t^2_{\mathrm{A}}$.

Discussion {#sec:discussion}
==========

Comparison with previous mean temperature determinations {#sec:comp-with-prev}
--------------------------------------------------------

Previous determinations of the mean value of the electron temperature in the Orion Nebula have been made using various methods (see [@2006agna.book.....O], Chapter 5 for an overview).  In this section, we compare our own results (§ \[sec:corr-value-electr\]) with some recent examples of these. 

@2003AJ....125.2590O used *HST* emission line imaging in the [\[\]]{} 4363 Å and 5007 Å lines to determine $T_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}} = 8890 \pm 496$ K, where the scatter represents real variations, rather than observational uncertainties.  Before comparing with our own result, this temperature needs to be corrected in two ways.  First, temperature fluctuations *within* the [\[\]]{}-emitting region will lead to this method overestimating the real temperature there, $T(\mathrm{O}^{++})$. Second, $T(\mathrm{O}^{++})$ (which corresponds to $T_\mathrm{A}$ from § \[sec:corr-diff-non\] above) is an underestimate of the mean temperature of the ionized gas, $T(\mathrm{H}^{+})$, since it does not include the low ionization zone along each line of sight (Zone B), which tends to have a higher temperature.  Taking estimates of the magnitudes of these corrections from @1998MNRAS.295..401E and @2003MNRAS.340..362R, we find that the @2003AJ....125.2590O results imply $T(\mathrm{O}^{++}) =
8400$ K and $T(\mathrm{H}^{+}) = 9100$ K, which is very similar to our result of 9190 K (Table \[tab:temp\]). 

Electron temperatures have also been derived from the brightness relative to the continuum of radio recombination lines , consistently finding values of around $8500$ K. Yet another technique is to measure the strength of the Balmer discontinuity at 3646 Å, which was used by @1995ApJ...450L..59L to measure a mean value of $T \simeq 9000$ K, whereas @1998MNRAS.295..401E used the same technique to find $T \simeq 8600$ K. In the presence of temperature fluctuations, these values are underestimates and should be increased by $\simeq 300$ K. However, the experimental uncertainties of these results are probably higher than for the other methods and furthermore they only correspond to restricted regions of the nebula. 

In summary, once temperature fluctuations are accounted for, there is a very good agreement between the three independent optical determinations of the mean electron temperature (line ratios, Balmer discontinuity, and line widths), giving an average value of $\simeq
9100 \pm 100$ K. The radio-derived temperatures are significantly lower by about $600$ K and we have no satisfactory explanation for this discrepancy. 

Plane-of-sky variation in the electron temperature {#sec:plane-sky-variation}
--------------------------------------------------

The point-to-point variations in electron temperature across the face of the nebula have been studied in detail by two recent papers [@2003AJ....125.2590O; @2003MNRAS.340..362R] based on optical line ratios from *HST* observations.  The values obtained for $t^2_{\mathrm{A}}(\mathrm{O{\ensuremath{^{++}}}})$ for different regions of the nebula range from 0.005 to 0.016.  @2003MNRAS.340..362R also calculate $t^2_{\mathrm{A}}(\mathrm{N{\ensuremath{^{++}}}})$, finding broadly similar values.  After correcting for the effects of noise, @2003AJ....125.2590O find a value for their entire surveyed region (roughly one third the area of our own dataset) of $t^2_{\mathrm{A}}(\mathrm{O{\ensuremath{^{++}}}}) = 0.0079$. Given that our own value of $t^2_{\mathrm{A}} = 0.012 $ will also be somewhat affected by noise, the agreement may at first seem satisfactory.  However, the spatial resolution of our own observations is at least an order of magnitude lower than these two studies, which find that power spectrum of the fluctuations peaks at sub-arcsecond scales.  Such fluctuations would be entirely smoothed out by atmospheric seeing in our own dataset, so one would expect our results to show a much lower value of $t^2_{\mathrm{A}}$.

An analysis of the temperature variations along the slits used in the @2003MNRAS.340..362R study supports this point of view.  Their Slit 5 is the only one that does not cross high-velocity jet features, and we find $t^2_{\mathrm{A}} = 0.002$ for this slit, which is very much lower than their value of 0.018, as is to be expected from our lower spatial resolution.  It therefore seems likely that much of the variation visible in even our corrected temperature map (lower right panel of Fig.  \[fig:tmaps\]) is spurious, and that $t^2_{\mathrm{A}} = 0.002$ is a more accurate estimate of the plane-of-sky temperature fluctuations at scales $> 1$ arcsec.  This is also consistent with calculations of the variation in three small $40\arcsec \times 40\arcsec$ patches of our map, which were selected for their high signal-noise ratio and avoidance of high-velocity flows.  The resulting values are $t^2_{\mathrm{A}} = 0.002$, 0.004, and 0.002 for regions centered at $(-30,+60)$, $(-30,-100)$, and $(-80,-50)$, respectively. 

The nature of the non-thermal line broadening {#sec:the-nature-non}
---------------------------------------------

A fundamental assumption of our method of determing the temperature is that the non-thermal broadening depends only on the emitting volume (see § \[sec:corr-diff-non\]), and does not vary systematically between collisionally excited and recombination lines.  This is necessary so that $\sigma^2_{\mathrm{nt,{\ensuremath{\mathrm{H}\alpha}}}}$ can be determined from the kinematics of the [\[\]]{} and [\[\]]{} lines, as in equation (\[eq:twozone\]).  The fact that we determine a mean nebular temperature that closely agrees with that from other techniques (§ \[sec:comp-with-prev\]) is a strong argument that this assumption is valid. 

However, this is in conflict with the result of @2003AJ....125.2590O, who found non-thermal widths (FWHM) of $\simeq 19~{\ensuremath{\mathrm{km\ s}^{-1}}}$ for recombination lines ([$\mathrm{H}\alpha$]{}, [$\mathrm{H}\beta$]{}, and He I 5876 Å), but only $\simeq 11~{\ensuremath{\mathrm{km\ s}^{-1}}}$ for collisional lines ([\[\]]{}, [\[\]]{}, [\[\]]{}, [\[\]]{}, [\[\]]{}).  Taken at face value, these widths would imply a large difference of $(19^2 - 11^2) / (8 \ln 2) =
43~{\ensuremath{\mathrm{km^2\ s}^{-2}}}$ between $\sigma^2_{\mathrm{nt,{\ensuremath{\mathrm{H}\alpha}}}}$ and $\sigma^2_{\mathrm{nt,{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}}$, meaning that over half of the difference $\sigma^2_{{\ensuremath{\mathrm{H}\alpha}}} - \sigma^2_{{\mathrm{[{\setcounter{ionstage}{3}  \ensuremath{\mathrm{O\,\Roman{ionstage}}}}]}}}$ (Fig.  \[fig:sigmavssigma\]*a*) is due to non-thermal effects.  In this case, the mean electron temperature would only be $\sim
4000$ K, in stark contrast to all other temperature determinations for the nebula. 

Instead, we suggest that the difference in non-thermal widths found by @2003AJ....125.2590O is an artifact of their methodology, which consisted in fitting Gaussian components to the line shapes.  The basic problem with this is that the extra thermal broadening of the H and He lines tends to blend the fine details of the line shape, allowing a satisfactory fit to be obtained using fewer Gaussian components than are necessary for the metal lines.  These fewer components would then naturally tend to have larger individual widths (even after correcting for thermal broadening), which may explain the @2003AJ....125.2590O result.  Our own method of calculating the RMS linewidths does not suffer from this problem, since it avoids any subjective decision about how many components to fit.  Indeed, an application of our methodology to the *same data* as used by @2003AJ....125.2590O shows no evidence for extra non-thermal broadening of the recombination lines (Fig.  4 of [@1999AJ....118.2350H]).  A further problem with the @2003AJ....125.2590O study is that no correction was made for the H fine-structure broadening (although the larger He I fine structure was accounted for). 

A third methodology was employed by @1976MNRAS.174..105G, who convolved the observed [\[\]]{} profile with Gaussians of different widths and determined which width gave the best match to the observed [$\mathrm{H}\beta$]{} profile.  This technique gave temperatures that are very similar to those derived here, which is further evidence that the @2003AJ....125.2590O widths cannot be correct. 

One promising class of explanations for the non-thermal linewidths is that they are kinematic in nature, caused by velocity structure in the ionized gas.  This may be due to velocity gradients in ordered, large-scale champagne flows , or to multi-scale, turbulent motions [e.g. @2006ApJ...647..397M; @2007dmsf.book..103H].  However, as @2005ApJ...627..813H showed in their § 6.1.2, velocity gradients in an ordered flow are incapable of explaining the observed linewidths.  Furthermore, the disordered appearance of the maps of the line profile parameters (Figs.  \[fig:moments-mean\] to \[fig:moments-fwhm\]) would favor the turbulent explanation.  In the radiation-hydrodynamical simulations of @2006ApJ...647..397M, the turbulence in the ionized gas is driven primarily by photoevaporation flows from dense concentrations of molecular gas, which are themselves the product of supersonic turbulence in the molecular cloud from which the high-mass star formed.  The internal velocity dispersion of the region remains roughly sonic for the full duration of the simulations ($\simeq
0.5$ Myr), in agreement with the observed linewidths in Orion.  In addition, it is quite possible that other mechanisms may be contributing to the turbulent driving, such as ionization front instabilities [@1996ApJ...469..171G; @1999MNRAS.310..789W; @2007astro.ph..3463W] or supersonic jets from young stars. 

Conclusions {#sec:conclusions}
===========

We have presented an atlas of spatially resolved high resolution emission line spectra of the Orion Nebula, covering a wide range of degrees of ionization.  This has allowed us to investigate the global kinematics of the nebular gas in unprecedented depth.  We have studied the statistical correlations between the line profile parameters of the different emission lines.  Our principal conclusions are as follows:

1.   The mean velocities of different emission lines are always positively correlated.  The correlation is greatest when there is considerable spatial overlap between the emitting volumes of the two lines, but correlation also exists between lines from non-overlapping volumes.  This is indicative of the existence of an ordered mean flow that connects the spatially separate regions. 

2.   There is no significant correlation, on the other hand, between the linewidths of emission lines from non-overlapping volumes, such as [[[\[\]]{}6584Å]{}]{} and [[[\[\]]{}5007Å]{}]{}.  This is consistent with kinematic broadening of an essentially turbulent nature.  The kinetic energies associated with both ordered and unordered motions are of the same order as the thermal energy. 

3.   The mean electron temperature in the nebula is determined to be 9190 K from comparison of the widths of [$\mathrm{H}\alpha$]{} and [\[\]]{} lines, after correction for differences in non-thermal widths between the different emission zones.  Although apparent spatial variations in temperature are detected ($t^2_{\mathrm{A}} = 0.012$), most are probably due to uncertainties in these corrections, rather than real temperature variations.  In regions of the nebula where these correction terms are small, the spatial variations are much smaller ($t^2_{\mathrm{A}} \simeq 0.002$) at the scales accessible to our observations ($> 1~\arcsec$). 

4.   We find that there can be no systematic difference in the non-thermal linewidths between collisionally excited and recombination lines.  This contradicts a recent claim that the recombination lines are significantly broader [@2003AJ....125.2590O]. 

We acknowledge financial support from DGAPA-UNAM projects PAPIIT IN112006 and IN110108, and IN116908.  MTGD is supported by a postdoctoral research grant from CONACyT, Mexico.  We are extremely grateful to Bob O’Dell for many discussions and advice, and for freely sharing observational data.  We thank Michael Richer, John Meaburn, Alex Raga, and Jorge García-Rojas for helpful discussions.  The referee is thanked for a useful report. 

Some properties of the velocity moments {#sec:prop-veloc-moments}
=======================================

The velocity moments $M_k$ of an emission line, as defined in equation , are all linear functions of the line intensity $I(V)$. This means that when two or more optically thin emission regions (A, B, C, …) are superimposed along the same line of sight then the moments of the combined emission can be found by simple addition of the moments of the individual regions: $$\begin{gathered}
  \label{eq:momsum}
  M_0 = M_{0,A} + M_{0,B} + M_{0,C} + \cdots , \\
  M_1 = M_{1,A} + M_{1,B} + M_{1,C} + \cdots , \\
  M_2 = M_{2,A} + M_{2,B} + M_{2,C} + \cdots .\end{gathered}$$

The first three formal moments $M_k$ are related to the more physical quantities: $$\begin{gathered}
  \label{eq:intensity-def}
  \text{total intensity: } I = M_0 ,\\
  \label{eq:vmean-def}
  \text{mean velocity: } {\ensuremath{\langle V\rangle}}= M_1/M_0 ,\\
  \label{eq:sigma-def}
  \text{rms velocity width: } \sigma = \bigl(M_2/M_0 - {\ensuremath{\langle V\rangle}}^2\bigr)^{1/2} .\end{gathered}$$ These too can be calculated for the superposition of multiple optically thin regions.  In particular, for the superposition of two regions, A and B, one can combine equations – to find $$\begin{gathered}
  \label{eq:intensitycombo}
  I = I_{A} + I_{B} ,\\
  \label{eq:meancombo}
  {\ensuremath{\langle V\rangle}}= \frac{{\ensuremath{\langle V\rangle}}_{A} I_{A} + {\ensuremath{\langle V\rangle}}_{B} I_{B}}{I_{A} + I_{B}} ,
  \\
  \label{eq:sigcombo}
  \sigma^2 = \frac{\sigma^2_{A} I_{A} + \sigma^2_{B} I_{B}}{I_{A} +
    I_{B}} + \frac{I_A I_B ({\ensuremath{\langle V\rangle}}_{A} - {\ensuremath{\langle V\rangle}}_{B})^2 }{(I_{A} +
    I_{B})^2} .  \end{gathered}$$

Fine structure components of Hydrogen Balmer lines {#sec:fine-struct-comp}
==================================================

The fine structure components of the Balmer lines have velocity differences of $< 7~{\ensuremath{\mathrm{km\ s}^{-1}}}{}$, which are too small to cause observable line splitting at typical nebula temperatures, where the Doppler broadening is of order $20~{\ensuremath{\mathrm{km\ s}^{-1}}}{}$. However, they do make a significant contribution to the observed line width and, furthermore, must be considered in order to make an accurate determination of the line-center wavelength.  We use the calculations of , for densities $n \le 10^4~{\ensuremath{\mathrm{cm}^{-3}}}{}$ and interpolated at a temperature of 9100 K. The resulting fine structure contribution to the velocity dispersion, $\sigma^2_\textrm{fs}$, and rest wavelengths, $\lambda_0$, are given in the first two columns of Table \[tab:fine\]. 

  ----------- -------- ----------- -- -------- ------
                                               
  line                                         
  H$\alpha$     10.233   6562.7910      $-0.7$   0.62
  H$\beta$       5.767   4861.3201         0.2   0.81
  H$\gamma$      4.635   4340.4590         1.3   0.68
  ----------- -------- ----------- -- -------- ------

Many previous studies have used outdated and incorrect values for the rest wavelengths, with a value of 6562.82 Å being commonly assumed for [$\mathrm{H}\alpha$]{}, rather than the correct value of 6562.7910 Å. This leads to systematic errors of order 2 [$\mathrm{km\ s}^{-1}$]{} in the derived radial velocities of the Balmer lines, as can be seen, for example, in Fig.  6 of @2000ApJS..129..229B. However, as is shown in Table \[tab:fine\], these discrepancies disappear when the correct rest wavelengths are used.